- Email: [email protected]
An improved biased difference scheme for the solution of hyperbolic equations by the method of lines
005%1354/82/020101-lOSO3.lWO @ 1982 Pergamon Press Ltd.
Compulerz & Chcmiiol Enginecting Vol. 6, No. 2, pp. 101-110. 1982 Printed in Great Britain.
AN IMPROVED BIASED DIFFERENCE SCHEME FOR THE SOLUTION OF HYPERBOLIC EQUATIONS BY THE METHOD OF LINES J. C. HEYDWEILLER* and H. S. PATEL Departmentof ChemicalEngineeringand MaterialsScience,SyracuseUniversity,Syracuse,NY 13210,U.S.A. (Received
30 Murch 1981;received for publicntion 20 July 1981)
Abstract-An upstream biased differencing scheme is presented for the solution of hyperbolic partial differential equations by the method of lines. The parameters in the difference formula are functions of the grid spacing and can give accuracy between first and second order. Testing has indicated that an order of 1.5 leads to a reasonable balance between stability and accuracy for typical problems. The mmethod is successfully demonstrated on hyperbolic equations which do not have shock discontinuities and is shown to be an effective means of solving coupled sets of hyperbolic and parabolic equations. Scope-While parabolic partial differential equations can be efficiently solved by the method of lines by using centered, second-order finite differences to replace the spatial derivatives, substitution of these differences into hyperbolic equations leads to an unstable solution procedure. The use of first-order, backward differences renders a stable procedure but their relative inaccuracy results in an inefficient solution. The biased difference formula presented in this paper is obtained by taking a weighted average of forward and backward differences with the weighting chosen to give a balance between accuracy and stability. In order to ensure stability, the weights for the three point difference formula are selected so that the coefficient for the point which is upstream relative to the movement is larger than the coefficient for the downstream point. The coefficients in the new formula are prescribed functions of the grid-spacing used for the spatial discretization. Once the grid size is selected to give a truncation error of the desired order of accuracy, the coefficients in the formula are fixed. Since no numerical experimentation is required to determine the best values of the parameters, the biased difference formula is convenient for the user. Conclusionsand SIgnhIcance-The biased difference scheme proposed in this paper overcomes the instability associated with second-order, centered differences while at the same time giving accuracy between first and second order. Numerical testing on three problems with smooth solutions has indicated that an order of 1.5 leads to a reasonable balance between stability and accuracy. Many engineering problems involve only the first partial derivative of the spatial variable and thus have only one spatial boundary condition. A three point difference formula, such as the proposed biased difference scheme, cannot simply be substituted into the differential equation at the point described by a null boundary condition. The dilliculty is overcome by integrating the differential equation over the appropriate interval adjacent to the boundary point to derive the correct difference equation. This approach preserves both the material balance and the order of accuracy, and also provides a convenient method for using a three-point difference formula with any type of boundary condition. The biased difference procedure is especially attractive for solving a coupled system of hyperbolic and parabolic equations, a common occurrence in the mathematical models of many chemical processes described by distributed parameter models. If such a coupled system were to be solved by standard finite techniques, considerable programming effort and computer time would be required to develop a stable procedure of comparable accuracy. INTRODUCTION
A convenient and powerful method of solving a system of nonlinear partial differential equations is to discretize all but one of the independent variables and then solve the resulting set of ordinary differential equations numerically. The advantage of this procedure, known as the method of lines[l,2], is that sophisticated computer programs based on higher-order time integration formulas can be employed. Gear-type integration programs [3,4], which can solve stiff as well as non-stiff *Authorto whom correspondence should he addressed.
systems of equation, feature automatic step-size adjustments and integration-order selection to maintain a user-specified error tolerance and to solve the equations with near optimal efficiency. Sincovec & Madsen[S] have developed and successfully tested a computer program for solving sets of parabolic partial differential equations on one spatial domain by the method of lines. While they were able to use centered, second-order finite differences to replace the spatial derivatives, substitution of these differences into hyperbolic equations would lead to an unstable solution procedure. The upstream-biased difference formula presented in this paper leads to a 101
J. C. HEYDWEILLER and
102
stable solution for hyperbolic equations which have smooth solutions thus permitting sets of coupled hyperbolic and parabolic equations to be solved by the method of lines. Such coupling is common in the mathematical models of many chemical processes, especially when two phases are present. While chemical engineers rarely deal with equations of the form of the wave equation (utl = uXX),equations of the form of the advection equation (Us= - u,) are often encountered. One way to solve the advection equation by the method of lines would be to approximate the spatial derivative by the following backward difference:
(1) where the subscript i denotes the spatial grid point and the last term represents the truncation error. The resulting solution would be stable but relatively inefficient because a small value of Ax would be required for accuracy due to the first-order error term. If centered differences are used to reduce the truncation error, an artificial dissipation term must be added to provide stability: Ui+;i&l
+ DC ui+l
-2Ui +
l&l (2)
(uxh'
(Ax)’
.
A large dissipation coefficient will destroy the secondorder accuracy while too small a value of DC will lead to instability; the best value must be determined by numerical experimentation. The upstream-biased difference formula proposed by Heydweiller & Sincovec [6], AUi+l t (Mr)i'
H. S. PATEL
=&[(l
tZEUi_(1 t
-C)Ui+l
(4)
E)Ui-11.
While the solution of the advection equation with the parameter l set equal to 0.3 was accurate and stable, the computation time required was 50% greater than the time required using either backward or upstream biased differences. This increase was due to the fact that using Eq. (4) necessitates the solving of a matrix problem at each time step even though a non-stiff (explicit) integrator was employed. If Eq. (4) were used to solve a system of coupled hyperbolic and parabolic equations, which would require the use of a stiff (implicit) integrator, then the computational effort expended would be greater than if Eq. (1) or (3) were used since the bandwidth of the Jacobian matrix would be larger. It should also be noted that the error criterion, given by
which was employed by Carver and Hinds to evaluate the relative accuracies of various difference formulas, will be large if the calculated solution lags the exact solution. This error criterion is relatively insensitive to variation in p in the upstream-biased difference formula since a change in this parameter which improves the peak height also increases the time lag. New Formula The improved upstream-biased difference scheme to be evaluated in this paper was originally derived by combining an artificial dissipation term with the biased difference formula [6] given by Eq. (3). For the equation U, = cou, ,
CUi- BUi-1
(At B)Ax
(5)
the resulting approximation for the first derivative is:
-@$(U,)i;
C’E B-A=+
1,
(3) '
(&Ii
reduces the first-order error term by a factor of At B. When At B = 1, the formula reduces to a backward difference while it approaches a centered difference as At B becomes large. While increasing the parameter
(uX)ir[AA~t(At
’
CD
B
Ui+t-2s
t Ui-1
(Ax)’
’
(6) where C = (- co/lcoJ).This equation may be rewritten with the error terms in the form:
B)CDB]Ui+rl[CAX-~(A~B)CD~]U~-[BAX-(A+B)CD~]U~-~ (A + B)(Ax)*
Ax DB-2(AtB)
C
AUi+l t CUi- BUi-I+ (At B)Ax
-@$ [ 1(Ux.x)i
At B improves the accuracy, too large a value will lead to the instability associated with the centered difference. The best value of At B depends on the size of Ax and must be found by experimentation. The Lax-Wendroff [7] approach is often used to stabilize solutions by explicit finite difference techniques. This approach is not easily adapted to the method of lines because it involves manipulation of the time derivatives. Carver and Hinds[l] have proposed the following formula for use with the method of lines:
(U*x,)i
t F
(7)
(U,)i].
From the above form, it is obvious that the coefficient for the first error term (u,) can easily be controlled by the arbitrary choices of DB and At B. For DB = 0, the error is simply that of the original upstream-biased difference [6] and thus can be. controlled by the choice of the parameter A t B. When DB = Ax/2(A t B), the firstorder error term vanishes but Eq. (7) then reduces to a centered difference with its inherent instability. The accuracy can be set to an intermediate value by choosing Ax DI~-~(~+B)=-~;
(Ax)”
l
(8)
103
An improved biased difference scheme
The finite difference procedure for solving Eq. (5) is given by:
With this choice, Eq. (7) becomes:
au;,t CUi”‘+ hurt
uin+’ - un
ht=Co
(9)
a = 1 _ WAX)‘-
b
=
l
+
’
ECU-’
,,; = ,,n ei@x,
P
’
(12)
where j = d( - 1) and y is a complex constant. The von Neumann criterion for unconditional stability is that the magnitude of y must be less than one. For the implicit procedure, substitution of Eq. (12) into Eq. (11) with m = n + 1, leads to the requirement that
and c = 4C(Ax)‘-’ P
19(11)
where m = n for the explicit case and m = n t 1 for the implicit case. To use the Fourier method of analysis[ll], one assumes that the solution of the difference equation at the point Xi and time t. is given by:
where
P
2Ax
*
This equation has the same form as Eq. (3) since c = b - a and a t b = 2, and thus the stability analysis will be essentially the same. Before presenting an abbreviated form of that analysis, one should observe that the advantage of this new form is that the arbitrary parameter At B which is dependent on Ax has been replaced by the parameters p and P which are independent of Ax. Further, it will later be demonstrated that single values of these parameters (1.5 and 1.0, respectively) can be used for most problems on a normalized spatial domain. It is also interesting to note that Eq. (9) can be expressed in terms of a centered difference with a dissipation term as follows:
l-cr(l-cosqAx)>l,
(13)
where r = coAdAx. This condition will be satisfied for all non-zero values of q if c and co are opposite in sign. This implies that b > a for co < 0 and b < a for co > 0 and thus the upstream point is weighted more heavily than the downstream point in the difference formula. For the explicit procedure when m = n, the condition of marginal stability (Irl= 1) occurs when -C
r=
a2t b* t 2ab cos qhx *
(14)
Using the definitions of r and c, this may be rewritten as: (Ux)i
’
Ui+;hh-I
+
0,
Ui+j-2Ui
t
Ui-1
(Ax)*
’
(10) (15)
where 0, = - C(Ax)“/P and thus is not equal to DB. As it is difficult to analyze the stability of the method of lines procedure which uses a variable-order time integration, the analysis will be limited to the simplest explicit and implicit (in time) finite difference procedures. These procedures correspond to the explicit and implicit Euler’s methods which are the lowest order formulas in differential integrator the ordinary equation packages [3,4]. Stability for these procedures implies stability for the higher-order formulas because all of the formulas are stable for the same class of equations [9].
While this restriction on the size of the time step is more severe than for a simple backward difference (At5 Ax/Icol), the improved accuracy of the biased difference formula will lead to a more efficient solution procedure. This fact is demonstrated in Table I for p = 1.5 and P = 1.0, values which will later be shown to be reasonable. In the table, it is assumed that the amount of work for an explicit procedure is proportional to the number of equations multiplied by the number of time steps required, based on the maximum At allowed for stability.
Table 1. Comparison of backward difference to upstream biased difference (p = 1.5 and P = 1.0) for an explicit Upstream biased difference, eq. (9) procedure. Ratios Backward difference, Eq. (1)
1
Relative
Error Tolerance
/
Spatial
Grid,Ax
Maximum Time Step, At
Number of Equations
Work Rewire
10-l
1.08
l/2
0.928
1.856
1o-2
2.32
l/2
9.431
0.862
1o-3
5.00
l/2
0.200
0.400
1o-4
10.8
l/2
0.093
0.186
1o-5
23.2
l/2
0.043
0.086
lo-6
50.0
l/2
0.020
0.040
*Assumed to be number of equations multiplied by minimum number of time steps.
J. C. HEYDWEILLER and
104
When more than two digits of accuracy are required, the use of the upstream-biased difference formula will result in a considerable reduction in the computational effort. Results and Discussion The first problem on which the new formula was tested was the advection equation,
To compare the accuracy for various values of p and P, this equation was first solved for the movement of a triangular wave. The boundary condition used was u(O, t) = 0; t z 0,
(Mb)
and the initial condition was
WC)
u(x, 0) = uo(x).
The exact solution, u(x, t) = Uo(x-t), is shown as the triangular wave at time 0.4 in Fig. 1. As indicated in this figure, increasing p from 1.4 to 1.6 improves the peak height at the expense of increasing the time lag and having more and larger oscillations trailing the wave. This behavior is similar to the results presented by Heydweiller & Sincovec[6] showing the effect of increasing the parameter A + B in Eq. (3). The solutions using first-order backward differences and second-order centered differences are presented for comparison with the solutions using the upstream-biased differences which have accuracy O[(AX)~].Solutions for values of p less than 1.4 approach the backward difference solution while the solutions for values of p greater than 1.6 approach the centered difference solution. Although the value of 0.01 used for Ax in the solutions in Fig. 1 is not small enough to give much accuracy, it is sufficient for a visual comparison of the solutions. These results suggest that using p = 1.5 gives a reasonable balance between accuracy and stability. To show that this value of p gives reasonable solutions for any value of Ax, Eq. (16a) was also solved for a step input, i.e. u(O,t)=l;
t>o,
(16d)
u(x, 0) = 0;
t = 0.
We)
and
The results shown in Fig. 2 are for p = 1.5, P = 1 and Ax = 0.04, 0.02, 0.01 and 0.005. While the solutions for the larger values of Ax are not very accurate, they do remain stable as indicated by the amplitude of the oscillations decreasing with time. Solutions for larger values of P are not shown but had larger oscillations than the solutions for P equal to one. The solutions for the two cases of the advection equation considered show that using the upstream-biased difference given by Eq. (9) with p = 1.5 and P = 1 produces stable solutions with accuracy between first and second-order. This combination of values leads to the following simple formula:
tuxJi = (1 -~C~~(AX))S+I
+JCd/(AX)Ui
2Ax
- (I+
H. S. PATEL
where the sign of C is chosen to weight the upstream point more heavily than the downstream point. The advantage of fixing the values of p and P is that the choice of the size of Ax can be made solely on the basis of the desired accuracy. Once the size of Ax is determined from the truncation error estimate, the coefficients in Eq. (17) will provide a stable solution for problems with smooth solutions (this restriction will be discussed later). Before proceeding to the next test problem, it should be noted that while the results for the case with the triangular wave (Eqs. 16a, b and c) were not dependent on the treatment of the right boundary for the time period considered, the results for the case with the step input (Eqs. 16a, d and e) were highly dependent on how the null boundary condition was incorporated into the solution procedure. For example, if one incorrectly assumes that the first derivative is zero at the right boundary, then the wave would be reflected. Since the reflected wave would be traveling in the direction opposite to the one for which the upstream-biased difference procedure is stable, the solution would then become unstable. The correct treatment can be demonstrated by first observing that the upstream-biased differences for the interior points can be derived by integrating the partial differential equation as follows:
I
xi+(‘AJ2)
e
xi-(bAxl2)
dx = _
xi+(‘AdJz) e
I
Xi-(bAJ2)
a’
dx
(18)
ax
Using the trapezoidal rule for the integral on the 1.h.s. of the equation and making suitable numerical approximations, one obtains: a&+,
dui -J
+ (b - U)Ui
- b&-l
(a+b)Ax
IT--
1*
(19)
The r.h.s. is simply Eq. (9) since b - a = c = + 1 for this case and a t b = 2. At the right boundary, one should not integrate over a whole interval but rather should use the limits of integration indicated in Fig. 3. The integration given by: (20) produces the following difference formula to be used for the null boundary condition:
dur dr+-
(2 - ah - bu,_, bAx
1’
(21)
This reduces to a simple backward difference since 2a = b but it does preserve the O[(Ax)“]accuracy because the approximation is not used over an entire interval. In addition, the matrix structure using Eq. (21) preserves the so-called material balance. The integration approach to incorporating the boundary conditions makes the three-point, upstream-biased difference scheme as easy to use as the less accurate, two-point backward difference formula. The second example, due to Richtmyer[lO, 111, is
~CV\/(AX))U~-~ _ C(Axjl.5(u
)..
XXI.
c = r 1
(17)
105
An improved biased difference scheme
1 backward
difference) II
IO
t =7
t = .4
I II
II II 0
0
05
-Ai
L
0
00 L
0.2
(centered
06
1.0
II
0.2
06
difference)
Fig. I. Effect of p on solution of Eq. (16a-c)
1.0
+--k-+0
Ax
~0.02
Ax
q
0.005
04
O0
I##IIJ
02
0.6
IO
Fig. 2. Stable solutionsfor Eqs. 16(a,d, e) for various values of Ax.
given by:
au ai -u* -=--
( )
ax2
at
’
Wa)
with u(0, t) = 0
Wb)
and u(x,O)=x;
05x51.
(224
This nonlinear equation is very useful for measuring the accuracy of a numerical solution because it has the following analytical solution: U(X,t) = &.
(23)
The results, presented in Table 2 in terms of the maximum error at various solution times, are similar to those presented by Heydweiller & Sincovec[6] for Eq. (3) except that the treatment of the right boundary was modified. Rather than using a one-sided, second-order difference as was done previously, the null boundary condition was incorporated by integrating Eq. (22a) over the appropriate fraction of the last interval. This approach was also followed for the backward and centered difference formulas and thus the present results are somewhat different from those previously reported. For Ax = 0.01, the errors for the backward, biased and centered differences compare reasonably well with the respective orders of 1, 1.5 and 2. To test this method on solutions which develop shock discontinuities, the solution to the equations studied by Houghton & Kasahara[lZ, 131for the flow of a fluid over a barrier was attempted. The numerical solution using Eq. (17) was not successful, probably because the solu-
At improved biased
differencescheme
xI-2
xI-3
107
xI-I
Fig. 3. Interval of integrationfor null boundary condition.
Table 2. Comparisonof errors in solution for Eq. (22)
1 A; =' ;p;
1 .000286
.000359
.000302
.0002141
IA; 1 ;":
] .000828
.001031
.000860
.000626
.000042
.000055
.000139
.000221
AX =
.Ol
Centered difference
tion contains fronts moving in both the positive and negative x directions. This failure was similar to the unsuccessful attempt presented by Heydweiller & Sincovec[6] when they used Eq. (3) with various values of the parameter A + B. The reader is referred to their paper for a more complete description of the equations and the difficulties encountered in the numerical solution. The inability to solve this problem prompted the previous comment on limiting the application of this approach to equations which have smooth solutions, i.e. solutions without shock discontinuities. One of the main incentives for developing a scheme by which hyperbolic equations can be solved by the method of lines is that sets of coupled hyperbolic and parabolic equations can then be solved by this method. Such systems are frequently encountered in the dynamic simulation of chemical processes described by distributed parameter models. To demonstrate how easily coupled sets can be solved, the start-up of the chemical plant shown in Fig. 4 will be simulated. In this example, taken from Heydweiller ef al.[14], simple mathematical models are used for the processing units so that model development does not overshadow the presentation of the computational procedure. The equations and boundary conditions are given below in dimensionless form:
I I
Mixed Tank: -d+Am _dr
CSb$A’ + c9!he - (1-k G)b”l
(24) (25)
-d&z” _- C*I@ + c&z= - (l+ C9)9crnl. dr
(26)
Reactor: +&=c a7
c,$=h-4~“;
a”#‘A I-$-
$$-
h=O
C24Ab
(27a) (2%) WC) (28a) (28b)
108
J. C. HEYDWEILLER and H. S. PATEL
Fig. 4. Process flowsheet.
cretization of the absorber equations leads to a second
WC) set of ordinary differential equations which are also only (29a) (2%) WC) Absorber:
(W 0%) (314 (31b) (32a) OW (334 (W The values of the constants C, - C9 are 0.1,3,0.1,20, 10, 0.003, 75, 30.7 and 0.5, respectively. The parabolic partial differential equations describing the reactor were discretized using centered, second-order differences. As shown by Madsen and Sincovec[2], this discretization will result in a set of ordinary differential equations which are stiff. To avoid the instability associated with the use of centered differences of hyperbolic partial differential equations, the discretization of the absorber equations was effected with the improved formula for upstream-biased differencing. Since the flow in the absorber is countercurrent, the negative 5 direction is upstream relative to the gas phase but is downstream relative to the liquid phase. Equation (17) was thus used with C= t 1 for Eqs. (3Oa), (31a) and (32a) and with C = - 1 for Eq. (33a). The discretization of the reactor equations leads to a set of (dimensionless) time-dependent ordinary differential equations while the dis-
dependent on dimensionless time. These two sets which are coupled to each other through the boundary conditions at the entrance to the absorber are then easily linked to the ordinary differential equations describing the tank in which complete mixing of the feed and recycle streams occurs. The final step of the procedure is then to numerically integrate the three sets of ordinary differential equations simultaneously. For the simulation of a start-up, all the concentrations were initially set equal to zero. At r = 0, 4,.,’ and 4~’ were increased to 1.0 and 0.5, respectively. The results are presented as reactor and absorber concentration profiles in Fig. 5 and as time histories of the entrance and exit concentrations in Fig. 6. Component B, the limiting reactant, was almost completely consumed in the reactor and thus plots for B are not presented since they are relatively uninteresting. The solutions were obtained for Ah = 0.05 and A<=O.OZ; these values produce spatial truncation errors of nearly equal size. A slightly smaller value of 0.001 was specified for the relative time integration error in the Gear-type integration package written by Hindmarsh[lS] for systems of equations having a banded Jacobian. Because the stiff option employs an implicit procedure, the Jacobian matrix must be calculated. Rather than writing a lengthy subroutine to enter the analytic expressions for the partial derivatives, it was simpler to let the program generate numerical approximations internally. The entries outside the main band, which extended four entries above and below the main diagonal, caused by the recycle stream were smaller than the entries within the main band. Thus omitting the off-band terms in the Jacobian did not adversely affect the rate of convergence of the iterative procedure used to resolve the nonlinearity. The concentration profiles and histories presented in Figs. 5 and 6, respectively, are essentially the same as those previously presented [14]. One significant difference in the computation was that the value of A[ used in Eq. (17) to discretize the absorber equations was twice as large (0.02~s 0.01) as the value used in the previous study in which Eq. (3) was employed with A+B=3. The reason that Al could be increased without introducing oscillations into the solution was the correct incorporation of the null boundary conditions at the absorber exits into the difference approximations. Even though the increase cannot be attributed to the use
An improved biased difference
109
scheme
Rromaterf
Poromcter
r
04
0.6
06-
0
0.2
04
Parameter
*A
0.6
0.6
I.C
0
0.2
Parametw
T
0.6
T
0.0 t 0.61
04
0.2
rr_ 40
06-
17
6
0
0
0.2
04
0.6
0.6
E
Fig. 5. Dimensionless reactor/absorber
-
+:
06-
concentration
E profiles.
sc’
- 000 0 06
r
0;
9: -
0 04
Fig. 6. Dimensionless entrance/exit concentration histories.
of the new difference formula, the larger Al did result in a reduction of the number of discretized absorber equations from 400 to 200 and a subsequent reduction in the computational time. The ordinary differential equations were integrated to a dimensionless time of 40 with no
difficulty. The integrator took 210 time steps and evaluated the Jacobian matrix 24 times; these values compare to 450 and 31 in the previous study. The refinement of the upstream-biased difference formula has resulted in a simple method of producing
110
J. C. HEYDWEILLER and H. S. PATEL
stable solutions to hyperbolic partial differential equations by the method of lines. As demonstrated in the last
example, this method is especially attractive when the hyperbolic equations are coupled to parabolic equations. The coupling would make it difficult to obtain efficient solutions by standard finite difference techniques.
Acknowledgement-Acknowledgment is made to the Donors of the Petroleum Research Fund, administered by the American Chemical Society, for the partial support of this research. RJZFEReNCES 1. 0. A. Liskovets, The method of lines (Review). English translation in Difference Euuotions 1. 1308f1%5). 2. N. K. Madsen i R. F. Sincovec, The numerical method of lines for the solution of nonlinear partial differential equations. Computational Methods in Nonlinear Mechanics (Eds. J. T. Oden et al.), Texas Institute of Computational Mechanics, Austin, TX (1974). 3. C. W. Gear, The automatic integration of ordinary differential eauations. Commun. ACM 3. 176(1971). 4. A. C. *Hindmarsh,GEAR: Ordinary differential equation system solver. Rep. UCID30001, Rev. 3, Lawrence Livermore Laboratory (1974). 5. R. F. Sincovec & N. K. Madsen, Software for nonlinear partial differential equations. ACM Trans. Math. Software 1, 232 (1975).
6. J. C. Heydweiller & R. F. Sincovec, A stable difference scheme for the solution of hyperbolic equations using the method of lines. L Comput. Phys. 3,377 (1976). 7. P. D. Lax & B. Wendroff, Systems of conservation laws. Commun. Pure Appl. Math W, 217 (1960). 8. M. B. Carver & H. W. Hinds, The method of lines and the advection equation. Simulation 31,59 (1978). 9. C. W. Gear, Numerical Initial Value Problems in Ordinary Diferential Equations. Prentice-Hall, Englewood Cliffs, New Jersey (1971). 10. R. D. Richtmyer, The stability criterion of Godunov and Ryabenkii for difference schemes. Rep. NYO-1480-4, Courant Inst. of Math. Sciences, New York University, NY (1%4). II. R. D. Richtmyer & K. W. Morton, Diflerence methods for Initial Value Problems. 2nd Edn. Wiley, New York (1%7). 12. D. D. Houghton & A. Kasahara, Nonlinear shallow fluid flow over an isolated ridge. Commun. Pure Appi. Math. 21, 603 (1972). 13. A. Kasahara & D. D. Houghton, An example of nonunique, discontinuous solutions in fluid dynamics. J. Comput. Phys. 4, 377 (1%9). 14. J. C. Heydweiller, R. F. Sincovec & L. T. Fan, Dynamic simulation of chemical processes described by distributed and lumped parameter models. Comput. and Chem. Engng 1, 125(1977). 15. A. C. Hindmarsh, GEARB: Solution of ordinary differential equations having banded Jacobian. Rep. UCID-30059,Rev. 1, Lawrence Livermore Laboratory (1975).
Compulerz & Chcmiiol Enginecting Vol. 6, No. 2, pp. 101-110. 1982 Printed in Great Britain.
AN IMPROVED BIASED DIFFERENCE SCHEME FOR THE SOLUTION OF HYPERBOLIC EQUATIONS BY THE METHOD OF LINES J. C. HEYDWEILLER* and H. S. PATEL Departmentof ChemicalEngineeringand MaterialsScience,SyracuseUniversity,Syracuse,NY 13210,U.S.A. (Received
30 Murch 1981;received for publicntion 20 July 1981)
Abstract-An upstream biased differencing scheme is presented for the solution of hyperbolic partial differential equations by the method of lines. The parameters in the difference formula are functions of the grid spacing and can give accuracy between first and second order. Testing has indicated that an order of 1.5 leads to a reasonable balance between stability and accuracy for typical problems. The mmethod is successfully demonstrated on hyperbolic equations which do not have shock discontinuities and is shown to be an effective means of solving coupled sets of hyperbolic and parabolic equations. Scope-While parabolic partial differential equations can be efficiently solved by the method of lines by using centered, second-order finite differences to replace the spatial derivatives, substitution of these differences into hyperbolic equations leads to an unstable solution procedure. The use of first-order, backward differences renders a stable procedure but their relative inaccuracy results in an inefficient solution. The biased difference formula presented in this paper is obtained by taking a weighted average of forward and backward differences with the weighting chosen to give a balance between accuracy and stability. In order to ensure stability, the weights for the three point difference formula are selected so that the coefficient for the point which is upstream relative to the movement is larger than the coefficient for the downstream point. The coefficients in the new formula are prescribed functions of the grid-spacing used for the spatial discretization. Once the grid size is selected to give a truncation error of the desired order of accuracy, the coefficients in the formula are fixed. Since no numerical experimentation is required to determine the best values of the parameters, the biased difference formula is convenient for the user. Conclusionsand SIgnhIcance-The biased difference scheme proposed in this paper overcomes the instability associated with second-order, centered differences while at the same time giving accuracy between first and second order. Numerical testing on three problems with smooth solutions has indicated that an order of 1.5 leads to a reasonable balance between stability and accuracy. Many engineering problems involve only the first partial derivative of the spatial variable and thus have only one spatial boundary condition. A three point difference formula, such as the proposed biased difference scheme, cannot simply be substituted into the differential equation at the point described by a null boundary condition. The dilliculty is overcome by integrating the differential equation over the appropriate interval adjacent to the boundary point to derive the correct difference equation. This approach preserves both the material balance and the order of accuracy, and also provides a convenient method for using a three-point difference formula with any type of boundary condition. The biased difference procedure is especially attractive for solving a coupled system of hyperbolic and parabolic equations, a common occurrence in the mathematical models of many chemical processes described by distributed parameter models. If such a coupled system were to be solved by standard finite techniques, considerable programming effort and computer time would be required to develop a stable procedure of comparable accuracy. INTRODUCTION
A convenient and powerful method of solving a system of nonlinear partial differential equations is to discretize all but one of the independent variables and then solve the resulting set of ordinary differential equations numerically. The advantage of this procedure, known as the method of lines[l,2], is that sophisticated computer programs based on higher-order time integration formulas can be employed. Gear-type integration programs [3,4], which can solve stiff as well as non-stiff *Authorto whom correspondence should he addressed.
systems of equation, feature automatic step-size adjustments and integration-order selection to maintain a user-specified error tolerance and to solve the equations with near optimal efficiency. Sincovec & Madsen[S] have developed and successfully tested a computer program for solving sets of parabolic partial differential equations on one spatial domain by the method of lines. While they were able to use centered, second-order finite differences to replace the spatial derivatives, substitution of these differences into hyperbolic equations would lead to an unstable solution procedure. The upstream-biased difference formula presented in this paper leads to a 101
J. C. HEYDWEILLER and
102
stable solution for hyperbolic equations which have smooth solutions thus permitting sets of coupled hyperbolic and parabolic equations to be solved by the method of lines. Such coupling is common in the mathematical models of many chemical processes, especially when two phases are present. While chemical engineers rarely deal with equations of the form of the wave equation (utl = uXX),equations of the form of the advection equation (Us= - u,) are often encountered. One way to solve the advection equation by the method of lines would be to approximate the spatial derivative by the following backward difference:
(1) where the subscript i denotes the spatial grid point and the last term represents the truncation error. The resulting solution would be stable but relatively inefficient because a small value of Ax would be required for accuracy due to the first-order error term. If centered differences are used to reduce the truncation error, an artificial dissipation term must be added to provide stability: Ui+;i&l
+ DC ui+l
-2Ui +
l&l (2)
(uxh'
(Ax)’
.
A large dissipation coefficient will destroy the secondorder accuracy while too small a value of DC will lead to instability; the best value must be determined by numerical experimentation. The upstream-biased difference formula proposed by Heydweiller & Sincovec [6], AUi+l t (Mr)i'
H. S. PATEL
=&[(l
tZEUi_(1 t
-C)Ui+l
(4)
E)Ui-11.
While the solution of the advection equation with the parameter l set equal to 0.3 was accurate and stable, the computation time required was 50% greater than the time required using either backward or upstream biased differences. This increase was due to the fact that using Eq. (4) necessitates the solving of a matrix problem at each time step even though a non-stiff (explicit) integrator was employed. If Eq. (4) were used to solve a system of coupled hyperbolic and parabolic equations, which would require the use of a stiff (implicit) integrator, then the computational effort expended would be greater than if Eq. (1) or (3) were used since the bandwidth of the Jacobian matrix would be larger. It should also be noted that the error criterion, given by
which was employed by Carver and Hinds to evaluate the relative accuracies of various difference formulas, will be large if the calculated solution lags the exact solution. This error criterion is relatively insensitive to variation in p in the upstream-biased difference formula since a change in this parameter which improves the peak height also increases the time lag. New Formula The improved upstream-biased difference scheme to be evaluated in this paper was originally derived by combining an artificial dissipation term with the biased difference formula [6] given by Eq. (3). For the equation U, = cou, ,
CUi- BUi-1
(At B)Ax
(5)
the resulting approximation for the first derivative is:
-@$(U,)i;
C’E B-A=+
1,
(3) '
(&Ii
reduces the first-order error term by a factor of At B. When At B = 1, the formula reduces to a backward difference while it approaches a centered difference as At B becomes large. While increasing the parameter
(uX)ir[AA~t(At
’
CD
B
Ui+t-2s
t Ui-1
(Ax)’
’
(6) where C = (- co/lcoJ).This equation may be rewritten with the error terms in the form:
B)CDB]Ui+rl[CAX-~(A~B)CD~]U~-[BAX-(A+B)CD~]U~-~ (A + B)(Ax)*
Ax DB-2(AtB)
C
AUi+l t CUi- BUi-I+ (At B)Ax
-@$ [ 1(Ux.x)i
At B improves the accuracy, too large a value will lead to the instability associated with the centered difference. The best value of At B depends on the size of Ax and must be found by experimentation. The Lax-Wendroff [7] approach is often used to stabilize solutions by explicit finite difference techniques. This approach is not easily adapted to the method of lines because it involves manipulation of the time derivatives. Carver and Hinds[l] have proposed the following formula for use with the method of lines:
(U*x,)i
t F
(7)
(U,)i].
From the above form, it is obvious that the coefficient for the first error term (u,) can easily be controlled by the arbitrary choices of DB and At B. For DB = 0, the error is simply that of the original upstream-biased difference [6] and thus can be. controlled by the choice of the parameter A t B. When DB = Ax/2(A t B), the firstorder error term vanishes but Eq. (7) then reduces to a centered difference with its inherent instability. The accuracy can be set to an intermediate value by choosing Ax DI~-~(~+B)=-~;
(Ax)”
l
(8)
103
An improved biased difference scheme
The finite difference procedure for solving Eq. (5) is given by:
With this choice, Eq. (7) becomes:
au;,t CUi”‘+ hurt
uin+’ - un
ht=Co
(9)
a = 1 _ WAX)‘-
b
=
l
+
’
ECU-’
,,; = ,,n ei@x,
P
’
(12)
where j = d( - 1) and y is a complex constant. The von Neumann criterion for unconditional stability is that the magnitude of y must be less than one. For the implicit procedure, substitution of Eq. (12) into Eq. (11) with m = n + 1, leads to the requirement that
and c = 4C(Ax)‘-’ P
19(11)
where m = n for the explicit case and m = n t 1 for the implicit case. To use the Fourier method of analysis[ll], one assumes that the solution of the difference equation at the point Xi and time t. is given by:
where
P
2Ax
*
This equation has the same form as Eq. (3) since c = b - a and a t b = 2, and thus the stability analysis will be essentially the same. Before presenting an abbreviated form of that analysis, one should observe that the advantage of this new form is that the arbitrary parameter At B which is dependent on Ax has been replaced by the parameters p and P which are independent of Ax. Further, it will later be demonstrated that single values of these parameters (1.5 and 1.0, respectively) can be used for most problems on a normalized spatial domain. It is also interesting to note that Eq. (9) can be expressed in terms of a centered difference with a dissipation term as follows:
l-cr(l-cosqAx)>l,
(13)
where r = coAdAx. This condition will be satisfied for all non-zero values of q if c and co are opposite in sign. This implies that b > a for co < 0 and b < a for co > 0 and thus the upstream point is weighted more heavily than the downstream point in the difference formula. For the explicit procedure when m = n, the condition of marginal stability (Irl= 1) occurs when -C
r=
a2t b* t 2ab cos qhx *
(14)
Using the definitions of r and c, this may be rewritten as: (Ux)i
’
Ui+;hh-I
+
0,
Ui+j-2Ui
t
Ui-1
(Ax)*
’
(10) (15)
where 0, = - C(Ax)“/P and thus is not equal to DB. As it is difficult to analyze the stability of the method of lines procedure which uses a variable-order time integration, the analysis will be limited to the simplest explicit and implicit (in time) finite difference procedures. These procedures correspond to the explicit and implicit Euler’s methods which are the lowest order formulas in differential integrator the ordinary equation packages [3,4]. Stability for these procedures implies stability for the higher-order formulas because all of the formulas are stable for the same class of equations [9].
While this restriction on the size of the time step is more severe than for a simple backward difference (At5 Ax/Icol), the improved accuracy of the biased difference formula will lead to a more efficient solution procedure. This fact is demonstrated in Table I for p = 1.5 and P = 1.0, values which will later be shown to be reasonable. In the table, it is assumed that the amount of work for an explicit procedure is proportional to the number of equations multiplied by the number of time steps required, based on the maximum At allowed for stability.
Table 1. Comparison of backward difference to upstream biased difference (p = 1.5 and P = 1.0) for an explicit Upstream biased difference, eq. (9) procedure. Ratios Backward difference, Eq. (1)
1
Relative
Error Tolerance
/
Spatial
Grid,Ax
Maximum Time Step, At
Number of Equations
Work Rewire
10-l
1.08
l/2
0.928
1.856
1o-2
2.32
l/2
9.431
0.862
1o-3
5.00
l/2
0.200
0.400
1o-4
10.8
l/2
0.093
0.186
1o-5
23.2
l/2
0.043
0.086
lo-6
50.0
l/2
0.020
0.040
*Assumed to be number of equations multiplied by minimum number of time steps.
J. C. HEYDWEILLER and
104
When more than two digits of accuracy are required, the use of the upstream-biased difference formula will result in a considerable reduction in the computational effort. Results and Discussion The first problem on which the new formula was tested was the advection equation,
To compare the accuracy for various values of p and P, this equation was first solved for the movement of a triangular wave. The boundary condition used was u(O, t) = 0; t z 0,
(Mb)
and the initial condition was
WC)
u(x, 0) = uo(x).
The exact solution, u(x, t) = Uo(x-t), is shown as the triangular wave at time 0.4 in Fig. 1. As indicated in this figure, increasing p from 1.4 to 1.6 improves the peak height at the expense of increasing the time lag and having more and larger oscillations trailing the wave. This behavior is similar to the results presented by Heydweiller & Sincovec[6] showing the effect of increasing the parameter A + B in Eq. (3). The solutions using first-order backward differences and second-order centered differences are presented for comparison with the solutions using the upstream-biased differences which have accuracy O[(AX)~].Solutions for values of p less than 1.4 approach the backward difference solution while the solutions for values of p greater than 1.6 approach the centered difference solution. Although the value of 0.01 used for Ax in the solutions in Fig. 1 is not small enough to give much accuracy, it is sufficient for a visual comparison of the solutions. These results suggest that using p = 1.5 gives a reasonable balance between accuracy and stability. To show that this value of p gives reasonable solutions for any value of Ax, Eq. (16a) was also solved for a step input, i.e. u(O,t)=l;
t>o,
(16d)
u(x, 0) = 0;
t = 0.
We)
and
The results shown in Fig. 2 are for p = 1.5, P = 1 and Ax = 0.04, 0.02, 0.01 and 0.005. While the solutions for the larger values of Ax are not very accurate, they do remain stable as indicated by the amplitude of the oscillations decreasing with time. Solutions for larger values of P are not shown but had larger oscillations than the solutions for P equal to one. The solutions for the two cases of the advection equation considered show that using the upstream-biased difference given by Eq. (9) with p = 1.5 and P = 1 produces stable solutions with accuracy between first and second-order. This combination of values leads to the following simple formula:
tuxJi = (1 -~C~~(AX))S+I
+JCd/(AX)Ui
2Ax
- (I+
H. S. PATEL
where the sign of C is chosen to weight the upstream point more heavily than the downstream point. The advantage of fixing the values of p and P is that the choice of the size of Ax can be made solely on the basis of the desired accuracy. Once the size of Ax is determined from the truncation error estimate, the coefficients in Eq. (17) will provide a stable solution for problems with smooth solutions (this restriction will be discussed later). Before proceeding to the next test problem, it should be noted that while the results for the case with the triangular wave (Eqs. 16a, b and c) were not dependent on the treatment of the right boundary for the time period considered, the results for the case with the step input (Eqs. 16a, d and e) were highly dependent on how the null boundary condition was incorporated into the solution procedure. For example, if one incorrectly assumes that the first derivative is zero at the right boundary, then the wave would be reflected. Since the reflected wave would be traveling in the direction opposite to the one for which the upstream-biased difference procedure is stable, the solution would then become unstable. The correct treatment can be demonstrated by first observing that the upstream-biased differences for the interior points can be derived by integrating the partial differential equation as follows:
I
xi+(‘AJ2)
e
xi-(bAxl2)
dx = _
xi+(‘AdJz) e
I
Xi-(bAJ2)
a’
dx
(18)
ax
Using the trapezoidal rule for the integral on the 1.h.s. of the equation and making suitable numerical approximations, one obtains: a&+,
dui -J
+ (b - U)Ui
- b&-l
(a+b)Ax
IT--
1*
(19)
The r.h.s. is simply Eq. (9) since b - a = c = + 1 for this case and a t b = 2. At the right boundary, one should not integrate over a whole interval but rather should use the limits of integration indicated in Fig. 3. The integration given by: (20) produces the following difference formula to be used for the null boundary condition:
dur dr+-
(2 - ah - bu,_, bAx
1’
(21)
This reduces to a simple backward difference since 2a = b but it does preserve the O[(Ax)“]accuracy because the approximation is not used over an entire interval. In addition, the matrix structure using Eq. (21) preserves the so-called material balance. The integration approach to incorporating the boundary conditions makes the three-point, upstream-biased difference scheme as easy to use as the less accurate, two-point backward difference formula. The second example, due to Richtmyer[lO, 111, is
~CV\/(AX))U~-~ _ C(Axjl.5(u
)..
XXI.
c = r 1
(17)
105
An improved biased difference scheme
1 backward
difference) II
IO
t =7
t = .4
I II
II II 0
0
05
-Ai
L
0
00 L
0.2
(centered
06
1.0
II
0.2
06
difference)
Fig. I. Effect of p on solution of Eq. (16a-c)
1.0
+--k-+0
Ax
~0.02
Ax
q
0.005
04
O0
I##IIJ
02
0.6
IO
Fig. 2. Stable solutionsfor Eqs. 16(a,d, e) for various values of Ax.
given by:
au ai -u* -=--
( )
ax2
at
’
Wa)
with u(0, t) = 0
Wb)
and u(x,O)=x;
05x51.
(224
This nonlinear equation is very useful for measuring the accuracy of a numerical solution because it has the following analytical solution: U(X,t) = &.
(23)
The results, presented in Table 2 in terms of the maximum error at various solution times, are similar to those presented by Heydweiller & Sincovec[6] for Eq. (3) except that the treatment of the right boundary was modified. Rather than using a one-sided, second-order difference as was done previously, the null boundary condition was incorporated by integrating Eq. (22a) over the appropriate fraction of the last interval. This approach was also followed for the backward and centered difference formulas and thus the present results are somewhat different from those previously reported. For Ax = 0.01, the errors for the backward, biased and centered differences compare reasonably well with the respective orders of 1, 1.5 and 2. To test this method on solutions which develop shock discontinuities, the solution to the equations studied by Houghton & Kasahara[lZ, 131for the flow of a fluid over a barrier was attempted. The numerical solution using Eq. (17) was not successful, probably because the solu-
At improved biased
differencescheme
xI-2
xI-3
107
xI-I
Fig. 3. Interval of integrationfor null boundary condition.
Table 2. Comparisonof errors in solution for Eq. (22)
1 A; =' ;p;
1 .000286
.000359
.000302
.0002141
IA; 1 ;":
] .000828
.001031
.000860
.000626
.000042
.000055
.000139
.000221
AX =
.Ol
Centered difference
tion contains fronts moving in both the positive and negative x directions. This failure was similar to the unsuccessful attempt presented by Heydweiller & Sincovec[6] when they used Eq. (3) with various values of the parameter A + B. The reader is referred to their paper for a more complete description of the equations and the difficulties encountered in the numerical solution. The inability to solve this problem prompted the previous comment on limiting the application of this approach to equations which have smooth solutions, i.e. solutions without shock discontinuities. One of the main incentives for developing a scheme by which hyperbolic equations can be solved by the method of lines is that sets of coupled hyperbolic and parabolic equations can then be solved by this method. Such systems are frequently encountered in the dynamic simulation of chemical processes described by distributed parameter models. To demonstrate how easily coupled sets can be solved, the start-up of the chemical plant shown in Fig. 4 will be simulated. In this example, taken from Heydweiller ef al.[14], simple mathematical models are used for the processing units so that model development does not overshadow the presentation of the computational procedure. The equations and boundary conditions are given below in dimensionless form:
I I
Mixed Tank: -d+Am _dr
CSb$A’ + c9!he - (1-k G)b”l
(24) (25)
-d&z” _- C*I@ + c&z= - (l+ C9)9crnl. dr
(26)
Reactor: +&=c a7
c,$=h-4~“;
a”#‘A I-$-
$$-
h=O
C24Ab
(27a) (2%) WC) (28a) (28b)
108
J. C. HEYDWEILLER and H. S. PATEL
Fig. 4. Process flowsheet.
cretization of the absorber equations leads to a second
WC) set of ordinary differential equations which are also only (29a) (2%) WC) Absorber:
(W 0%) (314 (31b) (32a) OW (334 (W The values of the constants C, - C9 are 0.1,3,0.1,20, 10, 0.003, 75, 30.7 and 0.5, respectively. The parabolic partial differential equations describing the reactor were discretized using centered, second-order differences. As shown by Madsen and Sincovec[2], this discretization will result in a set of ordinary differential equations which are stiff. To avoid the instability associated with the use of centered differences of hyperbolic partial differential equations, the discretization of the absorber equations was effected with the improved formula for upstream-biased differencing. Since the flow in the absorber is countercurrent, the negative 5 direction is upstream relative to the gas phase but is downstream relative to the liquid phase. Equation (17) was thus used with C= t 1 for Eqs. (3Oa), (31a) and (32a) and with C = - 1 for Eq. (33a). The discretization of the reactor equations leads to a set of (dimensionless) time-dependent ordinary differential equations while the dis-
dependent on dimensionless time. These two sets which are coupled to each other through the boundary conditions at the entrance to the absorber are then easily linked to the ordinary differential equations describing the tank in which complete mixing of the feed and recycle streams occurs. The final step of the procedure is then to numerically integrate the three sets of ordinary differential equations simultaneously. For the simulation of a start-up, all the concentrations were initially set equal to zero. At r = 0, 4,.,’ and 4~’ were increased to 1.0 and 0.5, respectively. The results are presented as reactor and absorber concentration profiles in Fig. 5 and as time histories of the entrance and exit concentrations in Fig. 6. Component B, the limiting reactant, was almost completely consumed in the reactor and thus plots for B are not presented since they are relatively uninteresting. The solutions were obtained for Ah = 0.05 and A<=O.OZ; these values produce spatial truncation errors of nearly equal size. A slightly smaller value of 0.001 was specified for the relative time integration error in the Gear-type integration package written by Hindmarsh[lS] for systems of equations having a banded Jacobian. Because the stiff option employs an implicit procedure, the Jacobian matrix must be calculated. Rather than writing a lengthy subroutine to enter the analytic expressions for the partial derivatives, it was simpler to let the program generate numerical approximations internally. The entries outside the main band, which extended four entries above and below the main diagonal, caused by the recycle stream were smaller than the entries within the main band. Thus omitting the off-band terms in the Jacobian did not adversely affect the rate of convergence of the iterative procedure used to resolve the nonlinearity. The concentration profiles and histories presented in Figs. 5 and 6, respectively, are essentially the same as those previously presented [14]. One significant difference in the computation was that the value of A[ used in Eq. (17) to discretize the absorber equations was twice as large (0.02~s 0.01) as the value used in the previous study in which Eq. (3) was employed with A+B=3. The reason that Al could be increased without introducing oscillations into the solution was the correct incorporation of the null boundary conditions at the absorber exits into the difference approximations. Even though the increase cannot be attributed to the use
An improved biased difference
109
scheme
Rromaterf
Poromcter
r
04
0.6
06-
0
0.2
04
Parameter
*A
0.6
0.6
I.C
0
0.2
Parametw
T
0.6
T
0.0 t 0.61
04
0.2
rr_ 40
06-
17
6
0
0
0.2
04
0.6
0.6
E
Fig. 5. Dimensionless reactor/absorber
-
+:
06-
concentration
E profiles.
sc’
- 000 0 06
r
0;
9: -
0 04
Fig. 6. Dimensionless entrance/exit concentration histories.
of the new difference formula, the larger Al did result in a reduction of the number of discretized absorber equations from 400 to 200 and a subsequent reduction in the computational time. The ordinary differential equations were integrated to a dimensionless time of 40 with no
difficulty. The integrator took 210 time steps and evaluated the Jacobian matrix 24 times; these values compare to 450 and 31 in the previous study. The refinement of the upstream-biased difference formula has resulted in a simple method of producing
110
J. C. HEYDWEILLER and H. S. PATEL
stable solutions to hyperbolic partial differential equations by the method of lines. As demonstrated in the last
example, this method is especially attractive when the hyperbolic equations are coupled to parabolic equations. The coupling would make it difficult to obtain efficient solutions by standard finite difference techniques.
Acknowledgement-Acknowledgment is made to the Donors of the Petroleum Research Fund, administered by the American Chemical Society, for the partial support of this research. RJZFEReNCES 1. 0. A. Liskovets, The method of lines (Review). English translation in Difference Euuotions 1. 1308f1%5). 2. N. K. Madsen i R. F. Sincovec, The numerical method of lines for the solution of nonlinear partial differential equations. Computational Methods in Nonlinear Mechanics (Eds. J. T. Oden et al.), Texas Institute of Computational Mechanics, Austin, TX (1974). 3. C. W. Gear, The automatic integration of ordinary differential eauations. Commun. ACM 3. 176(1971). 4. A. C. *Hindmarsh,GEAR: Ordinary differential equation system solver. Rep. UCID30001, Rev. 3, Lawrence Livermore Laboratory (1974). 5. R. F. Sincovec & N. K. Madsen, Software for nonlinear partial differential equations. ACM Trans. Math. Software 1, 232 (1975).
6. J. C. Heydweiller & R. F. Sincovec, A stable difference scheme for the solution of hyperbolic equations using the method of lines. L Comput. Phys. 3,377 (1976). 7. P. D. Lax & B. Wendroff, Systems of conservation laws. Commun. Pure Appl. Math W, 217 (1960). 8. M. B. Carver & H. W. Hinds, The method of lines and the advection equation. Simulation 31,59 (1978). 9. C. W. Gear, Numerical Initial Value Problems in Ordinary Diferential Equations. Prentice-Hall, Englewood Cliffs, New Jersey (1971). 10. R. D. Richtmyer, The stability criterion of Godunov and Ryabenkii for difference schemes. Rep. NYO-1480-4, Courant Inst. of Math. Sciences, New York University, NY (1%4). II. R. D. Richtmyer & K. W. Morton, Diflerence methods for Initial Value Problems. 2nd Edn. Wiley, New York (1%7). 12. D. D. Houghton & A. Kasahara, Nonlinear shallow fluid flow over an isolated ridge. Commun. Pure Appi. Math. 21, 603 (1972). 13. A. Kasahara & D. D. Houghton, An example of nonunique, discontinuous solutions in fluid dynamics. J. Comput. Phys. 4, 377 (1%9). 14. J. C. Heydweiller, R. F. Sincovec & L. T. Fan, Dynamic simulation of chemical processes described by distributed and lumped parameter models. Comput. and Chem. Engng 1, 125(1977). 15. A. C. Hindmarsh, GEARB: Solution of ordinary differential equations having banded Jacobian. Rep. UCID-30059,Rev. 1, Lawrence Livermore Laboratory (1975).