Journal of Scientific Computing, Vol. 5, No. 1, 1990
Parallel Multischeme Computation with the Switch of Computation History Yi-ling F. Chiang 1'3 and Jih-Fu Lai 2 Received February 28, 1990
A parallel multischeme computation in the solutions of differential equation initial-value problems has been studied. The mathematical switch of computation history is used successfully in the identification of the best approximation among all available ones at a computing step. A solution correction factor is also developed to achieve an extra four digits in solution accuracy. Based on our results, if the truncation error of computation history as we defined it can be properly utilized, then a computation engaging high-order schemes or using fine grids may be unnecessary. KEY WORDS: Parallel computing; initial value problems; differential equations. 1. I N T R O D U C T I O N A parallel m u l t i s c h e m e c o m p u t a t i o n for the solutions of differential equation initial-value p r o b l e m s has been studied. A m a t h e m a t i c a l switch (Chiang, 1986) designed from the c o m p u t a t i o n history of a scheme has been e x a m i n e d for its use in the identification of the best a p p r o x i m a t i o n a m o n g several available ones at a c o m p u t i n g step. The c o m p u t a t i o n h i s t o r y of a scheme refers to the a p p r o x i m a t i o n s of various mesh sizes c o m p u t e d from the same scheme. In a m u l t i s c h e m e c o m p u t a t i o n , at every c o m p u t i n g step, m o r e t h a n one a p p r o x i m a t i o n results. Intuitively, one of the a p p r o x i m a t i o n s m u s t b e a r a smaller a b s o l u t e e r r o r to the exact s o l u t i o n t h a n the others. If the a p p r o x i m a t i o n of the smallest a b s o l u t e e r r o r can be identified a m o n g all c o m p u t e d ones, then the c o m b i n e d solution of several schemes achieves an a c c u r a c y at least as one of an i n d i v i d u a l scheme The William Paterson College of New Jersey, Wayne, New Jersey 07470. 2 New Jersey Institute of Technology. 3 Present address: West Virginia University, Morgantown, West Virginia 26506. 35 854/5/1-3"
0885-7474/90/0300-0035506.00/0 9 1990 Plenum Publishing Corporation
36
Chiang and Lai
engaged in the computation. The use of the switch of computation history is successful. The best approximation can be identified correctly at almost every computing step in all cases that we ran. A solution correction factor is also developed from the computation history of a scheme to achieve an extra four digits in solution accuracy. The corrected solution of a lower-order scheme approaches the corrected higher-order ones, hence, a high-order scheme computation seems unnecessary. In a partial differential equation initial-value problem, the corrected solution on a coarse grid approaches the corrected ones on fine grids; therefore, a fine grid computation also seems unnecessary. Several mathematical switches have been designed and examined. The use of a switch of exact error confirms that a combined solution of several schemes is indeed more accurate than one of an individual scheme engaged in the computation (Chiang et al., 1988). However, the use of the switches of estimated errors were unsatisfactory. The truncation error estimate of a finite difference scheme is established by substituting Taylor series expansion into the scheme and by eliminating higher-order terms in the resulting series (Hamming, 1962). The existing finite difference schemes were derived mainly from the given equations in a problem, which consists of both equations and conditions. Hence, the establishment of the truncation error estimates is based on only partial information about a problem. A switch of estimated truncation error lacks the sensitivity necessary to identify correctly the best approximation among all available ones. Its use is unsatisfactory. A summation of the estimated truncation error at every computing step is the estimated accumulated error. However, this summation may not overcome the weakness of its elements; therefore, the switch of estimated accumulated error fails also. The computation error is a combination of the roundoff and the truncation errors. The machine-dependent roundoff error occurs at random, but the problem-dependent truncation error (Chiang, 1988; Richter and Falk, 1984) characterizes the scheme and the problem. In a computation, if the truncation error dominates, then the computed approximation should carry both properties of the scheme and the problem. The approximations of various mesh sizes should share the same property of the scheme in the given problem. A switch derived from the approximations of various mesh sizes should have the necessary information to identify correctly which scheme produces the best approximation (Petzoid, 1984). A two-term truncation error estimate of a finite difference scheme is suggested to replace the conventional one-term estimate. Based on the one-term estimate, a two-point truncation error of computation history of a scheme is derived in terms of the scheme approximations of two mesh
Parallel Multischeme
Computation
37
sizes. However, based on our two-term estimate, a three-point truncation error of computation history exists in terms of approximations of three mesh sizes. The switch of computation history compares the truncation errors of computation history of various schemes to determine which one bears the smallest magnitude. In the solutions of ordinary differential equation problems, the truncation error of computation history varies with the order of the scheme, but not the scheme itself. In partial differential equation problems, even linear ones, the truncation error of computation history has to be derived individually for each scheme. The truncation errors of computation history of equally ordered schemes may not be the same, but the ones of differently ordered schemes can be similar. In a multischeme computation environment, solutions of several linear initial-value problems have been studied. The implementation in NYU Ultra Computer has been arranged to carry out the computation simultaneously both in parallel and in array to achieve computation efficiency. In Section 2, a two-term truncation error estimate of a finite difference scheme is suggested to replace the commonly used one-term estimate. In Section 3, the derivation of the truncation errors of computation history is described. The implementation in Ultra Computer is shown in Section 4, and the numerical results are stated in Section 5. In the last section, our new discoveries and our future research plan are summarized.
2. TRUNCATION ERROR OF FINITE DIFFERENCE SCHEME Let {xn} be a set of equally spaced mesh points, Ax apart, and .d. a finite difference operator of order v derived from the differential operator of the given equations in a problem. Then the truncation error in the approximation of scheme .~. is defined at the mesh points only. Let y(x) be the problem solution. Then at x = x n , the estimated truncation error, E~.(x,) is given as _ E~'(x")-e~
ddxv +~+~ l y x = ~., ~ ~.o< ~ x~ o x ~ l _<
(2.1)
In (2.1), the superscript "o" indicates an ordinary differential equation problem, and the function co= c~O(x, y(t)(x), l = 0, 1,..., v) is independent of the mesh size Ax, where y(t)(x) is the/th derivative of the problem solution y(x), respectively and y(~ y(x). Equation (2.1) is derived by substituting a Taylor series expansion into the expression ~ . y , and by eliminating the y(t)(x~) terms in the resulting series for all l > v + 1, where y, = y(xn). In general, this y(l)(x) term can be written as e~ where ~o~=
38
Chiang and Lai
S~
y(m)(x) for m = 0, 1,..., l - 1) is independent of the mesh size Ax. Since the operator .d. is derived from the given equations in a problem alone, the establishment of (2.1) is based on only partial information about the problem. Furthermore, the computer is a finite machine, and the elimination assumption of the y(t)(x,) term is too strong to be satisfied even by smooth functions. Consequently, the truncation error estimate underestimates the true error. However, the estimated truncation error, E;(x,,), of (2.1) is still shown to be problem dependent. Let f(x) be a continuous or a piecewise continuous function defined in the finite interval [ - N Ax, N Ax]. Then on the set of 2N equally spaced points, {x~, -N<~n ~
f(x,)=
2
(
knTz
\c~ksin-~--+/~kc~
knrc~
k~0
for
n = -N, -N+
1..... 0, 1,..., N - 1
(2.2)
where ek and/~k are the Fourier coefficients as defined in many textbooks (Hamming, 1962). From (2.2), the derivatives of f(x) in magnitude can be written as N~I Q~)(2s) ( knT~ kmt'~[ If'2*)(x,,)[ -~ o~k sin ~ + ]~k cos --N-) I k=0
and
(2.3)
I f (z'+ l)(x~) I =
~k cos ~
+ / ~ sm - ~ - )
for
s = 0, I,...
k~0
In the establishment of the truncation error, E~(x,,) in (2.1), the assumption that the higher-order terms, c~~ for all l > v + 1 are considerably smaller in magnitude than e~ at x = x n cannot alway be held true. Consider the following simple function:
y(x) = f(x) = ek sin krcx
for any nonzero k
(2.4)
From (2.3), we have
ly(2S~(x,,)l - [~k(krc) (a~ sin kTtxl and
(2.5) [y(Z'+l)(xn) I = Icck(krc)(2s+l)coskrcxl
for
s = 0 , 1,...
Parallel Moitischeme Computation
39
Consider an integer v. For a small positive number e, a positive c5< 1 exists such that at a mesh point x = x , , we have
but
(2.6)
ly(V+2)(x,,)[> A ( I
-fi)
where A = c~k(kn)~+z. Condition (2.6) can be satisfied if x~ is near the roots of sin k~x for odd v, and of cos knx for even v. The computer is a finite machine; therefore, a small positive number e' exists such that the mesh size Ax satisfies lAx[/> e ' > 0. This gives that
[ o~OAx~+ld~+lY dxV+1 x=r~ ] <~B~ AxV+le
(2.7)
but
a~ where B1 = Ic~~ to satisfy
dV+2y2 :,=~0 [ >~B2AxV+2A(1-fi) o and B2 = [7~162
and
Ax
(2.8)
is finite. If e is chosen
< B2 Ax A(1 - ~)/Bt
(2.9)
then the y(~+2)(x,) term is larger than the y(~+~)(x,) term in magnitude, and the estimated truncation error E~-(x,) of (2.1) underestimates the true error. The order of a finite difference scheme misleads. Let the problem solution in the interval of consideration be either continuous or piecewise continuous (Richtmyer and Morton, 1969). Then from (2.3) a n d from the special properties of the sine and the cosine functions, the following truncation error estimate can be more meaningful than the one in (2.1): E~(x,)=
dV+ly d~+2y) ~~Ax"+' d-~r+ fl~ AxV+Z~xT-27j X~
~o'
x,,_~<~r176 (2.10)
In (2.10), both functions, c~~ and/3o are independent of the mesh size, Ax. The established truncation error, E]-(x,) of (2.10) is problem dependent. Since the derivation of (2.10) is independent of the given conditions in a problem, the estimated truncation error E~(xn) of (2.10) is applicable to an ordinary differential equation boundary-value problem as well.
40
Chiang and Lai
Correspondingly, in the solution of a partial differential equation problem of one space dimension, the truncation error of a finite difference scheme of order v should be written as
EP(x"' tj+ l ) = ( O;pAXv+ l ~v+ lu ~x~+~
~v+ 2 f ~ x
+ ~ ~ Ax~+~ ~-f~-r-~]
=~ ,
x._l <~~P <~x.
(2.11)
Xn_l~P~Xn
(2.12)
instead of EP(xn,
tj+l)=
O~p AjXv+l __Ov+lu ~xV+ 1 x=~-Pn
where U = U(x, t) is the problem solution and the superscript "p" indicates a PDE problem. In both (2.11) and (2.12) the functions, eP and ~P are not independent of the mesh size Ax; however, the truncation error, EP(x,) is still shown to be problem-dependent even though its derivations are independent of the given conditions in the problem. 3. TRUNCATION ERROR OF COMPUTATION HISTORY The computation errors in the solutions of differential equation problems are problem dependent, but the estimated truncation errors of (2. t0) and (2.11) are derived only from the given equations in the problem. Their establishment consumes only partial information about the problem; therefore, a direct engagement of (2.10) and (2.11) into the design of sensitive switches may not achieve its goal. Other approaches are necessary. Consider the problem approximation of a scheme. It must share the properties of the problem and the scheme. A combination of approximations of various mesh sizes must carry more closely the property of the scheme in the given problem than a single approximation. A switch designed from this combination must have the necessary information to identify correctly which scheme produces the best approximation. Based on this assumption, the truncation error of computation history is established, and the switch of computation history is developed. 3.1. Ordinary Differential Equation Problem
In the solution of an ordinary differential equation initial-value problem, solution accuracy is the preliminary concern. Let y~ = y(x~) be the exact solution of the given problem at the mesh point, x = xn. Let the scheme employed be of order v and its approximations, ya, y~, and y3n be
Paraild Multischeme Computation
41
computed at x = x, with mesh sizes zlx, 2Ax, and 3Ax, respectively. Then from (2.1) by applying Richardson's technique (Morris, 1985), an estimated truncation error is 1
2
In (3.1), E~ is called the two-point truncation error of computation history. Correspondingly, the application of Richardson's technique to (2.10) leads to the three-point truncation error of computation history; E~3(x, ) -- 2-3 v+ ~(yl - Y nz) - 2 v + l (Yn--Y3n) 1 2.3v+1_2~+1.3v+l_2v+l
(3.2)
As shown in (3.1) and (3.2), the truncation errors of computation history, and E~ depend on only the computed approximations of various mesh sizes and the order of the employed scheme, but not the scheme itself.
E~
3.2. Partial Differential Equation Problem In the solution of a partial differential equation initial-value problem, discontinuities may occur even when the initial condition is continuous (Chiang, 1984). Hence, a mathematical switch should have the capability of choosing a combined solution, which is accurate, and in which the solution discontinuities do not smear (Ames, 1976). A single technique cannot accomplish this. In this article, in the design of sensitive switches, only solution accuracy will be considered. The discussion on discontinuity smearing is left for a later paper. We studied the advection equation of one space variable, 80
80
8t
8x
U(x, O) = Uo(x)
(3.3) is given
where the wave speed c is a constant (Vichnevetsky, 1981). As mentioned in the previous section, both functions s ~ and /?o in (2.11) and (2.12) are not independent of the mesh size Ax; therefore, the truncation errors of computation history in a general form as given in (3.1) and (3.2) are not possible. Three schemes are engaged in a multischeme computation to solve problem (3.3), namely, the Lax, the upwind, and the Lax-Wendroff, respectively (Roache, 1977). The application of
42
Chiang and Lai
Richardson's technique to (2.12) gives the following two-point truncation errors of computation history:
Ax 2 / (U2 'j'c- U2'j'L) T
~'mn'tj+l)= k
-'~X / ,--n
--On''
)
(3.4)
(3.5)
and
1 ( ca At2"~ (U~,j, LW Un2,j,Lw) r , ~ , tj+~) = ~ 1 - ~ j EP2,Lw(•,
(3.6)
where the superscripts, "L," "UW," and "LW" indicate the Lax, the upwind, and the Lax-Wendroff schemes, respectively. In (3.4), the symbol, U~'j'I- and U~'j'L represent the numerical solutions of the Lax method at the mesh point (x,, tj) computed by assuming a constant At and variable mesh sizes, Ax and 2Ax, respectively. Correspondingly, the other symbols are defined. The application of Richardson's technique to (2.11) leads to the three-point truncation errors of computation history, i.e.,
1( c2_~t2~ EP3"L(x~,tj+~)=~ 1 Ax 2 j(U~,J,L-U~ 'i,L) cAt'~ [3( U~1,j,vw 1 +-A-x-x)
Ep3'UW(x., tj+l)= -
(3.7)
U~2,j,vw )
U w - uz, j, uw)3
(3.8)
and
2 (V2j, LW_ Un2,j, Lw ) EP3'LW(x,, t j + , ) = ~1 ( 1 c~-~--//
(3.9)
In (3.8), the symbol U3'j'uw is the numerical solution of the upwind scheme computed at (xn, tj) with the time interval At and a mesh size, 3Ax, respectively. It is interesting to mention that
E~z'L(x,, tj+ t) -- E~3"L(x,,, Q+ ,) and
(3.10) EP2,LW(y
T
_ pp3,Ll y
,'~n, t j + l ) ' = ~ T
,J~n, t j + l )
Relation (3.10) is not satisfied by the truncation errors of computation history of the upwind scheme.
Parallel
Multischeme Computation
43
4. PARALLEL MULTISCHEME COMPUTATION Three approximations for each scheme are computed as shown in Fig. 1. In the figure, the approximations y,~, y], and y3 are computed from the same scheme at the mesh point, x = x, with mesh sizes, Ax, 2Ax, and 3Ax, respectively. The substitution of y~, and y] into (3.1) gives the two-point truncation error of computation history, and all three into (3.2), the three-point. The truncation errors of various schemes are then evaluated and compared to determine the one of minimum magnitude; therefore, the solution of minimum absolute truncation error is the best approximation. Only the best approximation at every computing step is saved and used in later computation. 4.1. Ordinary Differential
Equation Problem
The following ODE initial-valued problem in one space variable is studied, i.e. (Richtmyer and Morton, 1969),
~x= f( x, Y) (4.1) y(O) = Yo
is given
where both y(x) and fix, y) are column vectors of order N. Since any linear or nonlinear Nth-order ODE can be written as a set of N first-order
[ .......
I
I ~ .......
I ! .......
.!
J
y~. I
I
l
i
1
2
Yn*3
............... i
.
.
.
.
.
.........
.....
>
Yn*l
.
.
gn*3
< .....
Yn*)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .+ .. . . . . . . , . . . . . . . ~ . . . . . . .
X~_ 2
X~.~ 9 .....
.......
>
X~
Xn. 1
Yn
< .....
Xn§ 2 >
3
T: ......
. . . . . . . . . . . . . . . . . . . . .
......7.......i.......1
l Yn*2 2
, Figure 1.
~gVn~-'
Xn.]
44
Chiang and Lai
equations, problem (4.1) is rather general. Four schemes were engaged in a multischeme computation to solve problem (4.1), namely, the Euler, the modified Euler, the Runge-Kutta and the predictor-corrector methods. The schemes are all second order except the first-order Euler. The multischeme computation is implemented in the NYU Ultra Computer of eight processors. Because in Ultra the execution of the arithmetical operations is rather slow, a parallel multischeme computation as described in Chiang et al. (1988) does not pay for the overhead in processor parallelism (Ultra Computer, 1987); therefore, the computation process is arranged as shown in Fig. 2 to be carried out simultaneously both in parallel and in array to attend computation efficiency. In Fig. 2, the superscripts, "EU," "ME," "RK," and "PC" represent the Euler, the modified Euler, the Runge-Kutta, and the predictor-corrector methods, respectively. The notations will be used again in the rest of the paper. The symbols ~-2,EU and F3'EU represent the two-point and Tn ~ Tn three-point truncation errors of computation history at x = x n in the solution of the Euler's method. For comparison purposes, the exact truncation EU error, Een , of the Euler method is also computed. Correspondingly, the other symbols are defined. Furthermore, at every step, the switches, S~ and S 3, of the two-point and the three-point computation history are used to identify accordingly the best approximations, Yn(2) and y(3), to be output at that step. This gives a computation in array. The assignment s to the processors are different. In a shared memory environment, the processors are all busy at every step, except that at the first and the last two some of them are idle. In the implementation, the load balancing is achieved by counting the number of arithmetic operations executed by each processor, and the processor synchronization is achieved by using the Ultra parallel FORTRANcommands "parbegin" and "parend." This gives a computation in parallel. 4.2. Partial Differential Equation Problem Three schemes were employed in the multischeme computation to solve problem (3.3), namely, the Lax, the upwind, and the Lax-Wendroff, respectively. Both Lax and upwind are first order, but Lax-Wendroff is second order. The implementation is arranged as in the ODE problems to be carried out simultaneously both in parallel and in array. However, as shown in Fig. 3, the computation process represents only the one in the ( J + 1)th time interval and seven processors are sufficient to achieve the computation efficiency. For both Lax and Lax-Wendroff, the truncation errors of two-point and three-point computation history are identical to each other; therefore,
T4
T3
T~
T1
-
,
I,EU
Y5
I,EU
I
Tn+21
~Y2 '
3 #~C
Y4 ' ) Y4
l
1Y33,PC
]Y3 ,
Y4
2,PC Y3
Y2
ET4 ,
3,E1T
I
I
I
I
I
[
I
2,EU 2,ME ETn , ETn 3,EU 3 , ~ ETn , ETn
3,p(
2 PC
I
2,RK 2~pC ETn , ET n 3,R-K 3,PC ETn , ETn
3,RK
2,RK
%~-i. %,i-1
I
EU ME Een, Ee n RK PC Een, Een
EU ME Een-l,Een-ll RK PC Een-l,Een-i
Ee4
PC
EU ME Ee3, Ee3 RK PC Ee3, Ee3
Ee2,
RX ~'~ ~ , ~,,~'~ Ee4,
ET3 ,
ET3 ,
P7
EU ME Ee2, Ee 2
2,ME ~T4 2, R/~ 2,PCIEe4, , %4 i
ET4 , ET4
ET3 ,
3,E'U
ET3 ,
I
I
1,ME 2,ME 1,ILK 2,RK 1,PC 2,PC Yn , Yn ETn-I, Yn , Yn Yn ' Yn 3,ME 3,RK 3,PC Yn Yn Yn ETn-I,
1,ME 2,ME I,RK 2,RK I,PC 2,PC Y5 ' Y5 Y5 ' Y5 Y5 , Y5 3,ME 3,RK 3,PC YS Y5 !Y5
Y4
7~
2 , RKt--~, PC
Y2
I
P6
I
P5
7(3)
<3] (2} (3)
(~)
(2) Sn, Yn
3
2
Sn-l" Yn-i
3
Sn-l, Yn-i
2
Sn-2 , Yn-2
3
2 (2) Sn-2, ~'tl- 2
S 3,
3
2 y(21 $2,
I
P8
Fig. 2. The computation process of an ODE initial-valued problem in a multischeme computation of four schemes. Both two-point and three-point switches and the correction factors have been used in the computation.
I
"
2,EU
2,EU
Yn
Yn
Tn§
Y33,RK
Y3 ,
1,RK
Y2 '
I,RK Yl
I,EU 2,ZU I,ME 2,ME Y4 ' Y4 Y4 ' Y4 ' Y4 3,EU 3 , ~ 3,ME Y4 Y4 Y4
I/n
I
Y2
I,EU ~,EU ' 1,ME 2 , M ~ ' Y3 ly2 , Y3 3,EU 3,ME Y~ Y3
u
IY2 '
Y~ '
~ ...... P2
I Yl
Y2
P~
Yl
Tn I 3,E~
--
-
--I
[
9 "O
r~
E
I,UW u2"UW3+3 UJ+3, I,LW 2,LW 3,iR4 UJ+3' UJ+3 UJ+3
I,L 2,L UJ+3, UJ+3
I,L 2,L UJ+4, UJ+4
X3+2
Xj+3
Xj+4
Xj§
X n-J
T, UW E,UW En- J, En- J
T,L E,L En-J, En-j T,LW E,LW En- J , En- J
I
I
Xn-J+2
I
I
I
E
E Sn-j
T J+l Sn_j, Un_ J
I --I
Sn- J -
r"
E Sn -3 -.
E $3~ ]
S3§
I
I
I
P7
T J+l Sn-j- I, Un-j- 1
J+l Un-J-2
J+l Uj§
T Sn-J-2,
T Sj+2,
T 3+I SJ+ 2, UJ+ 2
J§ I,LW UJ+I " UJ+I
I
I
P6
Fig. 3. The computation process in the (J + 1)th time interval of a PDE initial-valued problem in a multischeme computation of three schemes. The switch and the correction factor used in the computation are as follows: The Lax and Lax-Wendroff: the trunedon error of two points computation history. The upwind: the trunction error of three points computation history.
I
I
Xn-J+l
I
T,UW E,UW En-j-l, En-j-I
T,L E,L En-J-1, En-J-l T,LW E,LW En-j- 1 , En-j- l
I,LW 2,LW Un-J' Un-j
I,UW 2,UW Un-J, Un-j 3,D~4 Un- J
I,L U n-J,
2,L U n-J
T, L E, L EJ+4' EJ+4 T,LW E,LW E J+ 4 , EJ+ 4
l, UW u2"UWj+5 UJ+5' I,LW 2 , L W 3,UW Uj+5' U/+5 UJ+5
1,L 2,L UJ+5, U5+5
T,UW E,UW Ej+ 4, Ej+ 4
T,UW E,UW EJ+ 3, EJ+ 3
T,L E,L E3+3, Ej+3 T,LW E,LW EJ+3, E5+3
Uj+4, 3,~'W Uj+4
I,LW 2,LW U/~4, Uj+4
I
I
P5
T,UW E,UW EJ+2, E3+2
I
I
P4
T,L E,L EJ§ EJ+2 T,LW E,LW EJ+ 2 , EJ§
I,LW 2 , L W UJ+2, Uj+ 2
1,UW 2,UW Uj§ 2, Uj+ 2
I,L 2,L UJ+2,, U3+2
Uj+ I
Uj+ 1
P~ 1,LW
P2
1,UW
XJ+l
Pl
I,L UO§
m
ore
J~
Parallel MultischemeComputation
47
only two approximations of different mesh sizes are computed. Of the upwind, the use of the switch of three-point computation history requires an additional approximation. The symbols U I'L and U~'L represent the approximations of the Lax scheme computed at x = x n with the mesh size, Ax and 2Ax, respectively. Correspondingly, the symbols, t7 l'Uw t& 'Uw Ua, uw t71,Lw and U~ 'Lw are defined. The switch S e of the exact truncation errors, E ne,L, E e,uw, and E ne,Lw at x = xn is used to check the correctness of the problem solution U~ +1 chosen by the switch S~r of the estimated errors, E nT,L , EnT, U W and -~,T, n L W , respectively. The occurrence of the problem solution, U~ +1 shows a two-space step delay, and the solution tr s+l at the point x=xj+~ is set as UsJ+l1 = H l'Lw for all J. No boundary ~J+ 1 ~J+ 1 conditions are assumed; therefore, in each advancement of time, the solutions at the boundaries are missing. 5. N U M E R I C A L R E S U L T S The truncation error of computation history is applied also to the computed approximations of the schemes as a solution correction factor. Its application improves the solution accuracy by adding four extra accurate digits. The corrected approximation of a lower-order scheme approaches those of higher orders, and in P D E problems the corrected solution on a coarse grid approaches the ones on fine grids.
5.1. Ordinary Differential Equation Problem Two problems have been studied in detail. In one, the solution slope is given as
f(x, y) = sin 2x + e -y
(5.1)
and the initial condition is y(0) = 0 therefore, an exact solution in closed form is possible and
y(x)= -89
2x-e
x+
3
(5.2)
In the second one,
f(x, y) = y(cos x - 3e x + 1) and
(5.3) y(0) = 1
48
Chiang and Lai
An exact solution in closed form is also possible and
y(x) - -
e sin
x + 3e -x
+ x -- 3
(5.4)
In Fig. 4, the curves are drawn in a uniform scale. In (a), the slope, and the solution, y(x), of problem (5.1) are shown. In (b), the approximations of the schemes were computed without the treatment of a switch and a correction factor and the error curves are separated from each other. As shown, the predictor-corrector method produces the smallest error in magnitude. However, in (c) and (d), after the application of the switch and the correction factor, the errors are considerably reduced in magnitude. In (c), a two-digit solution accuracy improvement is observed and in (d), up to three digits. The most interesting observation is that the errors in the corrected solutions of different schemes approach each other. This shows that the utilization of the truncation error of computation history can be more effective in the improvement of solution accuracy than a scheme variation. The results suggested that the two-term truncation error estimate of finite difference scheme in (2.10) is more suitable than the one in (2.1).
f(x, y),
y = -2-Lcos2x -
e-X + ~2
No Switch and no Correction
..,..."-9
....~ .OOl
.~...f\
.........f ~
......
-.0o~ -.c~z
] ......... y L'-7-o
~y/'d~, = f ( x . y ) ~
2
-.003 3
4
~
6
(a)
(b)
Switch and Correction (2 points)
Switch and Correction (3 points)
....iI
.... iI L
I ......... [ - -
~
* rill dx
=
s
e
=]n2x
+
]
Euler
R u n g e - Kurtcl PredJctor-Corrector
.o0ooa
x 9 -~,
Mo(lJfied
y(0) = o
(c)
x 1
2
3
(d)
Figure 4.
4
5
e
Parallel Multischeme Computation
49
In Fig. 5, the curves are drawn in a semilogarithm scale. In (a), the slope f(x, y) and the solution y(x) of problem (5.3) are shown. In (b), the approximations of the schemes were computed without the treatment of a switch and a correction factor. The error curve of the first-order Euler method is far apart from those of the second order. This agrees with the concept of scheme ordering. However, after the application of the switch and the correction factor not only do the errors of equal-ordered schemes approach each other, as shown in Figure 4, but also the one of a lower order. In this problem, the treatment results in not only a significant solution accuracy improvement, but also a stabilized error distribution. As shown in (d), in the interval (1, 4), the error curves appear as almost straight lines, without the variations as shown in both (b) and (c). This suggests again that the utilization of computation history is more effective in solution accuracy improvement than scheme variation, and the two-term truncation error estimate of (2.10) is more suitable than the one in (2.1). Most importantly, the results raise the question whether high-order scheme computation is necessary.
y lo"
=:
.........
e
sinx
t
3e -*
+
x
-
No Switch and no Correction
3
io 0 -
y
Io'
i ....................
-
9............
~o-'
~~
,I:
ua _=
_1o -~
L-T_ P~"~__=c~"~A=J 2
(a)
3
Switch and correction (2 points)
R~-Kutto Pr~lktor~tof 10 -2
10 -~
w
2
~
~3~dx =
5
S
Switch and Correction (3 points)
1oo
~o ~
4
(b)
4
y (cosx
e
5
-
3e -x
+
1),
(c)
y(o) = 1 (d)
Figure 5.
50
C h i a n g and Lai
5.2. Partial Differential Equation Problem We studied a continuous wave propagation. In problem (3.3), the initial wave is assumed to be
(5.5)
Uo(x ) = 2 sin 2x + cos x and the wave speed, c = 3/2. The exact wave is then
(5.6)
U(x, t) = 2 sin 2(x + 3t) + cos(x - 3t)
Two sets of numerical experiments have been performed with the parameters of mesh size Ax = 1/8 and 1/32, respectively, and a time interval At = 1/64. In both Figs. 6 and 7, the curves are drawn in a uniform scale. In (a), the propagating wave of (5.6) in the first six time intervals is shown. In (b), the error curves of the untreated approximations of the schemes carry wave packets (Courant and Hilbert, 1953) as described in wave propagation theory. However, after the treatment of the switch and the correction
u(~.*) = 2.1,2(x-~)
+ co,(x-,~)
U(x,O) = 2sin2x + cosx The solution is :
il,V,,v,v V,
U(x,t)
= 2sin2(x-2~t ) + c o s ( x - ~ t )
Ax-
1 8
At -
64
(a) No S ~ t c h
Slvitch and Correction
and no Correction
.15
.002
.1 .001
.05
..
/-
:-X
/
...
-
.-
":
-j.
:"*.
.
...
*,'-",
/-
/-
~"-
/" i
-
.1
-.002 X ~za
(c)
(b) F i g u r e 6.
Parallel Multischeme Computation u(x.t) = 2,m.2(,-~)
+
51
~162 3t
-
U(x,0)
2
8x
=
2sin2x
The solution
i::[/ "J
X
+
U(x,t)
= 2sin2(x-~t)
Ax -
1 32
At
64
-
cos•
is : + cos(x-~t)
1
(a) b-~tch and Correction
No Switch and no CorrecUon +
- -
Lox-Wend~(f
....
Upwind
....
up~
ZE-C~
~.02
-2E-O
x
121$
129
130
131
1~2
12e
(b)
129
130
131
(c) Figure 7.
factor, as shown in (c), the wave packets disappear, and the errors are reduced in magnitude. In (c) of Figure 6, the error is reduced to be about only 1/50 its original magnitude, and 1/10,000 in Figure 7. In fact, as shown in (c) of Figure 7, the randomness of the error suggests that the computation error is dominated by the roundoff error, and the truncation error approaches zero. This raises the question whether a fine grid computation is necessary. 6, CONCLUSION The difficulty in the application of a multischeme computation exists in the design of a sensitive mathematical switch. Originally, the multischeme computation was introduced to improve solution accuracy by taking advantage of the irregularity occurring in the errors of various schemes. However, this approach has its limitation. Our results on solution accuracy improvement are out of this sequel. They are mainly the contribution of using the truncation error of computation history as a solution correction factor. This achievement should not devalue the importance of a multischeme computation because only in a multischeme computation
52
Chiang and Lai
environment can the schemes be compared and studied in an equal footing. Our discoveries cannot be observed if only single-scheme computations are performed. Based on our results, a new direction in scientific computing may be necessary. (a) The application of our correction factor achieves a solution accuracy improvement of four extra digits in a single precision computation. The error in the corrected solution of a three-point correction factor appears more stable than in one of a two-point. This suggests that the two-term truncation error estimates of finite difference schemes as given in (2.10) and (2.11) are more suitable than the ones in (2.1) and (2.12), respectively. (b) In ODE problems, not only are the corrected solutions of equally ordered schemes approaching each other, but also the one of a lower-order scheme. This suggests that in solution accuracy improvement, the utilization of computation history may play a more important role than scheme variation. It suggests also that a high-order scheme computation may be unnecessary. (c) In a PDE initial-valued problem, the wave packet disappears in the error curve of a corrected solution; therefore, the corrected solutions are more stable than those without a correction. Furthermore, the errors in the corrected solution are also considerably reduced in magnitude; therefore, a fine grid computation may not be needed either. (d) For both Lax and Lax-Wendroff schemes, the truncation errors of both two-point and three-point computation history are identical, but not so for the upwind scheme. What is the reason for this? (e) The assumption on the derivation of the truncation error of computation history is that at the computing step, the nth derivative of the problem solution exists for all n. It cannot be held true at the solution discontinuities. Is the correction factor still applicable to the scheme solution in the neighborhood of a discontinuity to improve solution accuracy and also to depress discontinuity smearing? (f) Are there other applications of the computation history of a scheme?
REFERENCES Ames, W. F. (1976). Numerical Method for Partial Differential Equations, Academic Press, New York. Chiang, Y. F. (1986). Use of mathematical switch to solve differential equation problems, ACM SIGNUM Newsletter, October, pp. 20-33.
Parallel Multischeme Computation
53
Chiang, Y. F, Ma, J. S., Hu, K. L., and Chang, C. Y. (1988). Parallel multischeme computation, J. Sei. Comput. 3(3), 289-306. Chiang, Y. F. (1988). Truncation and accumulated errors in wave propagation, J. Comput. Phys. 79(2), 353-372. Chiang, Y. F. (1984). Pointwise work convergence for solving linear hyperbolic system with discontinuous initial data, Research Report No. 6, CIS, New Jersey Institute of Technology. Courant, R., and Hilbert, D. (1953). Methods of Mathematical Physics, Interscience Publishers, New York. Hamming, R. W. (1962). Numerical Methods for Scientists and Engineers, McGraw-Hill, New York. Morris, J. LI. (1985). Computational Method in Elementary Numerical Analysis, Wiley-Interscience, New York. Petzoid, L. (1984). Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations, SIAM J. Sei. Stat. Comp. 4(1), 136-148. Richter, G. R., and Falk, R. S. (1984). An analysis of a finite element methods for hyperbolic equations, Adv. Comput. Methods Partial Differential Equations V, 297-301. Richtmyer, R. D., and Morton, K. W. (1969). Difference Methods for Initial Value Problems. Interscience Publishers, New York. Roache, P. J. (1977). Computational Fluid Dynamics, Hermosa Publishers. Ultra Computer (1987). Documentation Note No. 2: How to write parallel programs for the NYU Ultra Computer prototype. Vichnevetsky, R. (1981). Computer Methods for Partial Differential Equations, Prentice-Hall, Englewood Cliffs, New Jersey. Vichnevetsky, R., and Bowles, J. B. (1982). Fourier Analysis of Numerical Approximations of Hyperbolic Equations, SIAM, Philadelphia.
854/5/1-4