- Email: [email protected]
Program for stability and accuracy analysis of finite difference methods
Computer Physics Communications 34 (1984) 175—186 North-Holland, Amsterdam
175
PROGRAM FOR STABILITY AND ACCURACY ANALYSIS OF FINITE DIFFERENCE METHODS R. LISKA Faculty of Nuclear Science and Physical Engineering~Technical University of Prague, Bi~ehovà7, 115 19 Praha 1, Czechoslovakia Received 8 February 1984; in revised form 4 July 1984
PROGRAM SUMMARY
equations, described in ref. [1], is applied. The program calculates a modified equation corresponding to a difference scheme
Title ofprogram: MODST
which approximates a partial differential equation. The modified equation is derived by expanding each term of the dif-
Catalogue number: ACDF
ference scheme in a Taylor series and then eliminating time and mixed derivatives of a higher order than the order of the only
Program obtainable from: CPC Program Library, Queen’s Urnversity of Belfast, N. Ireland (see application form in this issue)
time derivative which may occur in the differential equation. From the modified equation, the program determines the order of accuracy and consistency of the difference scheme. For the
Computer: ICL 4-72; Installation; University Regional Cornputing Centre, Praha, Czechoslovakia
stability analysis the program uses two methods described in ref. [1]: the modified equation (Hirt) method and the von Neumann method.
Operating system: Multijob
No. of bits in a word: 32
Restriction on the complexity of the problem Generally, this very difficult problem is solved only for simple partial differential equations (as specified below), for which an analytically soluble method exists. For computing the modified equation, the differential equation must be in two independent (time and space) variables with one unknown (looked for) function of these variables and must contain only one derivative of this function with respect to the time variable as a linear
Peripheral used: disc drive
term. The time order of this derivative must be one and its total order must be the minimum order of all the derivatives appear-
Programming language used: REDUCE 2 High speed storage required: depends on complexity of input parameters, minimum 350 kbytes
No. of lines in combinedprogram and test deck: 771 Keywords: partial differential equation, difference scheme, modified equation, order of accuracy, consistency and stability of a difference scheme Nature ofphysical problem One of the main problems in solving partial differential equations using finite difference methods is the evaluation of vanous properties of a finite difference algorithm such as order of accuracy, consistency and particularly stability. These problems must almost always be solved by hand. Using the analytical programming language REDUCE 2 [2], the program determines the above properties of finite difference algorithms for simple partial differential equations. Method of solution For the solution of the problem the method of modified
ing in the differential equation. An exact stability analysis using both methods is performed in case the differential equation contains one time derivative of the first order, the differential equation is linear with constant coefficients and the difference scheme is two-level (in the time step). An heuristic stability analysis using the modified equation method can be made without the restrictions given above. In addition, the exact von Neumann stability analysis can be made separately (using only the procedure NEUMANN) for all two-level difference schemes approximating linear partial differential equations with constant coefficients. Typical running time Interpretation of the complete program MODST (without execution) for the ICL 4-72 with 350 kbytes takes 75 s (we use the CHKPOINT-RESTORE facility of the LISP language), while with 450 kbytes it takes only 69 s. Execution of test run no. 1 with 400 kbytes takes 110 s, run no. 2 with 450 kbytes takes 45 a, and run no. 3 (using only the procedure NEUMANN) with
OO1O-4655/84/$03.OO © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
R. Liska
176
/
Accuracy analysis offinite difference methods
400 kbytes takes 15 s. The running times depend on the complexity of input and on the given storage. Unusual features of the program The program can be used in an analogous way to the test runs without detailed knowledge of the REDUCE 2 language.
References [1] R.F. Warming and B,J. Hyetl. J. Comput. Phys. 14 (1974) 159. [2] A.C. Heam, REDUCE 2 User’s Manual (University of Utah, Salt Lake City, 1973).
LONG WRITE-UP 1. Introduction
2. Theory
The program MODST is designed for the stability and accuracy analysis of finite difference approximations to simple partial differential equations for initial value problems. The method used is the technique described in ref. [1] using the so called modified equations and von Neumann stability analysis. The modified equation is the actual partial differential equation, which is solved numerically, apart from round-off errors, by the application of a given difference scheme to solving an initial value problem. The modified equation is derived by expanding each term of a finite difference equation in a Taylor series and then eliminating the time, and mixed time and space, derivatives of a higher order than the order of the only either time or mixed derivative, which may occur in the original partial differential equation.
2.1. Modified equation
The modified equation contains infinitely many terms, but only the first few terms of lowest order have to be computed, so that we can use a truncated Taylor series. The terms appearing in the modified equation which are not in the original partial differential equation represent a type of truncation error. These error terms allow us to determine the order of accuracy and consistency of a finite-difference scheme. The stability conditions may be determined from the full modified equation. Hirt was the first to introduce this idea in ref. [3] and thus the modified equation stability analysis method is sometimes called the Hirt method. Von Neumann stability analysis is the second method for stability analysis used in the program. Section 2 gives the basis of the theory; some details are omitted and most of them can be found in ref. [1], from which almost all the theory is taken.
In this section the partial differential equation a~lu/~taxx + Mji~) 0 (2.1) =
is considered, where M~(U) represents a spatial differential operator containing derivatives of higher order than ix and ii is the dependent variable, which is a function of a spatial variable x and a temporal variable t. To calculate a solution of the differential equation (2.1) by a finite difference method, a grid of mesh points is introduced in x—t space with increments LIx and ~Xt. The partial differential equation (2.1) is replaced by a finitedifference analogue D(T0, T1, ~X, L~t)u7 0, (2.2) whose solution on the discrete mesh will be denoted by u’ u(jtXx, nL~t). In order to use a Taylor expansion in deriving the modified equation, we assume the existence of a continuously differentiable function u(x, t), coinciding at the mesh points with the exact solution of the difference equation (2.2) (this function differs from the solution it of the differential equation (2.1)). 7~) and T1 are shift operators =
=
7~ u5’
=
u~+1,
T~ u
=
~
(2.3)
D is the difference operator. Eq. (2.2) is called the
difference scheme. By applying such a scheme at each mesh point, it is possible to compute the solution u7. The main question with finite-difference methods is how near the approximate solution u7 calculated from the difference scheme (2.2) is to the exact solution it(j~x, nL~t)of the differential
R. Liska
/ Accuracy analysis offinite difference methods
equation (2.1). The convergence of the approximate solution to the exact one is determined by certain properties of the difference scheme, such as the order of accuracy, consistency and stability, We can compute a differential equation that is actually solved by the difference scheme and from this, the so called modified equation, we can then specify such properties of the difference scheme. The modified equation is derived by a two-step procedure outlined as follows. The first step is to expand each term of a given difference algorithm (2.2) in a Taylor series about the point (jL~x, n~t). By this expansion, from the difference scheme (2.2) a partial differential equation + ‘ U/at
ax x +
~c, ( u)
=
0
(2.4)
is obtained, including an infinite number of both space and time derivatives. Nxt represents a time—space differential operator, which must not or aixu/axix, include any derivative a~~1u/atax~ where jx ~ ix; otherwise it means that the difference scheme (2.2) is not a finite-difference approximation of the differential equation (2.1) and we can finish our analysis by saying: the difference scheme is completely wrong. The second step in deriving the modified equation is to eliminate all time derivatives from N~~( u) in (2.4). Considering what has been written above, every time derivative appearing in the differential operator N 5~can be expressed as some derivative of eq. (2.4) and in such a way it can be eliminated from operator N~, by substitution. The time derivatives must be eliminated in the following sequence 2u/at2ax’~, a’~2u/atax”~, a”~ at~3u/at3ax~,a’~3u/at2ax~1 (2.5) Every derivative in this sequence, when expressed using (2.4), does not contain any derivative appearing earlier in this sequence. Considering infinite eliminating, we get the modified equation
a~~lu,/atax1x + Mx(u)
=
TER,~(u,Z~x,st). (2.6)
177
round-off error the modified equation represents the exact partial differential equation that is solved by the finite-difference algorithm; various properties of a finite-difference scheme can be deduced quite simply by examining a truncated version of the modified equation. For the analysis of a scheme only some low order terms (according to the orders of derivatives of u and degree of powers of the increments L~xand z~t) of TERx in (2.6) are required; thus in practice it is sufficient to use a truncated Taylor series and finite elimination (2.5) in deriving the modified equation. 2.2. Order of accuracy and consistency The order of accuracy of a finite-difference scheme is the lowest power of either the increments i~xor ~t appearing in the error terms TERX of the modified equation (2.6). A finite-difference scheme (2.2) is consistent with a partial differential equation (2.1) if TERX tends to zero as E~xand ~t approach zero in an arbitrary manner. 2.3. Modified equation stability analysis This section considers only linear partial differential (2.1) with constant coefficients, where ix equations = 0, and only two-level (in the time step) difference schemes (2.2). For such cases the equivalence between the modified equation and the von Neumann stability analysis methods is shown in ref. [1]. Without these restrictions it is possible to use the modified equation stability analysis method as an heuristic. In refs. [3,4] the validity of such an heuristic approach is shown by practical examples. Now the modified equation (2.6) with the above restrictions can be rewritten as
au/at= ~ ,.&(p)a’°u/ax” =
(2.7)
0
and the difference scheme in the form mr
nr
The spatial differential operator TER 5 on the right side of the modified equation does not appear in the original differential equation (2.1) and constitutes a form of truncation error, introduced by the finite-difference analogue. Apart from the
~
q= —m1
BqUJ~+~ =
~
q=
—n1
A qU7+ q’
(2.8)
Let ~s(2r) be the lowest non-zero coefficient of an even order spatial derivative appearing in the modified equation (2.7). Then a necessary condition
178
R. Liska
/ Accuracy analysis
for the stability of a finite-difference scheme is (—1Y’~s(2r)>O.
offinite difference methods
(2.8) it follows that
(2.9) I
n~
~ Denote =
If s If s
max(m1
+
mr, n1
+ nr)
—
r.
(2.10)
0, then condition (2.9) is already sufficient. 1 and if we denote 4)((~x2/3)~t(2) I (8~t/i~x —~t~t2(2)—~(4))forr=1,
I X
=
2
A~exp(iq0~
-2
m.
~
Iq=-m
exp(iq9)~
Bq
(2.15)
,
=
I
(_4Y+l~t/(2~x(2r+2))
((3
+
(2.11)
r)~x2~(2r)/12 forr~2
where 9 kz~x.Using the trigonometric identity sinxsiny=cos(x—y)—cosxcosy, (2.16) =
we get I
2 Aq~~~(~0)
~
=
AqAp
~
q-nip—n
q=—n
1
then the condition v> 0
x cos(q — p)ø. (2.12)
is a sufficient condition for the stability of the finite-difference scheme. Case s> I generally indicates that more spatial points are used in the difference algorithm than is necessary for the order of accuracy achieved. 2.4. Von Neumann stability analysis This section considers a linear partial differential equation (2.1) with constant coefficients and with a general spatial operator M~and a two-level difference scheme (2.8). The von Neumann stability analysis consists in replacing each term u1~ ~f the difference equation (2.8) by the kth Fourier component of a harmonic decomposition of u~,i.e. by v~(k)exp(ikjz~x), where v~(k)denotes the kth Fourier coefficient, and in subsequent examination of the amplification factor 1(k)/v~(k), (2.13) g(k) derived v’~ from (2.8). The amplification factor g(k) =
(2.17)
Using cos 9
1
=
cos nO
=
—
2z,
(2.18)
~ D.zi,
D0
1, D,
=
(—
=
1)’2(2~1),
0
2(O/2), (2.17) gives where z ~
sin
=
2
Aq exp(iq0)
in
I
~ ~
=
(2.19)
J°
q= ~
where m
=
max(n 1
+ nr,
m1
+
me). In the same
way 2
~
B~exp(iqO)~
=
q= —,n~
~ /31z’. J=o
(2.20)
Hence from (2.15), (2.19) and (2.20) it follows that g(k)~2 1 + ~m (y~, (2.21) =
—
/~ ~‘.
/
j=0
is generally complex. Since the two-level algorithm has only one dependent variable, a necessary and sufficient condition for stability is
Let r be an integer, so that
~ 1, (2.14) for all values of k. Now let us examine the value of g(k). From
This r is the same as the r in section 2.3. (2.21) can then be rewritten as
=
for 0
g(k)J2
=
1
~j
—
<
r
and
zrS(z)/P(z),
Yr *
$r’
(2.22)
(2.23)
R. Liska
/ Accuracy analysis offinite difference methods
where
left sides of substitution equations before performing substitutions.
S(z)= ~ ct~zJ, a
0=S(0)*0
(2.24)
=
2 P(z)
=
179
~0$~z’ =
~Bq
exp(iqO)~ > 0.
Again, the s here and in section 2.3 are the same. As 1 ~ z ~ 0 and P( z)> 0, the necessary and sufficient condition for stability (2.14) can be expressed in the form s(0) > 0, S(z) ~ 0 for 0 <.~ ~ 1, (2.25) which is the quite simple analytic form of the stability conditions already looked for.
3. Program description The program can be used in two ways: 1) The use of the main procedure MODSTAB to calculate the modified equation, to determine the consistency and order of accuracy of the given finitedifference scheme and to make its stability analysis. 2) The use of the procedure NEUMANN to make a stability analysis by the von Neumann method. Restrictions on the input differential equations and difference schemes have been given in the program summary; the form of parameters of both these procedures will be shown in this part in sections describing separate procedures. 3.1. Symbolic mode procedures A new statement REMWEIGHT is performed. Using this statement with any argument, for exampie REMWEIGHT T; causes the removal of all declarations of the weight, which have appeared in statements WEIGHT in the program earlier. Thus the system no longer executes the neglecting of terms given by these statements. The clearing of the weight is done by the clearing of the global variable of the REDUCE system WTL! *, which contains declarations of the weight. Symbolic procedure SIMPSUBE defines a new operator SUBE, which is similar to the REDUCE operator SUB,but differs from it by evaluating the
3.2. Necessary declarations of variables, operators and arrays The top-level variables ORD, S and MODE are declared. These variables have to be top-level, because they are used simultaneously in more than one procedure to hold some value (similar to the FORTRAN COMMON variables). All operators used in the program are declared and also all auxiliary procedures have to be declared as operators, because the main procedure where they appeared preceded their definitions in the program. The LISP language functions MIN, MAX, FIXP and EJECT are declared as symbolic operators, in order to allow their use as operators in the algebraic mode. The operators MIN and MAX give the minimum and maximum of their arbitrary number of numeric arguments, respectively. The operator FIXP is a predicate, which is true only if its only argument is an integer. EJECTO; means a new page in the output file if this actual output file is the principal output file LISPOUT. The auxiliary arrays COF, MI, ALGSOL and FACT are declared. Array FACT is filled in by factorials. 3.3. Main procedure The main procedure MODSTAB has eight parameters. DFE is an expression representing the differential equation (2.1), DFE = 0 with two independent variables, i.e. indeterminates X and T, and one dependent variable, i.e. the unknown function U(X, T). Derivatives are represented by the REDUCE operator DF. Restrictions as to which differential equation can be used, have been given in the program summary. FDE is an expression representing the difference scheme (2.2), FDE = 0, which approximates the given differential equation in the grid of mesh points in X—T space with increments DX and DT, respectively. The values of the approximate solution at the mesh points are designated by FU(i, j), and the mesh point about which we want to perform Taylor expansions during the calculation of the modified equation is i =j = 0. DFT is the only either time
R. Liska
180
/ Accuracy analysis offinite difference methods
or mixed derivative appearing in the given differential equation. Parameters N, M and K are integers, which specify the neglect of terms DX’ for i ~ N, DT’ for i ~ M and DX’DT’ for i +j> K during the calculation of the modified equation. PHEUR is the flag, determining the possibility of using the heuristic method of stability analysis. If it differs from the indeterminate HEURISTIC, the procedure executes a stability analysis by the modified equation and von Neumann methods if and only if the differential equation is linear with constant coefficients, DFT = DF(U(X, T), T) and the difference scheme is a two-level one. If the flag PHEUR is the indeterminate HEURISTIC, the procedure executes a stability analysis even if these restrictions do not agree; then in this case it is performed by the heuristic modified equation stability analysis method. MSG is a string containing some information about the given differential equation and difference scheme, which will be printed in the head of the output. Expressions in the arguments have to be written in the usual REDUCE (FORTRAN like) syntax. The notation of the mathematical symbols used in the text of section 2 in the program is as follows:
derivatives of a lower order than the order of DFT, which would mean that the difference scheme used a bad finite-difference approximation to some derivative. The modified equation is written without the terms neglected due to N, M and K, and without the time derivatives and space derivatives of a higher order than ORD, which have been included in the modified equation during the elimination of time derivatives. The procedure examines the truncation error term TER to verify consistency. The necessary and sufficient condition for consistency is that TER does not contain any zero or negative powers of DX or DT. From TER, the procedure determines the orders of accuracy of the given difference scheme in the increments DX and DT. Depending on the given differential equation, the difference scheme and parameter PHEUR as mentioned above, the procedure analyses the stability of the difference scheme. For the von Neumann method a special procedure NEUMANN was devised, while the modified equation method is a part of the procedure MODSTAB. The modified equation stability method was described in section 2.3. Using the operators MAXPOINT and MINPOINT the procedure determines MR, ML,
X
NR and NL. The procedure MINEVENORDERDF determines R. Then the conditions of stability (2.9) and (2.12) are written to the output file. In the heuristic approach only the necessary
-
x
ML
-
m1
—
mr
U(X, T) IX
zIx NL
—
n1
—
z~t NR
—
~r
—
r
—
S.
T
—
t
DX
—
DT R
MR
S
—
u(x, t) ix
M
—
m
TER
—
TER~(u,~x,
-
~t)
Derivatives in the program are represented by the REDUCE operator DF (DF(U(X, T), T, i, X,
j)
=
~
t)/at’ax~).
After initial declarations and printing of the head of output including the given arguments, the procedure verifies if the restrictions given on the differential equation are satisfied (see program summary). In the first step of the modified equation calculation the procedure expands terms FU(i, j) of the difference scheme in a truncated Taylor series up to ORD = MIN(K, MAX(N, M)) order derivative terms. The elimination of time derivatives in the second step proceeds in the sequence (2.5). Then the procedure verifies that the calculated modified equation does not contain mixed
condition (2.9) is written. We have to mention that in the von Neumann stability analysis as well as here, the conditions of stability are not always written in the simplest form. They can be modified to a simpler form by considering that DX > 0, DT> 0, etc. 3.4. Auxiliary procedures The procedure NEUMANN executes the von Neumann method of stability analysis of a difference scheme, described in section 2.4. It is used from the procedure MODSTAB or it can be used alone for any two-level difference scheme, which approximates a linear partial differential equation with constant coefficients and with two independent and one dependent variables. The procedure has two parameters XX and YY, which give a
R. Liska / Accuracy analysis offinite difference methods
difference scheme XX =0 and a differential equation YY = 0 in the same form as parameters FDE and DFE of the procedure MODSTAB. Variables ORD and S indicate whether procedure NEUMANN is used alone (then ORD = 0 and 5 = —1000), or from procedure MODSTAB. First the difference scheme XX is broken into two parts according to the time step as in (2.8), then using procedure EXPSIN polynomial S(z) = MODE is determined; conditions (2.25) have been applied to this polynomial. S( z) is the polynomial in z of degree s. If s ~ 3, then the procedure finds other conditions from (2.25) by examining the polynomial S(z). For s=0 the necessary and sufficient condition of stability (referred to further only as the s-condition) is S(0)> 0. For s = 1 the s-conditions are S(0)> 0, S(1) ~ 0. For s = 2 and a solution of the linear algebraic equation DF(S(z), z) = 0, the s-conditions are S(0)> 0, S(1) s 0, 0
z2 solutions of the quadratic algebraic equation DF( S(z), z) = 0 the s-conditions are S(0)> 0, S(1)~0, 0 1 by using the REDUCE operator COEFF and the procedure ALGEQUATION; some conditions are written by using the procedure CONDWRITE. The procedures RE, IM and NORM2 have one argument, a complex expression. The value of the procedure is its real part, imaginary part and square of its modulus, respectively, The procedure EXPSIN has one argument in the form —*
~CqFU(q,
1),
(3.1)
q
it executes substitutions FU(q, 1) = exp(iqJ), (2.16) and (2.18) and returns
(3.2)
2
~C q
exp(iqO)
=~a ~ J
(3.3)
2(O/2) (see (2.16)— as a polynomial in zis =used sin in the von Neumann (2.19)). The procedure stability analysis.
181
The procedure ALGEQUATION has one argument, an integer either 1 or 2, and it solves either a linear or quadratic algebraic equation, which has its coefficients in the array ALGSOL. The solution of the equation is returned to the array ALGSOL. The procedure CONDWRITE has one argument and is used for repeated special printing. The procedure TAYLOR has two arguments, integers A and B, and calculates truncated Taylor series of a continuously differentiable function U(X, T) at the point (A * DX, B * DT) about the point (0, 0) up the ORD order derivative terms. The procedure MINORDER has two parameters X and DD and determines the minimum degree (possibly zero) of indeterminate DD occurring in the expression X. The procedure MINEVENORDERDF looks for the minimum r> 0 so that DF(U(X, T), X, 2r) appears in its only argument. The procedure returns r as its value and fills array MI by the coefficient of this minimum even order derivative of U(X, T) with respect to X and the next even orderone. 3.5. Operators The operator AD determines the order of a derivative, which is its first argument, with respect to an indeterminate, which is its second argument. The operator DFP has one argument and its value is 1 if this argument is a non-evaluated derivative and 0 in the other cases. The operator FUP has one argument and its value is 1, if this argument is the operator FU with two arguments, and 0 in the other cases. The operator LINEARP has two arguments, the first is any expression and the second should be T. Before using this operator we must declare that X depends on T by the statement DEPEND X, T;. Then the value of the operator is zero if its first argument represents a linear partial differential equation with constant coefficients with two independent variables X and T, and non-zero in the other cases. theThe formoperator SUBONE has one argument in ~
C 1~FU(i,j)
(3.4)
R. Liska
182
/
Accuracy analysis offinite difference methods
and adapts it to the expression ~ FU’
~
guage, i.e. the creation of a new file MODE, which (3
~\
I.]
whose form is needed for the argument of the following operators. The operator TWOSTEP has one argument in the form of (3.5) and if the value of j in this argument is only 0 or 1 then the value of the operator is 1, in the other cases 0. The operator verifies whether the difference scheme is a two-level one, The operators MINPOINT and MAXPOINT have two arguments, the first in the form of (3.5) and the second is an integer. They determine, which mesh points are used in (3.5), i.e. integers NR, NL, MR and ML. 3.6. Input and output of the program As input the program requires only a call of the procedure MODSTAB or NEUMANN with the corresponding arguments described above. Because REDUCE is an interpretative language, this call can be made only after interpreting all program MODST definitions (lines 1—738 of the program). The output of the program would be directed to the file LISPOUT and will contain information from the analysis of the given difference scheme.
4. Test runs Test run no. 1 must be the first of the three. It shows the CHKPOINT facility of the LISP lan-
contains the necessary part file of the REDUCE tem and our program. This occupies aboutsys50 kbytes on the disc drive. The other two test runs show the use of the file MODE. This access is suitable for the repeated use of the program, because interpreting all the REDUCE program MODST requires 69 s with storage 450 kbytes, while reading the file MODE into the storage by the use of the LISP function RESTORE requires 1.8 s only. For single use lines 740—742 of the program can be omitted and line 743 replaced by calling procedure MODSTAB or NEUMANN, without using the file INFLE. We have to mention that the procedure NEUMANN does not need the whole program and for using only this procedure the program can be modified to a smaller one (all requirements of the procedure NEUMANN are written in the comment to this procedure). Test run no. 2 shows the heuristic stability analysis and no. 3 shows the separate use of the procedure NEUMANN.
References [1] R.F. Warming and B.J. Hyett, J. Comput. Phys. 14 (1974) 159. [2] AC. Heam, REDUCE 2 User’s Manual (University of Utah, Salt Lake City, 1973). [3] C.W. Hirt, J. Comput. Phys. 2 (1968) 339 [4] L.D. Cloutman and L.W. Fullerton, Report LA-6885-MS, University of California, Los Alamos (1977).
R.
Liska
/
Accuracy analysis offinite difference methods
183
* (U .5-
a C) 3
+ • • P
5 3 4
4 P 4
I-. 5
3 * *
x
* * • • • • 3 3 • 3 $ 3 S 3 • 5 •
* * 3 S 3 4 * 3 3 4 • .5 * 4
* 4 3 * * 4 5 3 3 4
C.,
$ 3 .5 4 P * 5 3 * 3 3 * 3 3 5 * 3 3 4 4 S *
P..
In ~ 0 a P.’ ‘4
—
4 4 3 * 3 3 3 5 * 5 3 5 3 4 3 * 3
3 * 3 5 .5 4 5 * P 3 5 * 5 4 3 3 S P 5 3 3
a
•
•
* .5 * * 5 * .5 3 P
0 UI _l ‘4 IZ
•
UP
• 3 5 3 • * 3 • S * S * 3 4
& UI U. I~
5-’
p
In
•
0
S
UI
PC
•
~,.
II.
.4
NP
0
MI
U. U) 0 ‘4 ‘4
PC
$
0
•
1
9P
MI
us U. U..
U. U.
U. 0
5 S 5
p.-
lIP
~
)“ —1
UP
‘p
~.
I— ~ 0
I—
‘4 PC ‘4 U. 0 U.
* 3 * *
0 MI PC CD .~.
‘p
up
•
S * * 4
$ 4 4 3 * * 4 5 3 4 * 4 4 * 5 3 5 * 5 p 5 S 4 5 S 4 3 5 4 * S 4
w
0
UI ~..
w PC
I..)
U) MI IJ ~C MI ‘4 U) U. U.
S
5 5 .5
3
* ‘P * 3 *
$ 4 .5 S S S
j
MI C.) PC.
0
5
~
PC MI PC (a U)
UI U.
• •
U. ~ — In U. 0
• •
~
..C
3 3 S • 5 3 5 5 3 .5 S 3 S S S 5 5 5 *
—
a ‘4 a I— U. ‘4 U.
4
5
‘P .5 4 5 S 3 5 5 5 * * * $ S $ 3 4 4 3 3 5 5 * * 5 5 5 4
.5.
3 —
C.) 4
a . ~—
I
5—. U.
I -~ Cv
5-.
(U ‘C 0
$ 3 4 3 5 3 5 5 $ 4 3 3 3 5 4 5 * 5 3 * 5 P 4 3 S 3 3 .5 5 5 ‘P 3
I (U 5—
a 3 (U 5 0
-
a
U.
+ I—. 0 5i’s
U. U.
0 ‘4 0 PC PU ~ I
PC
-
“~
0
4 U.
4
‘4 —á I — • 0 PC PC ~ U.
p.~ a~ MI p.”
* 4
& 0
U
* S
U.
(U
5
~ 5
—
•N
5 .5 $ S
PC 0 — i-
U 5 N
$ 3 5 * .5 * S $ * .5 4 * S 5 * 5 * 4 $ 5 *
~ 0 UI
S
VP 5 A
a
3
0
PC a 3 I’0
U. U.
0
N
_I ‘4 .‘. IPC UI U.
UI
-
p.-
s
N N N N N N N N N N N N N N N
‘4
I—’
p
~
+
~a U)
PC
U) Ca PC
.S
U) U.
P.’ PC
~
‘ U. 0
UI U. U.
a a
nN 5 N N N N N N N I N N N I
U
•
0
‘P S 3 P.-
PC
0
p
S IJ
UI ..J W
-P 5 3
PC
UI
PC
a
PC
0
* 0 .. 5’
5 CI
P 0 U.
U.
(U
P.a *
4 5 * 5 3 * 4 * S * 3 3 * 3 3 4 * * 5 5 4 3 5 3 5 3 4 4 5 5 5 .5 3 4 3 S 5 3 5 * * 4 S 5 3
(U (a 5-
3 C) * —, PiP PC — P.. ...
PC ~
0 — ‘4
II.
0 4
0 PM
.
UI —
-~ (U
PC 0 ~
(U
p.. a
.J ‘4 ‘P
I’PC UP ‘4 UI U. U.
a
4 Cv 1.1
PC I-. ~-.
* I... 0 .5 (U
5
S
PC a I
0 0 PC
P~ * (a
N
(U
* * * S 5 * 5
a
‘P
—
.
U) ‘4
0
I ‘0
*
4
I
a
CI
4 4 PC
0
~
I.
I-
5 4 P” 0 ‘P
S..
U..
0
“~
I
UI PC UI
-
PC -~
5 5 ‘P 5 3 5 ‘P 5 3 5 S 3 $ 5 * S 5 5 3 S 3
PC.
0 3 (‘U
3 5 P
—
(a
4 ‘P •
PC
PC 0 —. I” ‘4 ~ 0 UI 0 UI IC
I 5 N U N N N N N N N II II II N N U N N
P..
~
p.-
U.
a
pp
UI U.
5 $
U.
(a ‘4 U. ~
5 * 5 $ * 5 * S * * S 5
PM PC UI PC ~ U)
CI CJ
* ‘P
‘4
UI (a PC UI U.
Cl) U.
3 ‘P * * * 3 S S ‘P $ *
PC UI p... lIP a U) PC 0 Ca UP
PC ~
S P S * 4 3 3 4 3 3 * 3 * * 3 4 4 3 4 4 * 3 4 * * 4 3 5 * 3 3 $ 3 S 3 3 3 .5 S 3 54 0 4 P * * (U 3 PC * 0 5 3 UI 4 ‘4 4 ‘4 5 5 UI 3 PC S MI S PC S p,~ 5 4 * ~I S (a 5 PC 3 UI * U. S
a
0 U. 0 ~.-
U. 0
U I-’ is
P.” -
x
U) U.
u. a
0 U. 0
PU
0 ‘4 0
S
UI PC
p.-
* 4 S * *
R. Liska
184
..J
‘4 PC ‘4
•
P’ I” a
N N N ft N N U
-a a In
‘4 I— U)
I-’ ‘4 fl
I
UI
I P
‘.
PC 5-
PC 0 .
N ft
a
N U N N II N N N
0 PU
pp
SI U.
0 0 PC
p
I I I I I I I I
* ‘P * * 5
S .5 * 5 5
U) PC 0
* * $ * *
$ * 5 S
* ‘P
* **
* *
— I-’ ‘4
S $ *
* 5P
S 5 ** $ *
‘p
5
‘p
5 *
‘4 U.
S * 5 * S S S S * S $ * S * * * * * * * $ 5 * * $ *
.4 ‘4
* 5 5
UI PC UI PC (a U)
‘p ‘p
5 $ 5 *
MI ‘-5 ‘4
* * $
.5* .5 4
—
‘p
5-. PC U)
0 U. 0
a
P” I’-
$ * * * * *
..J — U) ‘4
I
PC
0 o
CP
S
PC
I
—
P.—
0 PC 0 Ca
0 A
I— PC
U) ~
U) a
—“s
(‘U
.SI IC U.
~ U)
+ (SI
UP
0 * (SI ~
,-
U. ‘4 U)
I
U)
S
UI
S p.0 S
4 * * * 4 * *
‘p
0 II,
* * *
U) ~p
‘p
‘p
up
* ‘P 5 S
$ 5 $
LI
PU PC (U
5 * 5 * S ‘p
(a 5-
S $
LI UI PC
0 A
U)
PC PC 0 ~
$ S $ $ $ * 5
I—’
a’‘4 ‘4 U)
up
S $ ‘p $
U)
$ * $
.,j • • •
‘p
$ * S
I’.’
*
CD
(5. —. III
I
I-’
PC
0’
~‘
C> ,..
II.
5’
— ‘4
* * 5
$
~
‘p
_., III .-4 PC PC
a PC U)
I-.
*
. (.5’
$ * *
0
0. pp — ...I
‘4 PU 5-~ PC PSI
up PC
—.
‘4 p.’
(“1
0
* $ S $ 5 * S 5 ‘p
C”
~
* *
—
P.’
* $
PIP a
‘4
*
‘4
U) SC PS. .‘. 0
‘p *
~ PU
S
PC
‘p *
$ 5
II.
‘p ‘p
U.
0
(5.
~ I—
~) a
* ‘p ‘p ‘p
S
— (5) P” —I
‘4 PC ‘4
I-’ ‘4
Ca ~.
0
4
$
~
$
— (5) UI 0
$ ‘p $
PC Ca ‘4 Ca U. 0. 0, ‘4
~
PU PC
4 4
5 * 5 ‘p ‘P $ $ S P 4 * $ * 5
PC P55
~ *
‘p
LI
‘4
0
•
PC PU
C)
PC
*
pj
(5)
(51
p
* $ $
I— a *
PC
•
0 PC 0
(SI
U) U) L( MI
p-.
•
‘P $ 4 ‘P S
LI
*
a
+
0
‘~
4 *
PC 0
PC W PC “
5
PC 0
* * S * S * *
1U
‘4
* $
‘p 5
4
II,
‘p
..a ‘4
(I)
‘p
U)
0
* ‘p
S 5
‘4 0
‘4
‘4
is
a PC ‘4
‘p
‘4 Ui
‘4 ‘4
Cv
S S
U)
U,
Ui
(F)
a
PC
*
0
U. ~
PU -)
~ PC
5
PC
* * S * 5
0
‘4
MI
* * *
0
P.
* 5 .5 *
‘p
SI
5 5
-3
5 *
5-’ ‘4 ‘4 0.
I—
‘P $ $ * * $ * $ *
* 5 $
‘4
$
— (.5 .. IL
$ $
$ $ *
0
,0
* * .5 5 * *
p.PC
S * $ $ *
‘p
(.5
‘p
$
~
0
I I I I I I I I
5 .5 $ * S
* 5
‘4 ‘4 U)
0
PC .4 PC ~
-.5
,.,
I-’
PC
$ * * 5 5 S $
* .5 * *
—
I I I I I I
‘p
“ PC UI ‘4 UI U. U.
PC ~
“.
5-’ Us PC
a U.
0
PC 0
‘p P.’
PC ‘4
5 S *
$
* 5 5 $ S * $ * $
‘p
p... up
* $ $ S S 5 S 5 * * * * ‘P * $ * * 5 S $
(a U. ~ (5)
s
U) PC 0
Ca
p
U. a
‘~
S $ p $ * * • * * * 5 * * $ 5
PU 5PC MI U. UI II,
•-.
‘4 0 U.
I— ‘4
* $ * • *
* $ **
‘p 5
up
I I ( I I I I I 5 I I I I I I I I I I
S .5 * ‘p * $
$ 5 * ‘P *
~ 0
UI ‘4 5-,
0 0 PC I— W PC
* • $ * * *
‘p
* * *
*
P. P.—
U)
* * * * S *
* S ‘p
**
U, 0
)-
5 S * * $ $
UI
O
N N N I N I N N N N N 5 N N I II
* $ S * $ * S S * $
* 5
Us C) PC PU ‘4 Mi II. II,
—‘
* S $ * * S S S $ ‘P
UI
MI PC UI PC (a NP
C’)
/ A ccuraey analisis offinite difference methods
5 5 *
5 ‘P $
C) Uj
.5 S
‘P
$ S
I
0, PC
‘p
$ 4
‘4 UP IU
5 S *
S * $ $ 5 $ ‘p
I (5) •
0 PC
* $
$ 4 $ $ * $
$ 5 $ *
PC ~
‘p
‘4
$ *
5-
$ * $ 5 ‘I’ *
4 * $
U) PU I”
4 ‘p $
‘p
$ 5
*
5
.5
5
R. Liska
/ Accuracy analysis offinite difference methods
185
I-’
5-
II.
PC
+
5U.
is
0
((I PC
I
a .5
50 5’ 515
PC O .5 I— 0 * 3) * Cs) I 5-
*
is
0 -
0 U..
.5 is
PC 0 S I-’ O 55is
I 5-
* is 5• 0 5U.
+ is (5)
PC 0 * 5-0
o
PC 0
PC 0 — ~—
“
‘4 ~
a U) .4
.4 — p.PC PU ‘4 U) U. II., — 0
N N N U Cl ft 5 U ft II N II N U ft ft ft
• p..’ •
P.0 *
is
(.5 5-
PC
UP +
Ui PC UI 14 U)
is
(U •
PC is
I— • PC 5-
~ 555.
0 * (a
‘4 0 U. “’p
is
* * I— 0 *
U) Ca
PC U) U. PU U. IS. ,_.
0
N N U N ft U I II N N N II (I II N N N
0
SI
pp
II
• — 5’
..a CD
‘P ~ 5
MI
PC
*
PC
0
5 * * *
•‘
5
U..
+ is
PC 0 $ I”. 0 5’ 55” ‘ 5U.
(51 pp C 5is
Ca $ is
0
• is ,.‘
I
U) ‘4 U) 0 ‘4 1’)
PC
* $ PC 0
‘4’ 5 5 P.C
SI
5-’ C)
•
PC 0 5 P.’ 0 * ‘5
I 5~
* (a S
* * * * * 3 S * * * * * * * *
* * * * * * * * $
is (I’) •
PC • -5
5-. •
.4 5I U.
0 +
‘0 5is
Cv
PC 0
PC 0 00
PC
5‘4 0 UI
* 5-. 0 *
_I ‘4
(a
* pp
a
I-’
5-. PC UI ‘4 UI
$
U. IS,
+
— 0
(5) 14
* ((4 I
5-. ‘4 ~ 0 Ui 0 U) — U. —
0 0 PC
PC
LI U)
UI PC U)
UI
UI U. U. —
.4’
I-’ U)
•
.-.
.4
U) PC 0
U.
Ca
P.’
• is ~.
a
MI PC Ui
‘4
is
PC 0
PU ~ ‘4
(.5
I— PC
$
ft II U N II ft II N N ft N II N N II It N ft ft
PC 0
—
K
*
•.
5-’ 0
PC P.-
5’
‘p
Ui
0
~
* * $ * * * * S S * * * * * *
5-
* * S * S $ * * .5 * * $ 5 * $ * $ $ * $ * * $ * S * S * S $ * * * $ * * * * * * * * * .5 5 * * $ 5 * * * $ * P * * $ * * 5 * $ * S * 5 5 * 5 * ‘p $ * $ * 5 * $ $ *
*
• 5’ ~ 5’ U.
(5*
I
is PC
ft PU II U
‘0 5 A -) *
5— 5is
U
Cv
* ‘P 3 * $ * * * * $ * 5 * * 3
LI
PC
$
n
is
—
~
(%5
•-
•
U)
U. 0
PC
PC Us PC
N
0 5-
5’
•
0
‘4 ‘4 ~P (a (4 ‘4
-‘ U)
U.
PC
U)
U)
•
14
‘4 U) 0 ‘4 C)
0
is
P.is
5-. •
PC 5’ U.
0
PC U) ‘4 Ui U. U. — 0
UI
PC P.
$ $ 5 S * S * * * * $ * * $ * * * * 5 * * 5 * * S 5 * 5 $ 5 * * * 5 * * $ * * $ * S $ S S * * * * * $ ‘P 5 .5 * * * S 5 S $ $ $ $ 5 $ $ $ * * $ $ * 5 $ ‘p
$
~ — UI PC III PC CI
5) U) C)
PC UI ‘4 (II U. U. PSI
0 U. 0 P 5a _I
— In ‘4 5-’ U) ‘4 0 IC
PC 0 5a
0 PC 0 Ca )_
U. .4 U) U) U) (a
U) PC 14 I’S
I-’ Cl)
0
— ‘4
A
~ MI PC
(.4
$ $ $ * * 5 * 3 .5 * $ * * 5 * 5 * * $ $ * * * $ * * * * $ * * S * * $ * * 5 $ * * * $ $ $ * * 5 * * ‘p $ * * * $ 5 $ $ * $ ‘P ‘p * * ‘p $ * * * * * * * * $ 5 * * $ * * * S
I
* *
up
U)
.4 ‘4 P
p •
5). (I) ISI
.4 CD PC a
‘4
U) PC -3 0 ‘0
‘4 P.‘4 0
UP
‘P
55.
0
UI
OP.” PC PC
.~ PSI I’s
UI PU
P. PC I—
C) PC U) *
S *
R, Liska
186
/ Aciuracv ana4’sis offinite difference
methods
* 0
C’ U.
* I0 * PSI
5 $ $
p
* ‘p * $ * * * * ‘p * S $ * * $ S * $ .5 $ * $ * * $ $ * 5 * .5 * * * ‘p * * S 5 * * $ 5 5 * ‘p 5
•
4 IS ¼* U) PC PC is 0 0 is * 5- (U (‘.1 is 5-is PC 5-C) • is•is is ~
5-. 0 I 5’• 5 U. C’ PC 5- C’ S 5’ C’ is U., ICC> C I • $ isC) CD 0~ — • C’ U) s”U. 1 5- * -5 (‘4 PC U. I •
• •
I”.i ‘4 t 5”
C’ I—
5- ~“
‘CD
5 C’ I— (4
5-’ U)
+ s’. PC I’-.-- 5’ C) C’ C’
C’ P.
5- U. 5is 5’ U. is 5 C) (C> is * • (P (.5 0. — + I~)is
‘4 ‘4 C’)
C)
CD + 5U. (U •
‘4
0,
“
is s”
5(‘4 U)
.4
(4 C) 55) ‘4
• 5- CD CD’-’ PC CD ‘4 IL CD C’
P— I-’ C) • S PC (5) 5* C” * 5’ (-(PU 5- 0.
‘p 5 * $ S 5 $ * $ * $ * ** * 5 $ $ * $ * 5 $ $ ‘p * ‘p $ * $ 5 $ ‘p
* * * $ * $ ‘p S $ ‘p * * ‘p 5 5 $ $
US
4
PC
$
‘P $ *
is
0 • 5— 5-
0
C’ U. * PC * P._ C *
N
(U
PC C $ I-’ C S (V
C)
S is
0
5. is is —‘ •
s’ I 5’ 5-
0
C’
‘.‘
C’ U. S
U.
$ PC 0 * 5-
is
(“.1
PC 0 $ (U
C
* (4
+ I . is
is
0
C)
C>
•
•
C)
0
5’
5’
C’
C’
IS.
N
(U (U
,_,
PC I’(‘p $
•.
W U) (.I ~‘
W ~
PC U) ~ Ui U. IL 5’ C> IS.
0 (~ — ~j)
P .4
‘4 C’ ~ P (_ — .4 — 0
‘4 I— (~
ft ft U (( It
•
$ (Si
I—
I * + (‘.1 is
U’ CD — .5 is.
C) IT fl (I It ii I I II Ii (I It II II I fl (( (( II II II II (I II II ft II (I II II H (~
• C”
C’ U.
$ (‘J
“4
P.0 *
5-.
C
rJ
• CD 5-
•. PC I’C) — .5
+
.4
CD —•
‘4
IS.
I’..
p.. PC Ui ‘4
Ui
(4
4
(5. (5.
I
CD
0
is
— VP
CC
*
is ,—
is 5” I’ .‘. “-
I
‘4
C.
0.
(U
$
5), ‘4 = (4
I
(‘~ (4
I 5-
CD — U)
I’ PC 5-
I
$ I— $
I—
—
U.
p.-
• .4 • is
CD
CD
C)
(5)
‘5. 5”
(1. (5)
is
PU p.‘4 P.
CD $
U) S
((P
PC C)
‘-‘ 5-
• .4 5U. C’ 0. 5-
CD C) PU
CD * 50. $ -4’
$ * $ * * S * * * * * * * .5 $ S ‘p * * $ $ $ 5 * $ * S S $ * * $ $ 5 * $ * $ * * 5 * 5 * * *
V ft (U V
o
U.
C) ‘I’
U) ‘4 ‘4
is PC • is
P” • PC .-
C’ PC
5’ (5. 0.
* C-)
N A
MI PC Us PC
is
.4’
C) (F)
PC
0 *
Us
5-
(4 PC
0 5’ 5is CD S~ ((P
U)
‘4 U) IL U. 5’
*
0
(~
0
0
5
S
IS) a’-
S 5 * * $ * $ * 5 ** $ * $ * *
P.-
+
— .4 I-S 5) “4
is (54
CD a
I—
U) * -3
(5) PS. ‘4
0
(5.1
C’)
PC C)
*
$ S * * * * $ $ * * $ $ * $ $ 5 5 $ * 5 $ $ * $ * * * $ * 5 S 5 * $
PC
U.
I PC 0 *
(‘4
—
5-.
(4
* * *
* $ *
$ $ * S * * $
** * * S ‘p * $ * 5 * $
‘P $ * * * $ * * * ‘p * * $ * * $ * S * * 5 * * * $ * $ * $ $ * $ $ $ * * *
‘4 0
$
*
U,
5-
+
PC
C’
(P
(a
—
(5) * P.’ C ‘p
I—
CC PC I— LU
P. PC CD
‘4 P. C’ PU CD PC
CD ~
I I I I I I I I I I I I I I I I I I
C’ Sn
(54 (
$ ‘4’
C’
I
A
$ * * * ‘p * * * 5 ‘p 4 S * * * $
C PC ‘4 (5’
‘4 ‘4
(5) VP PU C.’ Ui PC
is
.4’
5 I—
* 5 $ ‘p * $ $ $ $ $ * 5
(5
$
S
*
*
(‘.4 5’
(U 5-
S
is
5)
P.5>
o * I-
$ .‘S C-’
I is
CD — (5)
$
.4 0 $ I” C 5’ 5is
CD —
U) * (U
PC 0 $ (U +
•
5PC UI —
$
— U. Ii.,
CD
* * * * $ ‘p 5 ‘p * 5 $ 5 5 $$ $ 5 * S 5 $ * 5 * $ * ‘p ‘p * $ * * $ S $ S 5 * * * $ $ $ *
LU
$
CD —
* 50 $ -3 1 (‘4
(5)
PC CD
(‘4
0 PC ‘4 a’-
U) ‘4 ‘4 U) PU (.1
PC
U)
$ 50. *
a’.’ 0
(.4
—
C’
“4
U) ‘4 .4
I—
+
‘4
P-I C
CD
U)
PC C) (-5
(U
PC 0 -I —~
0. * (U (.1
(.5.
= C) ‘4 ‘5’
PC 0 — — —
* ‘.5
0. CD 0
I
C-’
(‘C
PC
‘4 ‘4 I,’)
I—
C PC (4 — PC
(5) 5
5C) $ -5 (4
I
‘p
14
C’ U)
(U
SI
C)
UP
N A
* 5 * * 5 * $ * 5 $ $ 5 * 5 * $ * * $ $
—
UI — I.) — U(S.
0
* * $ * .5 5 $ * * * * $ * S * * ‘p 5 ‘p $ * * * * * S $ $ $ S *
I— I-’ — ‘4
PC
* * ‘p $ $ 5 ‘p
5-
*
I’
U) —
U) C’ —~ •
‘4 >
0. U’ — ‘.~1
pP
‘4 P.-
~
C’ CD
0 ‘.5 CD ‘0
5’
,—
“4
0.
(5.
‘4 ‘. U) 0 I-’ PC PC U) U)
(C> .4
Q
PU — CC PC
U)
P.’
* $ $
175
PROGRAM FOR STABILITY AND ACCURACY ANALYSIS OF FINITE DIFFERENCE METHODS R. LISKA Faculty of Nuclear Science and Physical Engineering~Technical University of Prague, Bi~ehovà7, 115 19 Praha 1, Czechoslovakia Received 8 February 1984; in revised form 4 July 1984
PROGRAM SUMMARY
equations, described in ref. [1], is applied. The program calculates a modified equation corresponding to a difference scheme
Title ofprogram: MODST
which approximates a partial differential equation. The modified equation is derived by expanding each term of the dif-
Catalogue number: ACDF
ference scheme in a Taylor series and then eliminating time and mixed derivatives of a higher order than the order of the only
Program obtainable from: CPC Program Library, Queen’s Urnversity of Belfast, N. Ireland (see application form in this issue)
time derivative which may occur in the differential equation. From the modified equation, the program determines the order of accuracy and consistency of the difference scheme. For the
Computer: ICL 4-72; Installation; University Regional Cornputing Centre, Praha, Czechoslovakia
stability analysis the program uses two methods described in ref. [1]: the modified equation (Hirt) method and the von Neumann method.
Operating system: Multijob
No. of bits in a word: 32
Restriction on the complexity of the problem Generally, this very difficult problem is solved only for simple partial differential equations (as specified below), for which an analytically soluble method exists. For computing the modified equation, the differential equation must be in two independent (time and space) variables with one unknown (looked for) function of these variables and must contain only one derivative of this function with respect to the time variable as a linear
Peripheral used: disc drive
term. The time order of this derivative must be one and its total order must be the minimum order of all the derivatives appear-
Programming language used: REDUCE 2 High speed storage required: depends on complexity of input parameters, minimum 350 kbytes
No. of lines in combinedprogram and test deck: 771 Keywords: partial differential equation, difference scheme, modified equation, order of accuracy, consistency and stability of a difference scheme Nature ofphysical problem One of the main problems in solving partial differential equations using finite difference methods is the evaluation of vanous properties of a finite difference algorithm such as order of accuracy, consistency and particularly stability. These problems must almost always be solved by hand. Using the analytical programming language REDUCE 2 [2], the program determines the above properties of finite difference algorithms for simple partial differential equations. Method of solution For the solution of the problem the method of modified
ing in the differential equation. An exact stability analysis using both methods is performed in case the differential equation contains one time derivative of the first order, the differential equation is linear with constant coefficients and the difference scheme is two-level (in the time step). An heuristic stability analysis using the modified equation method can be made without the restrictions given above. In addition, the exact von Neumann stability analysis can be made separately (using only the procedure NEUMANN) for all two-level difference schemes approximating linear partial differential equations with constant coefficients. Typical running time Interpretation of the complete program MODST (without execution) for the ICL 4-72 with 350 kbytes takes 75 s (we use the CHKPOINT-RESTORE facility of the LISP language), while with 450 kbytes it takes only 69 s. Execution of test run no. 1 with 400 kbytes takes 110 s, run no. 2 with 450 kbytes takes 45 a, and run no. 3 (using only the procedure NEUMANN) with
OO1O-4655/84/$03.OO © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
R. Liska
176
/
Accuracy analysis offinite difference methods
400 kbytes takes 15 s. The running times depend on the complexity of input and on the given storage. Unusual features of the program The program can be used in an analogous way to the test runs without detailed knowledge of the REDUCE 2 language.
References [1] R.F. Warming and B,J. Hyetl. J. Comput. Phys. 14 (1974) 159. [2] A.C. Heam, REDUCE 2 User’s Manual (University of Utah, Salt Lake City, 1973).
LONG WRITE-UP 1. Introduction
2. Theory
The program MODST is designed for the stability and accuracy analysis of finite difference approximations to simple partial differential equations for initial value problems. The method used is the technique described in ref. [1] using the so called modified equations and von Neumann stability analysis. The modified equation is the actual partial differential equation, which is solved numerically, apart from round-off errors, by the application of a given difference scheme to solving an initial value problem. The modified equation is derived by expanding each term of a finite difference equation in a Taylor series and then eliminating the time, and mixed time and space, derivatives of a higher order than the order of the only either time or mixed derivative, which may occur in the original partial differential equation.
2.1. Modified equation
The modified equation contains infinitely many terms, but only the first few terms of lowest order have to be computed, so that we can use a truncated Taylor series. The terms appearing in the modified equation which are not in the original partial differential equation represent a type of truncation error. These error terms allow us to determine the order of accuracy and consistency of a finite-difference scheme. The stability conditions may be determined from the full modified equation. Hirt was the first to introduce this idea in ref. [3] and thus the modified equation stability analysis method is sometimes called the Hirt method. Von Neumann stability analysis is the second method for stability analysis used in the program. Section 2 gives the basis of the theory; some details are omitted and most of them can be found in ref. [1], from which almost all the theory is taken.
In this section the partial differential equation a~lu/~taxx + Mji~) 0 (2.1) =
is considered, where M~(U) represents a spatial differential operator containing derivatives of higher order than ix and ii is the dependent variable, which is a function of a spatial variable x and a temporal variable t. To calculate a solution of the differential equation (2.1) by a finite difference method, a grid of mesh points is introduced in x—t space with increments LIx and ~Xt. The partial differential equation (2.1) is replaced by a finitedifference analogue D(T0, T1, ~X, L~t)u7 0, (2.2) whose solution on the discrete mesh will be denoted by u’ u(jtXx, nL~t). In order to use a Taylor expansion in deriving the modified equation, we assume the existence of a continuously differentiable function u(x, t), coinciding at the mesh points with the exact solution of the difference equation (2.2) (this function differs from the solution it of the differential equation (2.1)). 7~) and T1 are shift operators =
=
7~ u5’
=
u~+1,
T~ u
=
~
(2.3)
D is the difference operator. Eq. (2.2) is called the
difference scheme. By applying such a scheme at each mesh point, it is possible to compute the solution u7. The main question with finite-difference methods is how near the approximate solution u7 calculated from the difference scheme (2.2) is to the exact solution it(j~x, nL~t)of the differential
R. Liska
/ Accuracy analysis offinite difference methods
equation (2.1). The convergence of the approximate solution to the exact one is determined by certain properties of the difference scheme, such as the order of accuracy, consistency and stability, We can compute a differential equation that is actually solved by the difference scheme and from this, the so called modified equation, we can then specify such properties of the difference scheme. The modified equation is derived by a two-step procedure outlined as follows. The first step is to expand each term of a given difference algorithm (2.2) in a Taylor series about the point (jL~x, n~t). By this expansion, from the difference scheme (2.2) a partial differential equation + ‘ U/at
ax x +
~c, ( u)
=
0
(2.4)
is obtained, including an infinite number of both space and time derivatives. Nxt represents a time—space differential operator, which must not or aixu/axix, include any derivative a~~1u/atax~ where jx ~ ix; otherwise it means that the difference scheme (2.2) is not a finite-difference approximation of the differential equation (2.1) and we can finish our analysis by saying: the difference scheme is completely wrong. The second step in deriving the modified equation is to eliminate all time derivatives from N~~( u) in (2.4). Considering what has been written above, every time derivative appearing in the differential operator N 5~can be expressed as some derivative of eq. (2.4) and in such a way it can be eliminated from operator N~, by substitution. The time derivatives must be eliminated in the following sequence 2u/at2ax’~, a’~2u/atax”~, a”~ at~3u/at3ax~,a’~3u/at2ax~1 (2.5) Every derivative in this sequence, when expressed using (2.4), does not contain any derivative appearing earlier in this sequence. Considering infinite eliminating, we get the modified equation
a~~lu,/atax1x + Mx(u)
=
TER,~(u,Z~x,st). (2.6)
177
round-off error the modified equation represents the exact partial differential equation that is solved by the finite-difference algorithm; various properties of a finite-difference scheme can be deduced quite simply by examining a truncated version of the modified equation. For the analysis of a scheme only some low order terms (according to the orders of derivatives of u and degree of powers of the increments L~xand z~t) of TERx in (2.6) are required; thus in practice it is sufficient to use a truncated Taylor series and finite elimination (2.5) in deriving the modified equation. 2.2. Order of accuracy and consistency The order of accuracy of a finite-difference scheme is the lowest power of either the increments i~xor ~t appearing in the error terms TERX of the modified equation (2.6). A finite-difference scheme (2.2) is consistent with a partial differential equation (2.1) if TERX tends to zero as E~xand ~t approach zero in an arbitrary manner. 2.3. Modified equation stability analysis This section considers only linear partial differential (2.1) with constant coefficients, where ix equations = 0, and only two-level (in the time step) difference schemes (2.2). For such cases the equivalence between the modified equation and the von Neumann stability analysis methods is shown in ref. [1]. Without these restrictions it is possible to use the modified equation stability analysis method as an heuristic. In refs. [3,4] the validity of such an heuristic approach is shown by practical examples. Now the modified equation (2.6) with the above restrictions can be rewritten as
au/at= ~ ,.&(p)a’°u/ax” =
(2.7)
0
and the difference scheme in the form mr
nr
The spatial differential operator TER 5 on the right side of the modified equation does not appear in the original differential equation (2.1) and constitutes a form of truncation error, introduced by the finite-difference analogue. Apart from the
~
q= —m1
BqUJ~+~ =
~
q=
—n1
A qU7+ q’
(2.8)
Let ~s(2r) be the lowest non-zero coefficient of an even order spatial derivative appearing in the modified equation (2.7). Then a necessary condition
178
R. Liska
/ Accuracy analysis
for the stability of a finite-difference scheme is (—1Y’~s(2r)>O.
offinite difference methods
(2.8) it follows that
(2.9) I
n~
~ Denote =
If s If s
max(m1
+
mr, n1
+ nr)
—
r.
(2.10)
0, then condition (2.9) is already sufficient. 1 and if we denote 4)((~x2/3)~t(2) I (8~t/i~x —~t~t2(2)—~(4))forr=1,
I X
=
2
A~exp(iq0~
-2
m.
~
Iq=-m
exp(iq9)~
Bq
(2.15)
,
=
I
(_4Y+l~t/(2~x(2r+2))
((3
+
(2.11)
r)~x2~(2r)/12 forr~2
where 9 kz~x.Using the trigonometric identity sinxsiny=cos(x—y)—cosxcosy, (2.16) =
we get I
2 Aq~~~(~0)
~
=
AqAp
~
q-nip—n
q=—n
1
then the condition v> 0
x cos(q — p)ø. (2.12)
is a sufficient condition for the stability of the finite-difference scheme. Case s> I generally indicates that more spatial points are used in the difference algorithm than is necessary for the order of accuracy achieved. 2.4. Von Neumann stability analysis This section considers a linear partial differential equation (2.1) with constant coefficients and with a general spatial operator M~and a two-level difference scheme (2.8). The von Neumann stability analysis consists in replacing each term u1~ ~f the difference equation (2.8) by the kth Fourier component of a harmonic decomposition of u~,i.e. by v~(k)exp(ikjz~x), where v~(k)denotes the kth Fourier coefficient, and in subsequent examination of the amplification factor 1(k)/v~(k), (2.13) g(k) derived v’~ from (2.8). The amplification factor g(k) =
(2.17)
Using cos 9
1
=
cos nO
=
—
2z,
(2.18)
~ D.zi,
D0
1, D,
=
(—
=
1)’2(2~1),
0
2(O/2), (2.17) gives where z ~
sin
=
2
Aq exp(iq0)
in
I
~ ~
=
(2.19)
J°
q= ~
where m
=
max(n 1
+ nr,
m1
+
me). In the same
way 2
~
B~exp(iqO)~
=
q= —,n~
~ /31z’. J=o
(2.20)
Hence from (2.15), (2.19) and (2.20) it follows that g(k)~2 1 + ~m (y~, (2.21) =
—
/~ ~‘.
/
j=0
is generally complex. Since the two-level algorithm has only one dependent variable, a necessary and sufficient condition for stability is
Let r be an integer, so that
~ 1, (2.14) for all values of k. Now let us examine the value of g(k). From
This r is the same as the r in section 2.3. (2.21) can then be rewritten as
=
for 0
g(k)J2
=
1
~j
—
<
r
and
zrS(z)/P(z),
Yr *
$r’
(2.22)
(2.23)
R. Liska
/ Accuracy analysis offinite difference methods
where
left sides of substitution equations before performing substitutions.
S(z)= ~ ct~zJ, a
0=S(0)*0
(2.24)
=
2 P(z)
=
179
~0$~z’ =
~Bq
exp(iqO)~ > 0.
Again, the s here and in section 2.3 are the same. As 1 ~ z ~ 0 and P( z)> 0, the necessary and sufficient condition for stability (2.14) can be expressed in the form s(0) > 0, S(z) ~ 0 for 0 <.~ ~ 1, (2.25) which is the quite simple analytic form of the stability conditions already looked for.
3. Program description The program can be used in two ways: 1) The use of the main procedure MODSTAB to calculate the modified equation, to determine the consistency and order of accuracy of the given finitedifference scheme and to make its stability analysis. 2) The use of the procedure NEUMANN to make a stability analysis by the von Neumann method. Restrictions on the input differential equations and difference schemes have been given in the program summary; the form of parameters of both these procedures will be shown in this part in sections describing separate procedures. 3.1. Symbolic mode procedures A new statement REMWEIGHT is performed. Using this statement with any argument, for exampie REMWEIGHT T; causes the removal of all declarations of the weight, which have appeared in statements WEIGHT in the program earlier. Thus the system no longer executes the neglecting of terms given by these statements. The clearing of the weight is done by the clearing of the global variable of the REDUCE system WTL! *, which contains declarations of the weight. Symbolic procedure SIMPSUBE defines a new operator SUBE, which is similar to the REDUCE operator SUB,but differs from it by evaluating the
3.2. Necessary declarations of variables, operators and arrays The top-level variables ORD, S and MODE are declared. These variables have to be top-level, because they are used simultaneously in more than one procedure to hold some value (similar to the FORTRAN COMMON variables). All operators used in the program are declared and also all auxiliary procedures have to be declared as operators, because the main procedure where they appeared preceded their definitions in the program. The LISP language functions MIN, MAX, FIXP and EJECT are declared as symbolic operators, in order to allow their use as operators in the algebraic mode. The operators MIN and MAX give the minimum and maximum of their arbitrary number of numeric arguments, respectively. The operator FIXP is a predicate, which is true only if its only argument is an integer. EJECTO; means a new page in the output file if this actual output file is the principal output file LISPOUT. The auxiliary arrays COF, MI, ALGSOL and FACT are declared. Array FACT is filled in by factorials. 3.3. Main procedure The main procedure MODSTAB has eight parameters. DFE is an expression representing the differential equation (2.1), DFE = 0 with two independent variables, i.e. indeterminates X and T, and one dependent variable, i.e. the unknown function U(X, T). Derivatives are represented by the REDUCE operator DF. Restrictions as to which differential equation can be used, have been given in the program summary. FDE is an expression representing the difference scheme (2.2), FDE = 0, which approximates the given differential equation in the grid of mesh points in X—T space with increments DX and DT, respectively. The values of the approximate solution at the mesh points are designated by FU(i, j), and the mesh point about which we want to perform Taylor expansions during the calculation of the modified equation is i =j = 0. DFT is the only either time
R. Liska
180
/ Accuracy analysis offinite difference methods
or mixed derivative appearing in the given differential equation. Parameters N, M and K are integers, which specify the neglect of terms DX’ for i ~ N, DT’ for i ~ M and DX’DT’ for i +j> K during the calculation of the modified equation. PHEUR is the flag, determining the possibility of using the heuristic method of stability analysis. If it differs from the indeterminate HEURISTIC, the procedure executes a stability analysis by the modified equation and von Neumann methods if and only if the differential equation is linear with constant coefficients, DFT = DF(U(X, T), T) and the difference scheme is a two-level one. If the flag PHEUR is the indeterminate HEURISTIC, the procedure executes a stability analysis even if these restrictions do not agree; then in this case it is performed by the heuristic modified equation stability analysis method. MSG is a string containing some information about the given differential equation and difference scheme, which will be printed in the head of the output. Expressions in the arguments have to be written in the usual REDUCE (FORTRAN like) syntax. The notation of the mathematical symbols used in the text of section 2 in the program is as follows:
derivatives of a lower order than the order of DFT, which would mean that the difference scheme used a bad finite-difference approximation to some derivative. The modified equation is written without the terms neglected due to N, M and K, and without the time derivatives and space derivatives of a higher order than ORD, which have been included in the modified equation during the elimination of time derivatives. The procedure examines the truncation error term TER to verify consistency. The necessary and sufficient condition for consistency is that TER does not contain any zero or negative powers of DX or DT. From TER, the procedure determines the orders of accuracy of the given difference scheme in the increments DX and DT. Depending on the given differential equation, the difference scheme and parameter PHEUR as mentioned above, the procedure analyses the stability of the difference scheme. For the von Neumann method a special procedure NEUMANN was devised, while the modified equation method is a part of the procedure MODSTAB. The modified equation stability method was described in section 2.3. Using the operators MAXPOINT and MINPOINT the procedure determines MR, ML,
X
NR and NL. The procedure MINEVENORDERDF determines R. Then the conditions of stability (2.9) and (2.12) are written to the output file. In the heuristic approach only the necessary
-
x
ML
-
m1
—
mr
U(X, T) IX
zIx NL
—
n1
—
z~t NR
—
~r
—
r
—
S.
T
—
t
DX
—
DT R
MR
S
—
u(x, t) ix
M
—
m
TER
—
TER~(u,~x,
-
~t)
Derivatives in the program are represented by the REDUCE operator DF (DF(U(X, T), T, i, X,
j)
=
~
t)/at’ax~).
After initial declarations and printing of the head of output including the given arguments, the procedure verifies if the restrictions given on the differential equation are satisfied (see program summary). In the first step of the modified equation calculation the procedure expands terms FU(i, j) of the difference scheme in a truncated Taylor series up to ORD = MIN(K, MAX(N, M)) order derivative terms. The elimination of time derivatives in the second step proceeds in the sequence (2.5). Then the procedure verifies that the calculated modified equation does not contain mixed
condition (2.9) is written. We have to mention that in the von Neumann stability analysis as well as here, the conditions of stability are not always written in the simplest form. They can be modified to a simpler form by considering that DX > 0, DT> 0, etc. 3.4. Auxiliary procedures The procedure NEUMANN executes the von Neumann method of stability analysis of a difference scheme, described in section 2.4. It is used from the procedure MODSTAB or it can be used alone for any two-level difference scheme, which approximates a linear partial differential equation with constant coefficients and with two independent and one dependent variables. The procedure has two parameters XX and YY, which give a
R. Liska / Accuracy analysis offinite difference methods
difference scheme XX =0 and a differential equation YY = 0 in the same form as parameters FDE and DFE of the procedure MODSTAB. Variables ORD and S indicate whether procedure NEUMANN is used alone (then ORD = 0 and 5 = —1000), or from procedure MODSTAB. First the difference scheme XX is broken into two parts according to the time step as in (2.8), then using procedure EXPSIN polynomial S(z) = MODE is determined; conditions (2.25) have been applied to this polynomial. S( z) is the polynomial in z of degree s. If s ~ 3, then the procedure finds other conditions from (2.25) by examining the polynomial S(z). For s=0 the necessary and sufficient condition of stability (referred to further only as the s-condition) is S(0)> 0. For s = 1 the s-conditions are S(0)> 0, S(1) ~ 0. For s = 2 and a solution of the linear algebraic equation DF(S(z), z) = 0, the s-conditions are S(0)> 0, S(1) s 0, 0
z2 solutions of the quadratic algebraic equation DF( S(z), z) = 0 the s-conditions are S(0)> 0, S(1)~0, 0
~CqFU(q,
1),
(3.1)
q
it executes substitutions FU(q, 1) = exp(iqJ), (2.16) and (2.18) and returns
(3.2)
2
~C q
exp(iqO)
=~a ~ J
(3.3)
2(O/2) (see (2.16)— as a polynomial in zis =used sin in the von Neumann (2.19)). The procedure stability analysis.
181
The procedure ALGEQUATION has one argument, an integer either 1 or 2, and it solves either a linear or quadratic algebraic equation, which has its coefficients in the array ALGSOL. The solution of the equation is returned to the array ALGSOL. The procedure CONDWRITE has one argument and is used for repeated special printing. The procedure TAYLOR has two arguments, integers A and B, and calculates truncated Taylor series of a continuously differentiable function U(X, T) at the point (A * DX, B * DT) about the point (0, 0) up the ORD order derivative terms. The procedure MINORDER has two parameters X and DD and determines the minimum degree (possibly zero) of indeterminate DD occurring in the expression X. The procedure MINEVENORDERDF looks for the minimum r> 0 so that DF(U(X, T), X, 2r) appears in its only argument. The procedure returns r as its value and fills array MI by the coefficient of this minimum even order derivative of U(X, T) with respect to X and the next even orderone. 3.5. Operators The operator AD determines the order of a derivative, which is its first argument, with respect to an indeterminate, which is its second argument. The operator DFP has one argument and its value is 1 if this argument is a non-evaluated derivative and 0 in the other cases. The operator FUP has one argument and its value is 1, if this argument is the operator FU with two arguments, and 0 in the other cases. The operator LINEARP has two arguments, the first is any expression and the second should be T. Before using this operator we must declare that X depends on T by the statement DEPEND X, T;. Then the value of the operator is zero if its first argument represents a linear partial differential equation with constant coefficients with two independent variables X and T, and non-zero in the other cases. theThe formoperator SUBONE has one argument in ~
C 1~FU(i,j)
(3.4)
R. Liska
182
/
Accuracy analysis offinite difference methods
and adapts it to the expression ~ FU’
~
guage, i.e. the creation of a new file MODE, which (3
~\
I.]
whose form is needed for the argument of the following operators. The operator TWOSTEP has one argument in the form of (3.5) and if the value of j in this argument is only 0 or 1 then the value of the operator is 1, in the other cases 0. The operator verifies whether the difference scheme is a two-level one, The operators MINPOINT and MAXPOINT have two arguments, the first in the form of (3.5) and the second is an integer. They determine, which mesh points are used in (3.5), i.e. integers NR, NL, MR and ML. 3.6. Input and output of the program As input the program requires only a call of the procedure MODSTAB or NEUMANN with the corresponding arguments described above. Because REDUCE is an interpretative language, this call can be made only after interpreting all program MODST definitions (lines 1—738 of the program). The output of the program would be directed to the file LISPOUT and will contain information from the analysis of the given difference scheme.
4. Test runs Test run no. 1 must be the first of the three. It shows the CHKPOINT facility of the LISP lan-
contains the necessary part file of the REDUCE tem and our program. This occupies aboutsys50 kbytes on the disc drive. The other two test runs show the use of the file MODE. This access is suitable for the repeated use of the program, because interpreting all the REDUCE program MODST requires 69 s with storage 450 kbytes, while reading the file MODE into the storage by the use of the LISP function RESTORE requires 1.8 s only. For single use lines 740—742 of the program can be omitted and line 743 replaced by calling procedure MODSTAB or NEUMANN, without using the file INFLE. We have to mention that the procedure NEUMANN does not need the whole program and for using only this procedure the program can be modified to a smaller one (all requirements of the procedure NEUMANN are written in the comment to this procedure). Test run no. 2 shows the heuristic stability analysis and no. 3 shows the separate use of the procedure NEUMANN.
References [1] R.F. Warming and B.J. Hyett, J. Comput. Phys. 14 (1974) 159. [2] AC. Heam, REDUCE 2 User’s Manual (University of Utah, Salt Lake City, 1973). [3] C.W. Hirt, J. Comput. Phys. 2 (1968) 339 [4] L.D. Cloutman and L.W. Fullerton, Report LA-6885-MS, University of California, Los Alamos (1977).
R.
Liska
/
Accuracy analysis offinite difference methods
183
* (U .5-
a C) 3
+ • • P
5 3 4
4 P 4
I-. 5
3 * *
x
* * • • • • 3 3 • 3 $ 3 S 3 • 5 •
* * 3 S 3 4 * 3 3 4 • .5 * 4
* 4 3 * * 4 5 3 3 4
C.,
$ 3 .5 4 P * 5 3 * 3 3 * 3 3 5 * 3 3 4 4 S *
P..
In ~ 0 a P.’ ‘4
—
4 4 3 * 3 3 3 5 * 5 3 5 3 4 3 * 3
3 * 3 5 .5 4 5 * P 3 5 * 5 4 3 3 S P 5 3 3
a
•
•
* .5 * * 5 * .5 3 P
0 UI _l ‘4 IZ
•
UP
• 3 5 3 • * 3 • S * S * 3 4
& UI U. I~
5-’
p
In
•
0
S
UI
PC
•
~,.
II.
.4
NP
0
MI
U. U) 0 ‘4 ‘4
PC
$
0
•
1
9P
MI
us U. U..
U. U.
U. 0
5 S 5
p.-
lIP
~
)“ —1
UP
‘p
~.
I— ~ 0
I—
‘4 PC ‘4 U. 0 U.
* 3 * *
0 MI PC CD .~.
‘p
up
•
S * * 4
$ 4 4 3 * * 4 5 3 4 * 4 4 * 5 3 5 * 5 p 5 S 4 5 S 4 3 5 4 * S 4
w
0
UI ~..
w PC
I..)
U) MI IJ ~C MI ‘4 U) U. U.
S
5 5 .5
3
* ‘P * 3 *
$ 4 .5 S S S
j
MI C.) PC.
0
5
~
PC MI PC (a U)
UI U.
• •
U. ~ — In U. 0
• •
~
..C
3 3 S • 5 3 5 5 3 .5 S 3 S S S 5 5 5 *
—
a ‘4 a I— U. ‘4 U.
4
5
‘P .5 4 5 S 3 5 5 5 * * * $ S $ 3 4 4 3 3 5 5 * * 5 5 5 4
.5.
3 —
C.) 4
a . ~—
I
5—. U.
I -~ Cv
5-.
(U ‘C 0
$ 3 4 3 5 3 5 5 $ 4 3 3 3 5 4 5 * 5 3 * 5 P 4 3 S 3 3 .5 5 5 ‘P 3
I (U 5—
a 3 (U 5 0
-
a
U.
+ I—. 0 5i’s
U. U.
0 ‘4 0 PC PU ~ I
PC
-
“~
0
4 U.
4
‘4 —á I — • 0 PC PC ~ U.
p.~ a~ MI p.”
* 4
& 0
U
* S
U.
(U
5
~ 5
—
•N
5 .5 $ S
PC 0 — i-
U 5 N
$ 3 5 * .5 * S $ * .5 4 * S 5 * 5 * 4 $ 5 *
~ 0 UI
S
VP 5 A
a
3
0
PC a 3 I’0
U. U.
0
N
_I ‘4 .‘. IPC UI U.
UI
-
p.-
s
N N N N N N N N N N N N N N N
‘4
I—’
p
~
+
~a U)
PC
U) Ca PC
.S
U) U.
P.’ PC
~
‘ U. 0
UI U. U.
a a
nN 5 N N N N N N N I N N N I
U
•
0
‘P S 3 P.-
PC
0
p
S IJ
UI ..J W
-P 5 3
PC
UI
PC
a
PC
0
* 0 .. 5’
5 CI
P 0 U.
U.
(U
P.a *
4 5 * 5 3 * 4 * S * 3 3 * 3 3 4 * * 5 5 4 3 5 3 5 3 4 4 5 5 5 .5 3 4 3 S 5 3 5 * * 4 S 5 3
(U (a 5-
3 C) * —, PiP PC — P.. ...
PC ~
0 — ‘4
II.
0 4
0 PM
.
UI —
-~ (U
PC 0 ~
(U
p.. a
.J ‘4 ‘P
I’PC UP ‘4 UI U. U.
a
4 Cv 1.1
PC I-. ~-.
* I... 0 .5 (U
5
S
PC a I
0 0 PC
P~ * (a
N
(U
* * * S 5 * 5
a
‘P
—
.
U) ‘4
0
I ‘0
*
4
I
a
CI
4 4 PC
0
~
I.
I-
5 4 P” 0 ‘P
S..
U..
0
“~
I
UI PC UI
-
PC -~
5 5 ‘P 5 3 5 ‘P 5 3 5 S 3 $ 5 * S 5 5 3 S 3
PC.
0 3 (‘U
3 5 P
—
(a
4 ‘P •
PC
PC 0 —. I” ‘4 ~ 0 UI 0 UI IC
I 5 N U N N N N N N N II II II N N U N N
P..
~
p.-
U.
a
pp
UI U.
5 $
U.
(a ‘4 U. ~
5 * 5 $ * 5 * S * * S 5
PM PC UI PC ~ U)
CI CJ
* ‘P
‘4
UI (a PC UI U.
Cl) U.
3 ‘P * * * 3 S S ‘P $ *
PC UI p... lIP a U) PC 0 Ca UP
PC ~
S P S * 4 3 3 4 3 3 * 3 * * 3 4 4 3 4 4 * 3 4 * * 4 3 5 * 3 3 $ 3 S 3 3 3 .5 S 3 54 0 4 P * * (U 3 PC * 0 5 3 UI 4 ‘4 4 ‘4 5 5 UI 3 PC S MI S PC S p,~ 5 4 * ~I S (a 5 PC 3 UI * U. S
a
0 U. 0 ~.-
U. 0
U I-’ is
P.” -
x
U) U.
u. a
0 U. 0
PU
0 ‘4 0
S
UI PC
p.-
* 4 S * *
R. Liska
184
..J
‘4 PC ‘4
•
P’ I” a
N N N ft N N U
-a a In
‘4 I— U)
I-’ ‘4 fl
I
UI
I P
‘.
PC 5-
PC 0 .
N ft
a
N U N N II N N N
0 PU
pp
SI U.
0 0 PC
p
I I I I I I I I
* ‘P * * 5
S .5 * 5 5
U) PC 0
* * $ * *
$ * 5 S
* ‘P
* **
* *
— I-’ ‘4
S $ *
* 5P
S 5 ** $ *
‘p
5
‘p
5 *
‘4 U.
S * 5 * S S S S * S $ * S * * * * * * * $ 5 * * $ *
.4 ‘4
* 5 5
UI PC UI PC (a U)
‘p ‘p
5 $ 5 *
MI ‘-5 ‘4
* * $
.5* .5 4
—
‘p
5-. PC U)
0 U. 0
a
P” I’-
$ * * * * *
..J — U) ‘4
I
PC
0 o
CP
S
PC
I
—
P.—
0 PC 0 Ca
0 A
I— PC
U) ~
U) a
—“s
(‘U
.SI IC U.
~ U)
+ (SI
UP
0 * (SI ~
,-
U. ‘4 U)
I
U)
S
UI
S p.0 S
4 * * * 4 * *
‘p
0 II,
* * *
U) ~p
‘p
‘p
up
* ‘P 5 S
$ 5 $
LI
PU PC (U
5 * 5 * S ‘p
(a 5-
S $
LI UI PC
0 A
U)
PC PC 0 ~
$ S $ $ $ * 5
I—’
a’‘4 ‘4 U)
up
S $ ‘p $
U)
$ * $
.,j • • •
‘p
$ * S
I’.’
*
CD
(5. —. III
I
I-’
PC
0’
~‘
C> ,..
II.
5’
— ‘4
* * 5
$
~
‘p
_., III .-4 PC PC
a PC U)
I-.
*
. (.5’
$ * *
0
0. pp — ...I
‘4 PU 5-~ PC PSI
up PC
—.
‘4 p.’
(“1
0
* $ S $ 5 * S 5 ‘p
C”
~
* *
—
P.’
* $
PIP a
‘4
*
‘4
U) SC PS. .‘. 0
‘p *
~ PU
S
PC
‘p *
$ 5
II.
‘p ‘p
U.
0
(5.
~ I—
~) a
* ‘p ‘p ‘p
S
— (5) P” —I
‘4 PC ‘4
I-’ ‘4
Ca ~.
0
4
$
~
$
— (5) UI 0
$ ‘p $
PC Ca ‘4 Ca U. 0. 0, ‘4
~
PU PC
4 4
5 * 5 ‘p ‘P $ $ S P 4 * $ * 5
PC P55
~ *
‘p
LI
‘4
0
•
PC PU
C)
PC
*
pj
(5)
(51
p
* $ $
I— a *
PC
•
0 PC 0
(SI
U) U) L( MI
p-.
•
‘P $ 4 ‘P S
LI
*
a
+
0
‘~
4 *
PC 0
PC W PC “
5
PC 0
* * S * S * *
1U
‘4
* $
‘p 5
4
II,
‘p
..a ‘4
(I)
‘p
U)
0
* ‘p
S 5
‘4 0
‘4
‘4
is
a PC ‘4
‘p
‘4 Ui
‘4 ‘4
Cv
S S
U)
U,
Ui
(F)
a
PC
*
0
U. ~
PU -)
~ PC
5
PC
* * S * 5
0
‘4
MI
* * *
0
P.
* 5 .5 *
‘p
SI
5 5
-3
5 *
5-’ ‘4 ‘4 0.
I—
‘P $ $ * * $ * $ *
* 5 $
‘4
$
— (.5 .. IL
$ $
$ $ *
0
,0
* * .5 5 * *
p.PC
S * $ $ *
‘p
(.5
‘p
$
~
0
I I I I I I I I
5 .5 $ * S
* 5
‘4 ‘4 U)
0
PC .4 PC ~
-.5
,.,
I-’
PC
$ * * 5 5 S $
* .5 * *
—
I I I I I I
‘p
“ PC UI ‘4 UI U. U.
PC ~
“.
5-’ Us PC
a U.
0
PC 0
‘p P.’
PC ‘4
5 S *
$
* 5 5 $ S * $ * $
‘p
p... up
* $ $ S S 5 S 5 * * * * ‘P * $ * * 5 S $
(a U. ~ (5)
s
U) PC 0
Ca
p
U. a
‘~
S $ p $ * * • * * * 5 * * $ 5
PU 5PC MI U. UI II,
•-.
‘4 0 U.
I— ‘4
* $ * • *
* $ **
‘p 5
up
I I ( I I I I I 5 I I I I I I I I I I
S .5 * ‘p * $
$ 5 * ‘P *
~ 0
UI ‘4 5-,
0 0 PC I— W PC
* • $ * * *
‘p
* * *
*
P. P.—
U)
* * * * S *
* S ‘p
**
U, 0
)-
5 S * * $ $
UI
O
N N N I N I N N N N N 5 N N I II
* $ S * $ * S S * $
* 5
Us C) PC PU ‘4 Mi II. II,
—‘
* S $ * * S S S $ ‘P
UI
MI PC UI PC (a NP
C’)
/ A ccuraey analisis offinite difference methods
5 5 *
5 ‘P $
C) Uj
.5 S
‘P
$ S
I
0, PC
‘p
$ 4
‘4 UP IU
5 S *
S * $ $ 5 $ ‘p
I (5) •
0 PC
* $
$ 4 $ $ * $
$ 5 $ *
PC ~
‘p
‘4
$ *
5-
$ * $ 5 ‘I’ *
4 * $
U) PU I”
4 ‘p $
‘p
$ 5
*
5
.5
5
R. Liska
/ Accuracy analysis offinite difference methods
185
I-’
5-
II.
PC
+
5U.
is
0
((I PC
I
a .5
50 5’ 515
PC O .5 I— 0 * 3) * Cs) I 5-
*
is
0 -
0 U..
.5 is
PC 0 S I-’ O 55is
I 5-
* is 5• 0 5U.
+ is (5)
PC 0 * 5-0
o
PC 0
PC 0 — ~—
“
‘4 ~
a U) .4
.4 — p.PC PU ‘4 U) U. II., — 0
N N N U Cl ft 5 U ft II N II N U ft ft ft
• p..’ •
P.0 *
is
(.5 5-
PC
UP +
Ui PC UI 14 U)
is
(U •
PC is
I— • PC 5-
~ 555.
0 * (a
‘4 0 U. “’p
is
* * I— 0 *
U) Ca
PC U) U. PU U. IS. ,_.
0
N N U N ft U I II N N N II (I II N N N
0
SI
pp
II
• — 5’
..a CD
‘P ~ 5
MI
PC
*
PC
0
5 * * *
•‘
5
U..
+ is
PC 0 $ I”. 0 5’ 55” ‘ 5U.
(51 pp C 5is
Ca $ is
0
• is ,.‘
I
U) ‘4 U) 0 ‘4 1’)
PC
* $ PC 0
‘4’ 5 5 P.C
SI
5-’ C)
•
PC 0 5 P.’ 0 * ‘5
I 5~
* (a S
* * * * * 3 S * * * * * * * *
* * * * * * * * $
is (I’) •
PC • -5
5-. •
.4 5I U.
0 +
‘0 5is
Cv
PC 0
PC 0 00
PC
5‘4 0 UI
* 5-. 0 *
_I ‘4
(a
* pp
a
I-’
5-. PC UI ‘4 UI
$
U. IS,
+
— 0
(5) 14
* ((4 I
5-. ‘4 ~ 0 Ui 0 U) — U. —
0 0 PC
PC
LI U)
UI PC U)
UI
UI U. U. —
.4’
I-’ U)
•
.-.
.4
U) PC 0
U.
Ca
P.’
• is ~.
a
MI PC Ui
‘4
is
PC 0
PU ~ ‘4
(.5
I— PC
$
ft II U N II ft II N N ft N II N N II It N ft ft
PC 0
—
K
*
•.
5-’ 0
PC P.-
5’
‘p
Ui
0
~
* * $ * * * * S S * * * * * *
5-
* * S * S $ * * .5 * * $ 5 * $ * $ $ * $ * * $ * S * S * S $ * * * $ * * * * * * * * * .5 5 * * $ 5 * * * $ * P * * $ * * 5 * $ * S * 5 5 * 5 * ‘p $ * $ * 5 * $ $ *
*
• 5’ ~ 5’ U.
(5*
I
is PC
ft PU II U
‘0 5 A -) *
5— 5is
U
Cv
* ‘P 3 * $ * * * * $ * 5 * * 3
LI
PC
$
n
is
—
~
(%5
•-
•
U)
U. 0
PC
PC Us PC
N
0 5-
5’
•
0
‘4 ‘4 ~P (a (4 ‘4
-‘ U)
U.
PC
U)
U)
•
14
‘4 U) 0 ‘4 C)
0
is
P.is
5-. •
PC 5’ U.
0
PC U) ‘4 Ui U. U. — 0
UI
PC P.
$ $ 5 S * S * * * * $ * * $ * * * * 5 * * 5 * * S 5 * 5 $ 5 * * * 5 * * $ * * $ * S $ S S * * * * * $ ‘P 5 .5 * * * S 5 S $ $ $ $ 5 $ $ $ * * $ $ * 5 $ ‘p
$
~ — UI PC III PC CI
5) U) C)
PC UI ‘4 (II U. U. PSI
0 U. 0 P 5a _I
— In ‘4 5-’ U) ‘4 0 IC
PC 0 5a
0 PC 0 Ca )_
U. .4 U) U) U) (a
U) PC 14 I’S
I-’ Cl)
0
— ‘4
A
~ MI PC
(.4
$ $ $ * * 5 * 3 .5 * $ * * 5 * 5 * * $ $ * * * $ * * * * $ * * S * * $ * * 5 $ * * * $ $ $ * * 5 * * ‘p $ * * * $ 5 $ $ * $ ‘P ‘p * * ‘p $ * * * * * * * * $ 5 * * $ * * * S
I
* *
up
U)
.4 ‘4 P
p •
5). (I) ISI
.4 CD PC a
‘4
U) PC -3 0 ‘0
‘4 P.‘4 0
UP
‘P
55.
0
UI
OP.” PC PC
.~ PSI I’s
UI PU
P. PC I—
C) PC U) *
S *
R, Liska
186
/ Aciuracv ana4’sis offinite difference
methods
* 0
C’ U.
* I0 * PSI
5 $ $
p
* ‘p * $ * * * * ‘p * S $ * * $ S * $ .5 $ * $ * * $ $ * 5 * .5 * * * ‘p * * S 5 * * $ 5 5 * ‘p 5
•
4 IS ¼* U) PC PC is 0 0 is * 5- (U (‘.1 is 5-is PC 5-C) • is•is is ~
5-. 0 I 5’• 5 U. C’ PC 5- C’ S 5’ C’ is U., ICC> C I • $ isC) CD 0~ — • C’ U) s”U. 1 5- * -5 (‘4 PC U. I •
• •
I”.i ‘4 t 5”
C’ I—
5- ~“
‘CD
5 C’ I— (4
5-’ U)
+ s’. PC I’-.-- 5’ C) C’ C’
C’ P.
5- U. 5is 5’ U. is 5 C) (C> is * • (P (.5 0. — + I~)is
‘4 ‘4 C’)
C)
CD + 5U. (U •
‘4
0,
“
is s”
5(‘4 U)
.4
(4 C) 55) ‘4
• 5- CD CD’-’ PC CD ‘4 IL CD C’
P— I-’ C) • S PC (5) 5* C” * 5’ (-(PU 5- 0.
‘p 5 * $ S 5 $ * $ * $ * ** * 5 $ $ * $ * 5 $ $ ‘p * ‘p $ * $ 5 $ ‘p
* * * $ * $ ‘p S $ ‘p * * ‘p 5 5 $ $
US
4
PC
$
‘P $ *
is
0 • 5— 5-
0
C’ U. * PC * P._ C *
N
(U
PC C $ I-’ C S (V
C)
S is
0
5. is is —‘ •
s’ I 5’ 5-
0
C’
‘.‘
C’ U. S
U.
$ PC 0 * 5-
is
(“.1
PC 0 $ (U
C
* (4
+ I . is
is
0
C)
C>
•
•
C)
0
5’
5’
C’
C’
IS.
N
(U (U
,_,
PC I’(‘p $
•.
W U) (.I ~‘
W ~
PC U) ~ Ui U. IL 5’ C> IS.
0 (~ — ~j)
P .4
‘4 C’ ~ P (_ — .4 — 0
‘4 I— (~
ft ft U (( It
•
$ (Si
I—
I * + (‘.1 is
U’ CD — .5 is.
C) IT fl (I It ii I I II Ii (I It II II I fl (( (( II II II II (I II II ft II (I II II H (~
• C”
C’ U.
$ (‘J
“4
P.0 *
5-.
C
rJ
• CD 5-
•. PC I’C) — .5
+
.4
CD —•
‘4
IS.
I’..
p.. PC Ui ‘4
Ui
(4
4
(5. (5.
I
CD
0
is
— VP
CC
*
is ,—
is 5” I’ .‘. “-
I
‘4
C.
0.
(U
$
5), ‘4 = (4
I
(‘~ (4
I 5-
CD — U)
I’ PC 5-
I
$ I— $
I—
—
U.
p.-
• .4 • is
CD
CD
C)
(5)
‘5. 5”
(1. (5)
is
PU p.‘4 P.
CD $
U) S
((P
PC C)
‘-‘ 5-
• .4 5U. C’ 0. 5-
CD C) PU
CD * 50. $ -4’
$ * $ * * S * * * * * * * .5 $ S ‘p * * $ $ $ 5 * $ * S S $ * * $ $ 5 * $ * $ * * 5 * 5 * * *
V ft (U V
o
U.
C) ‘I’
U) ‘4 ‘4
is PC • is
P” • PC .-
C’ PC
5’ (5. 0.
* C-)
N A
MI PC Us PC
is
.4’
C) (F)
PC
0 *
Us
5-
(4 PC
0 5’ 5is CD S~ ((P
U)
‘4 U) IL U. 5’
*
0
(~
0
0
5
S
IS) a’-
S 5 * * $ * $ * 5 ** $ * $ * *
P.-
+
— .4 I-S 5) “4
is (54
CD a
I—
U) * -3
(5) PS. ‘4
0
(5.1
C’)
PC C)
*
$ S * * * * $ $ * * $ $ * $ $ 5 5 $ * 5 $ $ * $ * * * $ * 5 S 5 * $
PC
U.
I PC 0 *
(‘4
—
5-.
(4
* * *
* $ *
$ $ * S * * $
** * * S ‘p * $ * 5 * $
‘P $ * * * $ * * * ‘p * * $ * * $ * S * * 5 * * * $ * $ * $ $ * $ $ $ * * *
‘4 0
$
*
U,
5-
+
PC
C’
(P
(a
—
(5) * P.’ C ‘p
I—
CC PC I— LU
P. PC CD
‘4 P. C’ PU CD PC
CD ~
I I I I I I I I I I I I I I I I I I
C’ Sn
(54 (
$ ‘4’
C’
I
A
$ * * * ‘p * * * 5 ‘p 4 S * * * $
C PC ‘4 (5’
‘4 ‘4
(5) VP PU C.’ Ui PC
is
.4’
5 I—
* 5 $ ‘p * $ $ $ $ $ * 5
(5
$
S
*
*
(‘.4 5’
(U 5-
S
is
5)
P.5>
o * I-
$ .‘S C-’
I is
CD — (5)
$
.4 0 $ I” C 5’ 5is
CD —
U) * (U
PC 0 $ (U +
•
5PC UI —
$
— U. Ii.,
CD
* * * * $ ‘p 5 ‘p * 5 $ 5 5 $$ $ 5 * S 5 $ * 5 * $ * ‘p ‘p * $ * * $ S $ S 5 * * * $ $ $ *
LU
$
CD —
* 50 $ -3 1 (‘4
(5)
PC CD
(‘4
0 PC ‘4 a’-
U) ‘4 ‘4 U) PU (.1
PC
U)
$ 50. *
a’.’ 0
(.4
—
C’
“4
U) ‘4 .4
I—
+
‘4
P-I C
CD
U)
PC C) (-5
(U
PC 0 -I —~
0. * (U (.1
(.5.
= C) ‘4 ‘5’
PC 0 — — —
* ‘.5
0. CD 0
I
C-’
(‘C
PC
‘4 ‘4 I,’)
I—
C PC (4 — PC
(5) 5
5C) $ -5 (4
I
‘p
14
C’ U)
(U
SI
C)
UP
N A
* 5 * * 5 * $ * 5 $ $ 5 * 5 * $ * * $ $
—
UI — I.) — U(S.
0
* * $ * .5 5 $ * * * * $ * S * * ‘p 5 ‘p $ * * * * * S $ $ $ S *
I— I-’ — ‘4
PC
* * ‘p $ $ 5 ‘p
5-
*
I’
U) —
U) C’ —~ •
‘4 >
0. U’ — ‘.~1
pP
‘4 P.-
~
C’ CD
0 ‘.5 CD ‘0
5’
,—
“4
0.
(5.
‘4 ‘. U) 0 I-’ PC PC U) U)
(C> .4
Q
PU — CC PC
U)
P.’
* $ $