- Email: [email protected]
An extension of Lagrange's method to certain systems of first order partial differential equations
Jonmal of Atmorpheric
ad Termmtrld Phy&x. 1072,Vol. 84,pp. 160-164.%~'tPmon Prey.
SHORT PAPER An extensionof Lafmmge’s method to oewtain SJrstemsof ii& order ptutid difkalu equatiolla W. M. &TEN* and R. D. HAIWS Department of Electrical EngineeringUtah State University, U.S.A. and J. D.
WATSON
M&hem&ice Department Utah State University, U.S.A. (Received3 Mcry 1971) Ah&a&--A modificationof the method of characterieticehas been developedto solve systems of first orderhyperbolicpartial differentialequationawhere the characteristiccurves are essentially identical. This method ie applicable to fluid mixturea whose component fluid velocities are the same. Techniqueswere developed to find the solutionsat desiredpoints in space and time and to handle the conditionsat the boundaries. Results of the applicationof this method to the solution of ion density distributions in the Earth’s ionosphereare illustrated. 1. INTBODUCTION
of the interaction processes which occur in a fluid mixture require knowledge of the density distribution of each component of the mixture. These individual densities are given by the simultaneous solution of the set of continuity equations, one for each fluid component, which describe the component fluids. In general this system of equations is difficult to solve, especially if one or more of the fluid densities are changing rapidly. The purpose of this paper is to describe a numerical technique which is quite simple to use and has applicability to situations where the individual fluid components have the same or nearly the same velocities, and where the interaction and transport processes make diffusion an unimportant factor. The example cited in this paper is the calculation of ion densities in the E-region of the ionosphere (90-140 km). The primary fluids are the ion species of molecular oxygen O,+ and nitric oxide NO +. Each of the ion species are produced and destroyed by collisions with other atmospheric gases as well as by ionization from external radiation. In the E-region, wind transport processes are also an important factor in the ion density distribution. In order to theoretically model this ionospheric plasma, it is necessary to calculate the density distribution of each of the ion species under the influence of the reaction and transport processes. Since O,+ and NO+ have approximately the same mass, a modification of the method of characteristics has been developed to solve for the ion densities. THE STUDY
* Present address: Department of Physics, University of Prince Edward Island, Charlottetown, Prince Edward Island, Canada. 11
169
W. M. C!EBN, R. D. HARBIS end J. D. WATSON
160
The mathematical theory is considered in Section 2 and the results of actual calculations of the ionospheric ion density distribution are given in Section 3. 2. THEORY The numerical solution of the single continuity equation
qJ + v . (V(x, f) iv@,t))
= P(x, t, N)
(1)
with initial condition N(x, to) = h(x)
(2)
can be solved by Lagrange’s method of converting to a system of ordinary differential equations (COURANT and HILBEBT,1962) for which standard numerical programs are available to do the integration. To solve equation (1) by Lagra@s method we must solve the corresponding subsidiary equations dx - = V(x, t) dt
dN = P(x, t, N) - N(x, t)(V dt
l
V(x, t)).
The initial condition iV(x, to) = h(x) must be known at every point over the region of interest. The solution of (3) yields the so-called charaoterikic curve, which can be seen to be independent of N(x, t). Once the characteristic curve is known, then N(x, t) can be determined along it by solving (4). For a fluid mixture composed of ti components, a complete solution requires n continuity equations of the form
awx,t) + v at
l
(V,(x, wtx, f)) = -zp,(x, t, N,, NW- * * 9 N,)
i=l,2,...,n
and initial conditions N,(x, to) = M)
i = 1, 2, . . . , 12.
Extending Lagrange’s method to the system of continuity equations requires the solution of the corresponding subsidiary equations,
dx, - = dt
~*(x, t) dt
Vi(X, t)
= FAX,4 N,, iv*) - Nl(X, w
l
VAX,0)
i = 1,2,. . . ,n
(6)
where N,(x, to), i = 1, 2, . . . , n are initial conditions. In many applications the velocities of the fluid particles can be found in terms of the external forces, and are independent of the velocities and densities of the other species. When this is so the characteristic curves are independent. This means that throngh any point (z, t) it is neoessary to integrate (6) for each of the n speoies along its partienk ~&a&&tic
Certain eyetemf3of fht
order partial differentialequations
101
curve from 1, to time t = 7. Each of the component densities is evaluated on its own characteristic curve. Since the characteristic curves are different, the densities are evaluated at different spatial positions. To be correct, however, densities of interacting species must all be evaluated at the same position. This limitation of the method of characteristics can be eased somewhat by an iteration technique, but because of its difhculty it will have to be the subject of another paper. A numerical scheme to solve two simultaneous hyperbolic differential equations by Lagrange’s method has been reported by LISTER (1965). If the velocities of the individual fluid components of the mixture are nearly the same, then the chatacteristic curves are almost identical. For this condition one can then utilize a single characteristic curve and all the densities are found along this curve. The set of equations (6) reduces to the single equation dx - = V(x, t). dt
(7)
The shapes of the characteristic curves obviously depend upon the fluid velocity profile. For convergent or divergent flows the characteristic curves in general are not oonveniently oriented to allow density calculations over the complete region of interest. A technique to eliminate this problem allows the specification of points at which the final density values are desired. The solution technique then requires that the characteristic curves through each point of interest be integrated backwards in time to the initial starting time t,. But at t,, the densities iV(x, to) me known everywhere, so that a value of density is available to start the forward integration of (6). The densities are then found along the same characteristic curves, in particular at the positions that were initially specified. This method provides a check on the accuracy of the calculations because the backward integration started at the point (x, t) and the forward integration should return to the point (x, t). Very often it is impractical to model the entire fluid body because of limitations in size or complexity of the region. In these situations it is advantageous to choose a restricted region where the important reactions are taking place. But before numerical computations can be made, the conditions at the boundaries of this restricted region must be specified. If a characteristic curve either enters or leaves the region (i.e. passes through the boundary) at any time during the calculations this means physically that fluid is passing across the boundary, and in order to perform computations inside the region, the fluid flux must be known at the boundaries for all time. But this sort of information is not usually known. Thus it may be advantageous to choose or assume a closed system with no flux at the boundary. For certain cases this may also be a valid assumption if the flux is small or if the principal interactions of the fluid constituents take place inside the sub-region. The assumption that no fluid passes through the boundary implies that the characteristic curves do not pass through the boundaries. This can be implimented by allowing or forcing particle velocities to go to zero at the boundary. Zero velocities at the boundary, however, imply singular boundary conditions and instabilities in numerical calculations near the boundary. This can be seen by the following illustration in one space dimension: Assume z1 is a boundary point and that the velocity u(z) vanishes at z1 but that
W. M. CuEN,R. D. HAERISand J. D. WATSON
162 v(z)
2
0 and dv/dx > 0 near xi, but inside the region. The characteristic equation dx z = v(x)
which passes through point xg near xi at time t, can be approximated by the first two terms of a Taylor’s series dx dv z = v(x1) + z ;t_ 21(x - x1).
since2)(x1) =
0, the approximate velocity is given by the linear, constant-ooeflicient
equation
dx --dt
dv dx =rZ1
(x - Xi) = 0
whioh is valid for x near xi. This has solution x-x2,=Cexp
(t - t1) z = .%I
(8)
where C is a suitable constant. If we go backward in time from t,, x leaves x2 and approaohes x1 asymptotically but never reaches it. Thus the characteristic curves tend to bunch near the boundaries. Forward integration is then unstable in the sense that a small error in the position of the characteristic or the calculated velocity can cause large errors in the end positions of the charaoteristic curves. This is shown schematically in Fig. 1 where we have shown how the oharacteristics can bunch at a boundary where V(xi, t) = 0 but (dV/dx)I, _ 21# 0 at the boundary. One way to increase the stability of calculations near the boundaries is to increase the size of the region of interest so that the fluid velocity can be reduced gradually (dV/dx small) to zero at the boundaries. 3. APPLICATION The numerical calculation scheme discussed above was used to calculate the O,+ and NO+ densities during the nighttime in the E-region. Since the masses of O,+ and NO+ were almost the same, the vertical drift of these two ions species produced by the neutral winds were indistinguishable. The equation for the common characteristic curve was dz $ = V*(z) = C,(z)U,(4 + C,(x)U,(x) where U,(z) and U,(x) were the height varying north-south and east-west horizontal winds, z is the height coordinate, and C,(z), G&z) were height varying parameters calculated from the neutral atmospheric density. The ion densities were calculated
Certain zyzteme of first order partial diikential equetions
Vefticd
dnft
Vebaty,
m/xc
Jim, Chamctrist
163
KC ic
curves
Fig. 1. Characteristiccurves for sinusoidalvertical velocity. 0300 ES1 Wallops lsbnd I50
~-
-
140I
/'
-: : I30 --
\
E 1
120-
p
110 ~- _10
m.
too ---
90
- -I
VII
m&c I
’ /
I
1(111/ IO'
Electron
densely.
IO4
cmm3
Fig. 2. Height variation of (a) initial ion density; (b) measuredelectrondmity; and (c) calculated ionization density for 03 00 EST wind measured at Wallops Island, Virginia, 22 February 1968.
from the respective coupled continuity eqrmtions (CHENand HARRIS, 1971).
y y
= QA4 - ~Aw2+1~, - 45mwh+1
- ; w,+lv,(z))
= &A4 - %w[No+l~, + ~rJmrO,+l - ~([NO+lV,(z))
(10) (11)
164
W. M. CEEN, R. D. HARBIS and 6. D. WATSON
where [O,+] and [NO+] stand for the height varying O,+ and NO+ densities, Qi(z) and a,(z) are the ion production functions md recombination coefficients, and N, is the height varying electron density. The modeling assumes that only the height variations are important. The altitudes at which the ion densities were desired were initially specified. From each altitude the characteristic (9) was integrated backward to t = 0 from the time at which the densities were to be calculated. At t = 0 the initial densities, which were known at all altitudes for all species, were substituted into the subsidiary equations (equation (6)) and were integrated forward along the characteristic curves to yield the desired NO+ and O,+ densities. Figure 2 illustrates a typical electron density profile (NO+ + O,+ = NJ that resulted from the calculations. While it is not our purpose to discuss the proflle, it can be seen that reasonable agreement w&a found between the theoretical model and the experimental profile. SUMMARY
A numerical method for the calculation of simultaneous first order hyperbolic partial differential equations has been presented. The method is a modification of Lagrange’s method of converting to a system of ordinary differential equations and is restricted to the case when the characteristic curves are essentially identical. For a fluid medium this means that the velocities of all the component fluids are nearly the same. Important features of the calculation technique were the methods of finding the solution at desired points in space and time and the treatment of boundaries when a finite region of space was considered. The method was applied to the modeling of a partially ionized gas mixture. The numerical computations showed that when applicable, this method of solving simultaneous partial differential equations is stable, economical and accurate. REFERENCES COURANT
R.and HILBERT
D.
CHEN W. M. and H.ARRISR. D. LISTER bf.
1962
1971 1960
Methode of Mathematical Physics, Vol. II. Partial Differential Equations. Inter.
scienceNew York. J. Atmph. Terr. Phyu. 88, 1193.
Mathematical Methods for LGpitcd Computem (Edited by A. RAIBTON and
H. S. WIFF). Wiley, New York.
ad Termmtrld Phy&x. 1072,Vol. 84,pp. 160-164.%~'tPmon Prey.
SHORT PAPER An extensionof Lafmmge’s method to oewtain SJrstemsof ii& order ptutid difkalu equatiolla W. M. &TEN* and R. D. HAIWS Department of Electrical EngineeringUtah State University, U.S.A. and J. D.
WATSON
M&hem&ice Department Utah State University, U.S.A. (Received3 Mcry 1971) Ah&a&--A modificationof the method of characterieticehas been developedto solve systems of first orderhyperbolicpartial differentialequationawhere the characteristiccurves are essentially identical. This method ie applicable to fluid mixturea whose component fluid velocities are the same. Techniqueswere developed to find the solutionsat desiredpoints in space and time and to handle the conditionsat the boundaries. Results of the applicationof this method to the solution of ion density distributions in the Earth’s ionosphereare illustrated. 1. INTBODUCTION
of the interaction processes which occur in a fluid mixture require knowledge of the density distribution of each component of the mixture. These individual densities are given by the simultaneous solution of the set of continuity equations, one for each fluid component, which describe the component fluids. In general this system of equations is difficult to solve, especially if one or more of the fluid densities are changing rapidly. The purpose of this paper is to describe a numerical technique which is quite simple to use and has applicability to situations where the individual fluid components have the same or nearly the same velocities, and where the interaction and transport processes make diffusion an unimportant factor. The example cited in this paper is the calculation of ion densities in the E-region of the ionosphere (90-140 km). The primary fluids are the ion species of molecular oxygen O,+ and nitric oxide NO +. Each of the ion species are produced and destroyed by collisions with other atmospheric gases as well as by ionization from external radiation. In the E-region, wind transport processes are also an important factor in the ion density distribution. In order to theoretically model this ionospheric plasma, it is necessary to calculate the density distribution of each of the ion species under the influence of the reaction and transport processes. Since O,+ and NO+ have approximately the same mass, a modification of the method of characteristics has been developed to solve for the ion densities. THE STUDY
* Present address: Department of Physics, University of Prince Edward Island, Charlottetown, Prince Edward Island, Canada. 11
169
W. M. C!EBN, R. D. HARBIS end J. D. WATSON
160
The mathematical theory is considered in Section 2 and the results of actual calculations of the ionospheric ion density distribution are given in Section 3. 2. THEORY The numerical solution of the single continuity equation
qJ + v . (V(x, f) iv@,t))
= P(x, t, N)
(1)
with initial condition N(x, to) = h(x)
(2)
can be solved by Lagrange’s method of converting to a system of ordinary differential equations (COURANT and HILBEBT,1962) for which standard numerical programs are available to do the integration. To solve equation (1) by Lagra@s method we must solve the corresponding subsidiary equations dx - = V(x, t) dt
dN = P(x, t, N) - N(x, t)(V dt
l
V(x, t)).
The initial condition iV(x, to) = h(x) must be known at every point over the region of interest. The solution of (3) yields the so-called charaoterikic curve, which can be seen to be independent of N(x, t). Once the characteristic curve is known, then N(x, t) can be determined along it by solving (4). For a fluid mixture composed of ti components, a complete solution requires n continuity equations of the form
awx,t) + v at
l
(V,(x, wtx, f)) = -zp,(x, t, N,, NW- * * 9 N,)
i=l,2,...,n
and initial conditions N,(x, to) = M)
i = 1, 2, . . . , 12.
Extending Lagrange’s method to the system of continuity equations requires the solution of the corresponding subsidiary equations,
dx, - = dt
~*(x, t) dt
Vi(X, t)
= FAX,4 N,, iv*) - Nl(X, w
l
VAX,0)
i = 1,2,. . . ,n
(6)
where N,(x, to), i = 1, 2, . . . , n are initial conditions. In many applications the velocities of the fluid particles can be found in terms of the external forces, and are independent of the velocities and densities of the other species. When this is so the characteristic curves are independent. This means that throngh any point (z, t) it is neoessary to integrate (6) for each of the n speoies along its partienk ~&a&&tic
Certain eyetemf3of fht
order partial differentialequations
101
curve from 1, to time t = 7. Each of the component densities is evaluated on its own characteristic curve. Since the characteristic curves are different, the densities are evaluated at different spatial positions. To be correct, however, densities of interacting species must all be evaluated at the same position. This limitation of the method of characteristics can be eased somewhat by an iteration technique, but because of its difhculty it will have to be the subject of another paper. A numerical scheme to solve two simultaneous hyperbolic differential equations by Lagrange’s method has been reported by LISTER (1965). If the velocities of the individual fluid components of the mixture are nearly the same, then the chatacteristic curves are almost identical. For this condition one can then utilize a single characteristic curve and all the densities are found along this curve. The set of equations (6) reduces to the single equation dx - = V(x, t). dt
(7)
The shapes of the characteristic curves obviously depend upon the fluid velocity profile. For convergent or divergent flows the characteristic curves in general are not oonveniently oriented to allow density calculations over the complete region of interest. A technique to eliminate this problem allows the specification of points at which the final density values are desired. The solution technique then requires that the characteristic curves through each point of interest be integrated backwards in time to the initial starting time t,. But at t,, the densities iV(x, to) me known everywhere, so that a value of density is available to start the forward integration of (6). The densities are then found along the same characteristic curves, in particular at the positions that were initially specified. This method provides a check on the accuracy of the calculations because the backward integration started at the point (x, t) and the forward integration should return to the point (x, t). Very often it is impractical to model the entire fluid body because of limitations in size or complexity of the region. In these situations it is advantageous to choose a restricted region where the important reactions are taking place. But before numerical computations can be made, the conditions at the boundaries of this restricted region must be specified. If a characteristic curve either enters or leaves the region (i.e. passes through the boundary) at any time during the calculations this means physically that fluid is passing across the boundary, and in order to perform computations inside the region, the fluid flux must be known at the boundaries for all time. But this sort of information is not usually known. Thus it may be advantageous to choose or assume a closed system with no flux at the boundary. For certain cases this may also be a valid assumption if the flux is small or if the principal interactions of the fluid constituents take place inside the sub-region. The assumption that no fluid passes through the boundary implies that the characteristic curves do not pass through the boundaries. This can be implimented by allowing or forcing particle velocities to go to zero at the boundary. Zero velocities at the boundary, however, imply singular boundary conditions and instabilities in numerical calculations near the boundary. This can be seen by the following illustration in one space dimension: Assume z1 is a boundary point and that the velocity u(z) vanishes at z1 but that
W. M. CuEN,R. D. HAERISand J. D. WATSON
162 v(z)
2
0 and dv/dx > 0 near xi, but inside the region. The characteristic equation dx z = v(x)
which passes through point xg near xi at time t, can be approximated by the first two terms of a Taylor’s series dx dv z = v(x1) + z ;t_ 21(x - x1).
since2)(x1) =
0, the approximate velocity is given by the linear, constant-ooeflicient
equation
dx --dt
dv dx =rZ1
(x - Xi) = 0
whioh is valid for x near xi. This has solution x-x2,=Cexp
(t - t1) z = .%I
(8)
where C is a suitable constant. If we go backward in time from t,, x leaves x2 and approaohes x1 asymptotically but never reaches it. Thus the characteristic curves tend to bunch near the boundaries. Forward integration is then unstable in the sense that a small error in the position of the characteristic or the calculated velocity can cause large errors in the end positions of the charaoteristic curves. This is shown schematically in Fig. 1 where we have shown how the oharacteristics can bunch at a boundary where V(xi, t) = 0 but (dV/dx)I, _ 21# 0 at the boundary. One way to increase the stability of calculations near the boundaries is to increase the size of the region of interest so that the fluid velocity can be reduced gradually (dV/dx small) to zero at the boundaries. 3. APPLICATION The numerical calculation scheme discussed above was used to calculate the O,+ and NO+ densities during the nighttime in the E-region. Since the masses of O,+ and NO+ were almost the same, the vertical drift of these two ions species produced by the neutral winds were indistinguishable. The equation for the common characteristic curve was dz $ = V*(z) = C,(z)U,(4 + C,(x)U,(x) where U,(z) and U,(x) were the height varying north-south and east-west horizontal winds, z is the height coordinate, and C,(z), G&z) were height varying parameters calculated from the neutral atmospheric density. The ion densities were calculated
Certain zyzteme of first order partial diikential equetions
Vefticd
dnft
Vebaty,
m/xc
Jim, Chamctrist
163
KC ic
curves
Fig. 1. Characteristiccurves for sinusoidalvertical velocity. 0300 ES1 Wallops lsbnd I50
~-
-
140I
/'
-: : I30 --
\
E 1
120-
p
110 ~- _10
m.
too ---
90
- -I
VII
m&c I
’ /
I
1(111/ IO'
Electron
densely.
IO4
cmm3
Fig. 2. Height variation of (a) initial ion density; (b) measuredelectrondmity; and (c) calculated ionization density for 03 00 EST wind measured at Wallops Island, Virginia, 22 February 1968.
from the respective coupled continuity eqrmtions (CHENand HARRIS, 1971).
y y
= QA4 - ~Aw2+1~, - 45mwh+1
- ; w,+lv,(z))
= &A4 - %w[No+l~, + ~rJmrO,+l - ~([NO+lV,(z))
(10) (11)
164
W. M. CEEN, R. D. HARBIS and 6. D. WATSON
where [O,+] and [NO+] stand for the height varying O,+ and NO+ densities, Qi(z) and a,(z) are the ion production functions md recombination coefficients, and N, is the height varying electron density. The modeling assumes that only the height variations are important. The altitudes at which the ion densities were desired were initially specified. From each altitude the characteristic (9) was integrated backward to t = 0 from the time at which the densities were to be calculated. At t = 0 the initial densities, which were known at all altitudes for all species, were substituted into the subsidiary equations (equation (6)) and were integrated forward along the characteristic curves to yield the desired NO+ and O,+ densities. Figure 2 illustrates a typical electron density profile (NO+ + O,+ = NJ that resulted from the calculations. While it is not our purpose to discuss the proflle, it can be seen that reasonable agreement w&a found between the theoretical model and the experimental profile. SUMMARY
A numerical method for the calculation of simultaneous first order hyperbolic partial differential equations has been presented. The method is a modification of Lagrange’s method of converting to a system of ordinary differential equations and is restricted to the case when the characteristic curves are essentially identical. For a fluid medium this means that the velocities of all the component fluids are nearly the same. Important features of the calculation technique were the methods of finding the solution at desired points in space and time and the treatment of boundaries when a finite region of space was considered. The method was applied to the modeling of a partially ionized gas mixture. The numerical computations showed that when applicable, this method of solving simultaneous partial differential equations is stable, economical and accurate. REFERENCES COURANT
R.and HILBERT
D.
CHEN W. M. and H.ARRISR. D. LISTER bf.
1962
1971 1960
Methode of Mathematical Physics, Vol. II. Partial Differential Equations. Inter.
scienceNew York. J. Atmph. Terr. Phyu. 88, 1193.
Mathematical Methods for LGpitcd Computem (Edited by A. RAIBTON and
H. S. WIFF). Wiley, New York.