- Email: [email protected]
Mesh refinements for parabolic equations
JOURNAL
OF CQMPUTAT1ONAL
12, 513-525
PHYSICS
(1973)
Mesh Refinements for Parabolic Equations* MELVYN CIMENT+ Mathematical Analysis Division, Naval Ordnance Laboratory, White Oak, Silver Spring, Maryland 20910
AND ROLAND A. SWEET Department of Mathematics, University of Colorado, Denver Center, Colorado Denver, 80203 Received December 28, 1972
Consider a finite difference approximation of the Crank-Nicolson type, using a local mesh refinement to an initial-boundary value problem for the parabolic equation m(X)& = wbd,
+ f 6, t).
We consider the above difference approximation to u(n, t) on a grid pattern, where the middle third of the unit interval is refined with respect to the rest of the interval. Using the standard matrix arguments, we prove the matrix stability (spectral radius less than one) for the above class of problems. The local mesh refinement is achieved without loss of the second-order spatial accuracy, or of computational efficiency in solving the resulting equations. In two dimensions, we use our techniques to adapt the alternating direction schemes to a general mesh refinement pattern. Numerical experiments are presented which illustrate the computational effectiveness of the mesh refinement technique.
1. INTR~OUC~~N
Mesh refinement schemesfor the numerical solution of initial-boundary value problems of parabolic type were proposed by Ciment in [l]. After obtaining very general results on difference shemes for initial-boundary value problems, Osher [3] and Varah [6] used the mesh refinement representation developed by Ciment in [2] to show that these parabolic methods are convergent. However, * This research was partially supported by the Marathon Oil Company. + This research supported in part by the Oflice of Naval Research, contract number NR 044-377.
513 Copyright All rights
Q 1973 by Academic Press. Inc. of reproduction in any form reserved.
581/n/4-6
514
CIMENT AND SWEET
their proofs were for implicit schemeson an infinite interval, a noncomputable case,since it gives rise to a doubly infinite matrix system. In this paper we consider the actual computational problem that one faces in performing the implicit mesh refinement calculations on a finite interval. Failing to prove convergence, we present a proof of the matrix stability (eigenvalues all less than one in absolute value [7]) of the associatediteration matrix. In two spacedimensions, we show how to adapt the alternating direction schemes of Peacemanand Rachford [4] to a general mesh refinement pattern. For this case, we present the results of numerical experiments which indicate under what circumstances the various methods may be useful. Some of the computations in this paper were carried out by J. P. Brzozowski of the Denver Research Center of Marathon Oil Company, and we welcome this opportunity to express our thanks to him.
2. MESH REFINEMENT SCHEMES
We are interested in approximating the solution of the following initial-boundary value problem. Consider
on the semi-inhnite strip {(x, t) : 0 < x < 1, t > 0}, where m, S, and f are given functions satisfying
(9 44 ECP, 11
and
(ii)
and
S(X)E Cl[O, l]
pin, .I
m(x) = g > 0,
,~I& s(x) = 2 > 0. . .
Without loss of generality we take for boundary conditions u(0, t) = u(l, t) = 0
0 > 01,
(1’4
and assumethe initial condition is given by
4x, 0) = v,(x)
(0 d x < l),
(14
9 being at least continuous on [0, I]. We approximate u(x, t) with a grid using a uniform time step, but with two different spatial mesh patterns. More precisely, let h > 0 and k > 0 be given, and let A4 > 1 be a positive integer.
515
MESH REFINEMENTS
We define the grid R by t, = nk,
I
n = 0, 1, 2,...
O)L M’
p+l
q+l
Note that p, q, and N depend on the refinement desired. As an example, for h = &, and refining the grid in [$, $1 by a factor of M = 3 (hence,p = 2 and q = S), we would have the grid depicted as
5 6
2 --7 8 3 --1011 4 6 18 18 6 18 18 6
1
0
6
1
We define v~,, as the difference approximation to U(X~, t,J which solves the following standard system of equations vi.0 = d%> %(%?I+,
-
%?a)
=
w~+I,n+ll + (1 - @%1VM.n- (%-1+ 4 0j.n+ Wi+1,nl~ + @ia, Pa)
~ceh-l%l,n+l
-
6-l
+
Si> Vj,n+l
+
where h = (k/h2), mi = m(xj), si = s(xf + (h/2)), and fjsn = f(xj , t,J, for l,Cj,
516
CIMENT
AND SWEET
In our numerical experiments, we found that the third approximation is generally the most accurate method. In the following, we shall provide a proof which indicates the matrix stability [7] of our mesh refinement problem, using the approximation of type (iii). Our proof will clearly hold for the approximation of type (i). For type (ii), however, no such proof was found. For type (iii) we are able to maintain an overall accuracy of O(P), and we need very little extra computational effort, as will be demonstrated later. Hence, we take the remaining two difference equations as
m9(v9.n+l - 44 = WLb-l.n+l +
(1
-
- (kl + s9)v,.~+~ + w,+M.~+J
@b9-lV9-1.n
-
69-l
+
s9)
v9.n
+ S9V9+AI,nlI
+ Win
P-9
and mq(vq.n+l-
v q,n1 =
%%+I,n+ll + (1 - @bq-lvq-M,~ - (47-l+ %>van+ %%+l.nl~ (24 + k&n -
4&-1%~ma+1
-
6%-l
+
%I
%,n+1
+
We now study the stability of this family of difference schemes.To accomplish this, let us rewrite the difference schemein matrix notation as v, = @ (3)
(D + AdJ) vn+1 = [D - X(1 - B)J] V, + kF,, , where ),
@f%j
D = diag{m, , m2 ,..., mN}, and J is an almost tridiagonal matrix. In fact, its jth 1 and q + 1 < row is {O,..., 0, --sj-1, q-1 + sj, --s, 3 O,...,0}, for 1 < j d p j < N; and (0,..., 0, -M?sj~l ) Mysj_l + Sj), -M%, , 0)..., O},
for p-j-1
Qj
Each of the two remaining rows corresponding to the interfaces has three nonzero elements in them, one on the main diagonal, one beside it, and one a distance M away. Row p is (0,...) 0, -s,_1 ) s,-1 + s, , 0,...) 0, -s, ) 0,...) O},
517
MESH REFINEMENTS
and row q is
(0)...) 0, -s,-1 , 0)..., 0, s,-1 + s, , -s, , o,..., 01. It is easily seen that D-l.J is irreducible, and diagonally dominant (Varga [7]). Hence, by an extension of Gerschgorin’s theorem, its eigenvalues 7i satisfy Re(Ti) > 0,
lfi
Thus, the matrix D + XOJ= D(I + XBD-lJ) is nonsingular, insuring that V,,, is always uniquely detied. Solving for V,,, in (3), we obtain the iteration matrix C = (I + ;\kD-‘J)-l[l - A(1 - 8) D-lJ]. To prove convergence, one needs to show that ((Ck (I is uniformly bounded [5] as the grid sizes tend to zero. For the cases 8 = 0, 1, the standard maximum principle arguments yield such proofs. However, for general 8, we were unable to prove convergence. An obviously necessary condition for convergence, referred to as matrix stability (Varga [7]), requires that the spectral radius of C be less than or equal to one. This latter (weaker) type of stability holds for (2). THEOREM. The mesh refinement scheme(2) gives rise to a stable iteration matrix
when: (a) 0 < e < a, if
where S = ~II..~ s(x),
(b) ~~O~l,f0rdlh>O. Proof. The eigenvalues pi of C are given by
CL. - 0)Ti 2= 1- i w+ xeTi 3
1
Let 7 = 01+ $3, 01,/I real, be an eigenvalue of D-‘J. We know that 01> 0. By simple algebra we get, using (4) 1 > j pi/z=
p - h(i - e)q + [h(i [i + xe42 + wgi2
is true if and only if -2o1x + xyi - 2ep is true.
+ fly < 0
e)jp
518
CIMENT AND SWEET
Since (Y> 0, for 4 < 0 < 1, the above inequality is satisfied for all h > 0. When 0 < 0 < 4, the above inequality is true if
An upper bound on A, which insures our inequality, may be determined by finding the minimum of
where 01+ ip lies in the union of the Gershgorin discs of D-% 9j has the form
A typical disc
sj = {z : / z - rj j < rj , where rj = (P’J)~~}.
(Ol
-
rj)2
+
p”
<
rj2.
Hence, 1 d% PI d a2 + rj2 - (a - rJ2 = 2rj ’ 01
so
Finally, if we take h < 2(1 - %) SM2 ’ then Ipi/ < 1.
3. SOLUTION OF EQUATIONS
For 8 > 0, the solution of (3) at each time step involves solving a linear system of equations of the form
where A is almost a tridiagonal matrix and B is a vector formed from known quantities. Solving this system is best accomplished by initially factoring A as the
519
MESH REFINEMENTS
product of a lower unit triangular matrix L, and an upper triangular matrix R (Wilkinson [S]), and then, at each time step, solving two triangular systems of equations. Since A is nearly tridiagonal, one would expect L and R to have nearly the same form-namely, bidiagonal-as when A is tridiagonal. We show now that very few nonzero elements are added in L and R, due to the form of A. Let A = (Q) be a N x N matrix where aij = 0, if 1i - j 1 > 1, except aB,D+m# 0 and u*,~-~ # 0. Letting L = (tij) and R = (rii), one easily verifies that the nonzero elements of L are eiii = 1
(i = 1, 2,..., N)
cz,i-l= -.Lli-l,i-1
(i = 2,..., N-
G.i-1
1; i # q)
1 e&,-M
=
aqsq-M
ra-M,q-M
G9,q-i
=
&
(aq.q-i
-
e9p9-i-lr9-i-l,9-i)
(i = M - 1, M - 2,..., l),
and the nonzero elements of R are r I,1 = %l rii = Uii - ti,i-,ri-, ri-1.i = ai-l,i r n,v+M = rp+i,p,+M
=
(i = 2,..., N) (i = 2, 3,..., N;i#p+m)
i
aB,efM avii.D+M
-
~~+i.P+i-lrp+i-l,p+iu
(i = 1, 2,..., M - 1).
An operation count (one operation is either a multiplication or division) shows that a total of 2(N - 1) + 3(M - 1) operations are required to compute L and R, initially. The term 3(M - 1) is due to the “extra” elements of A. To solve the two tirangular systemsat each time step, 3N - 2 + 2(M - 1) operations are required, where the factor 2(M - 1) is, again, from the ,,extra” elements of A. Clearly, M will be small compared to N, and so, very little additional work is required to solve the implicit system of equations. 4. HIGHER DIMENSIONS Consider the model initial-boundary value problem Ut = u,, + u,, + f(x, Y, 0;
dx, Y, 0) = @,(x,Y>,
0
& 1
t>o (6)
520
CIMENT
I
‘+ 6+
II+
AND SWEET
2+ 7
3 +
+
4
30
31
8
32
33
9
34
.
.
.
.
.
40
41
1:
.
.
.
.
47
12 i
.
.
.
.
17 .
16
46
.
.
.
.
.
5.q
55
.
.
.
.
.
61
+
5
P
+t9
20+ 1
‘;rj-+26
2:
AxI
‘Ay,.;
+2e
Ax*=
Ay2=~
+ 29
1
FIGURE 1
with the solution prescribed on the sides of the square. We treat the situation of refining the mesh in a region of interest (seeFig. 1) while maintaining a uniformly coarser grid elsewhere. The difficulty usually encountered is a sharp increase in the complexity of the equations to be solved. Our method avoids this by adapting the alternating direction (ADI) scheme of Peaceman and Rachford [4] to our general mesh refinement pattern. The ADI equations are a”,
, k)n+1’2
+
+
+J;(;k
At ur”$’ - v;;1/2 = ; [(83Jk) %+1/Z + (Sy2vj,k)fi+l]+ -#l/Z, where h = At/Ax2, and = vxl& - 2v;f, + &&. With ADI, it will be seen that our one-dimensional method becomesthe basic component for higher dimensional mesh refinement problems. We were unable to prove the matrix stability of this procedure, but numerical experiments presented in Section 5 indicate that the method is accurate. (Sx2vj,k)n
MESH REFINEMENTS
521
For simplicity, we consider a mesh pattern where the region
D, = ((x, Y) : f f x < 3,i < Y d $1
(8)
has different (finer) mesh spacings than the remaining region. Namely, Ax,=Ay,=$inD,;
Ax, = Ay, = i outside D, .
The way one numbers (orders) the mesh points (unknowns) can significantly alter the easewith which one can solve the resulting equations. Figure 1 exhibits the resulting pattern, containing 65 points and their ordering. We now describe the approximation for the mesh pattern of Fig. 1. The generalizations to more complicated patterns will be straightforward and, thus, will be omitted from our discussion. Initially, say t = IZAt, we solve Eq. (7) for the grid values at the points l-29. Note that to obtain these values on any one line, we solve the associated system of equations along that line as in a one-dimensional implicit method of Section 2. For this example, on the lines (l-5}, {6-lo}, (20-24}, (2529}, we obtain the standard tridiagonal systems. In this simple pattern only one “coarse” line cuts through D, , namely, points (11-19). The system for this line can be solved by the one-dimensional methods of Section 2. Among values at (l-29}, we have also obtained the values, at time t = (n + 4) At, along all the “coarse” points of D, . The data at the points (7, 12, 21, 9, 18, 23) allows us to use any desired degree of interpolation to obtain the values at the points (34, 41, 48, 55, 40, 47, 54, 61). Now the x sweep inside D, is carried out, using uniform alternating direction on (34-61). Completing this, we now interpolate in the data on the top and bottom of D, , at points (30-33) and (62-65). Note that we describe this interpolation at the end to emphasize that in doing the x sweep in D, , one uses the data from t = n At on the top and bottom of D, in the explicit part of (7). This completes the x sweep, and all our grid values are at t = (n + 4.)At. The y sweep is carried out analogously, and completes the full time step. In the previous example, the method on each line is the uniform grid method, except on the lines cutting through D, . To investigate which condition is best to use at the interface points, we used all interface conditions listed in Section 2 in several computational experiments. V. COMPUTATIONAL RESULTS We solved Eq. (la) on the interval [0,4], using three different grid regions. The true solution U(X,t) = exp[x( -(37/6) + (13/6) x) + t] grows very rapidly near the right side. Starting with a Ax, = .20, we refined, at x, = 3.0, by a factor of 10,
522
FIG. A
CIMENT
2A.
Crank-Nicolson
AND
SWEET
Percent Error, IO:1 Mesh Ratio at X = 3.0, 20 time steps.
4PT; 03PT; DINT.
To.12
FIG.
2B. Crank-Nicolson Percent Error, 4:l Mesh Ratio at X = 3.5. A 4PT;
0
3PT; 0 INT.
523
MESH REFINEMENTS
so that dx, = -02, and again refined at about X, = 3.5, down to dx, = XXX, and dt = .Ol. The computing times for all methods were within a 2 y0 difference of each other. The results shown in Figs 2A and 2B are characteristic for all our calculations. Namely, the three-point method causes an abrupt oscillation of the approximate solution, while the other two methods are smoothly varying across the refinement interface. For smooth problems, even though this overshoot creates a larger error for the three-point method, it will eventually be damped out, and can actually causean improvement of the error profile after many time steps. However, in most situations the refined area is usually changed after several time steps, and in such casesone would expect that the three-point method might introduce many oscillations. TABLE I Ax, = Ay, = l/15; Ax, = Ay, = l/90; PI = P, = 8; At = 2 x 10-k Alternating direction Mesh refinement (a) Integer-fit” (b) 4-point (c) 3-point
t, time
IOAr 1OAt
1OAt
EC
Eint
3.7 x 30-2 3.1 x 10-Z 5.5 x 10-z
7.6 x lO-s 7.8 x 10-B 1.2 x 10-Z
a For the entries in (a), ECis a 5 % relative error while Eint is 0.12 %.
In Table I, we list the results of solving (6) by using our method, with the three interface conditions of Section 2. We took the analytical solution as
4~ Y,t>= (t +! a)p . ew i (x -
1/2y
+ (u -
-4(t + a)
w
1’
with j3 = 40 and OL= l/432. To describe the accuracy of the refined scheme away from the interface, we compare the error in the interior of D, at points at a distance of at least, a coarse mesh-width away from the boundary of D, . The maximum deviation of our approximation from the solution over all such points we denote by Eint . The maximum error over the remaining points in D, and over all the coarse points, we denote by E, . Note that, for all our calculations, we always used third-order Lagrange interpolation along the interface. For the time interval indicated in Table I, this function exhibits large gradients near D, . In fact, at t = 0, u(x., y, 0) has a maximum value of 10.8, and U, , U,
524
CIMENT AND SWEET
take on values greater than 80 in absolute value in D, . However, the region of large gradients begins to spread out, and our method remains an effective computational procedure only if we have a dynamically changing grid. In this case, the relative errors reflect the effectivenessof the refinement technique. TABLE II Integer-fit: 10 time steps; At = 0.006 -5 1. No refinements Ax = Ay = l/l5 2. No refinement Ax = Ay = l/60 3. Refinement Ax, = Ay, = l/l5 Pi = P2 = 4
&t
Work estimate
0.004608
0.010235
“(15)?
0.000277
0.000671
h16(15)2
0.00506
0.00201
~2.77(15)~
In Table II, we show the results for a case where the region of sharp gradients is more static. We solved for U(X,y, t) = 1 + exp[--2Ot[(x - i)>” + (y - $))“I] a solution of (6) on D, as in (8), using the integer-# method. With the mesh refinement, we achieve an increase of accuracy by a factor of, at least, five over the interior square (much more away from the interface), by using only 2.77 as many operations. As shown in Table II, the error, with no refinement, is inversely proportional to the number of grid points. Thus, the only other way to have obtained this increased accuracy in the interior region would have been to use a uniform grid
Bx=By=l--.
15 2/5
This would have required nearly twice as much storage and calculation as compared to the mesh refinement scheme. REFERENCES 1. M. CIMENT, Mesh Refinement Techniques-Fundamental Considerations, Internal Report No. 69-5, Marathon Oil Company, Denver Research Center, January, 1969. 2. M. CIMENT, Stable difference schemes with uneven mesh spacings, A&h. Comp. 25 (1971), 219-227. 3. S. OSHER, Mesh refinements for the heat equation, SIAM J. Numer. Anal. 7 (1970), 199-205. 4. D. W. PEACEMAN AND H. H. RACHFORD, JR., The numerical solution of parabolic and elliptic differential equations, SIAM J. 3 (1955), 28-41.
MESH REFINEMENTS
525
5. R. D. RICHTM~ERAND K. W. MORTON, “Difference Methods for Initial-Value Problems,” 2nd Edition, Interscience, New York, 1967. 6. J. VARAH, Stability of difference approximations to the mixed initial boundary value problem for parabolic systems, SIAM J. Namer. Anal. 8 (1971), 598-615. 7. R. S. VARGA, “Matrix Iterative Analysis,” Prentice-Hall, New York, 1962. 8. 3. H. WILRINSON,“The Algebraic Eigenvalue Problem,” Oxford University Press, London/ New York, 1965.
OF CQMPUTAT1ONAL
12, 513-525
PHYSICS
(1973)
Mesh Refinements for Parabolic Equations* MELVYN CIMENT+ Mathematical Analysis Division, Naval Ordnance Laboratory, White Oak, Silver Spring, Maryland 20910
AND ROLAND A. SWEET Department of Mathematics, University of Colorado, Denver Center, Colorado Denver, 80203 Received December 28, 1972
Consider a finite difference approximation of the Crank-Nicolson type, using a local mesh refinement to an initial-boundary value problem for the parabolic equation m(X)& = wbd,
+ f 6, t).
We consider the above difference approximation to u(n, t) on a grid pattern, where the middle third of the unit interval is refined with respect to the rest of the interval. Using the standard matrix arguments, we prove the matrix stability (spectral radius less than one) for the above class of problems. The local mesh refinement is achieved without loss of the second-order spatial accuracy, or of computational efficiency in solving the resulting equations. In two dimensions, we use our techniques to adapt the alternating direction schemes to a general mesh refinement pattern. Numerical experiments are presented which illustrate the computational effectiveness of the mesh refinement technique.
1. INTR~OUC~~N
Mesh refinement schemesfor the numerical solution of initial-boundary value problems of parabolic type were proposed by Ciment in [l]. After obtaining very general results on difference shemes for initial-boundary value problems, Osher [3] and Varah [6] used the mesh refinement representation developed by Ciment in [2] to show that these parabolic methods are convergent. However, * This research was partially supported by the Marathon Oil Company. + This research supported in part by the Oflice of Naval Research, contract number NR 044-377.
513 Copyright All rights
Q 1973 by Academic Press. Inc. of reproduction in any form reserved.
581/n/4-6
514
CIMENT AND SWEET
their proofs were for implicit schemeson an infinite interval, a noncomputable case,since it gives rise to a doubly infinite matrix system. In this paper we consider the actual computational problem that one faces in performing the implicit mesh refinement calculations on a finite interval. Failing to prove convergence, we present a proof of the matrix stability (eigenvalues all less than one in absolute value [7]) of the associatediteration matrix. In two spacedimensions, we show how to adapt the alternating direction schemes of Peacemanand Rachford [4] to a general mesh refinement pattern. For this case, we present the results of numerical experiments which indicate under what circumstances the various methods may be useful. Some of the computations in this paper were carried out by J. P. Brzozowski of the Denver Research Center of Marathon Oil Company, and we welcome this opportunity to express our thanks to him.
2. MESH REFINEMENT SCHEMES
We are interested in approximating the solution of the following initial-boundary value problem. Consider
on the semi-inhnite strip {(x, t) : 0 < x < 1, t > 0}, where m, S, and f are given functions satisfying
(9 44 ECP, 11
and
(ii)
and
S(X)E Cl[O, l]
pin, .I
m(x) = g > 0,
,~I& s(x) = 2 > 0. . .
Without loss of generality we take for boundary conditions u(0, t) = u(l, t) = 0
0 > 01,
(1’4
and assumethe initial condition is given by
4x, 0) = v,(x)
(0 d x < l),
(14
9 being at least continuous on [0, I]. We approximate u(x, t) with a grid using a uniform time step, but with two different spatial mesh patterns. More precisely, let h > 0 and k > 0 be given, and let A4 > 1 be a positive integer.
515
MESH REFINEMENTS
We define the grid R by t, = nk,
I
n = 0, 1, 2,...
O)L M’
p+l
q+l
Note that p, q, and N depend on the refinement desired. As an example, for h = &, and refining the grid in [$, $1 by a factor of M = 3 (hence,p = 2 and q = S), we would have the grid depicted as
5 6
2 --7 8 3 --1011 4 6 18 18 6 18 18 6
1
0
6
1
We define v~,, as the difference approximation to U(X~, t,J which solves the following standard system of equations vi.0 = d%> %(%?I+,
-
%?a)
=
w~+I,n+ll + (1 - @%1VM.n- (%-1+ 4 0j.n+ Wi+1,nl~ + @ia, Pa)
~ceh-l%l,n+l
-
6-l
+
Si> Vj,n+l
+
where h = (k/h2), mi = m(xj), si = s(xf + (h/2)), and fjsn = f(xj , t,J, for l,Cj,
516
CIMENT
AND SWEET
In our numerical experiments, we found that the third approximation is generally the most accurate method. In the following, we shall provide a proof which indicates the matrix stability [7] of our mesh refinement problem, using the approximation of type (iii). Our proof will clearly hold for the approximation of type (i). For type (ii), however, no such proof was found. For type (iii) we are able to maintain an overall accuracy of O(P), and we need very little extra computational effort, as will be demonstrated later. Hence, we take the remaining two difference equations as
m9(v9.n+l - 44 = WLb-l.n+l +
(1
-
- (kl + s9)v,.~+~ + w,+M.~+J
@b9-lV9-1.n
-
69-l
+
s9)
v9.n
+ S9V9+AI,nlI
+ Win
P-9
and mq(vq.n+l-
v q,n1 =
%%+I,n+ll + (1 - @bq-lvq-M,~ - (47-l+ %>van+ %%+l.nl~ (24 + k&n -
4&-1%~ma+1
-
6%-l
+
%I
%,n+1
+
We now study the stability of this family of difference schemes.To accomplish this, let us rewrite the difference schemein matrix notation as v, = @ (3)
(D + AdJ) vn+1 = [D - X(1 - B)J] V, + kF,, , where ),
@f%j
D = diag{m, , m2 ,..., mN}, and J is an almost tridiagonal matrix. In fact, its jth 1 and q + 1 < row is {O,..., 0, --sj-1, q-1 + sj, --s, 3 O,...,0}, for 1 < j d p j < N; and (0,..., 0, -M?sj~l ) Mysj_l + Sj), -M%, , 0)..., O},
for p-j-1
Qj
Each of the two remaining rows corresponding to the interfaces has three nonzero elements in them, one on the main diagonal, one beside it, and one a distance M away. Row p is (0,...) 0, -s,_1 ) s,-1 + s, , 0,...) 0, -s, ) 0,...) O},
517
MESH REFINEMENTS
and row q is
(0)...) 0, -s,-1 , 0)..., 0, s,-1 + s, , -s, , o,..., 01. It is easily seen that D-l.J is irreducible, and diagonally dominant (Varga [7]). Hence, by an extension of Gerschgorin’s theorem, its eigenvalues 7i satisfy Re(Ti) > 0,
lfi
Thus, the matrix D + XOJ= D(I + XBD-lJ) is nonsingular, insuring that V,,, is always uniquely detied. Solving for V,,, in (3), we obtain the iteration matrix C = (I + ;\kD-‘J)-l[l - A(1 - 8) D-lJ]. To prove convergence, one needs to show that ((Ck (I is uniformly bounded [5] as the grid sizes tend to zero. For the cases 8 = 0, 1, the standard maximum principle arguments yield such proofs. However, for general 8, we were unable to prove convergence. An obviously necessary condition for convergence, referred to as matrix stability (Varga [7]), requires that the spectral radius of C be less than or equal to one. This latter (weaker) type of stability holds for (2). THEOREM. The mesh refinement scheme(2) gives rise to a stable iteration matrix
when: (a) 0 < e < a, if
where S = ~II..~ s(x),
(b) ~~O~l,f0rdlh>O. Proof. The eigenvalues pi of C are given by
CL. - 0)Ti 2= 1- i w+ xeTi 3
1
Let 7 = 01+ $3, 01,/I real, be an eigenvalue of D-‘J. We know that 01> 0. By simple algebra we get, using (4) 1 > j pi/z=
p - h(i - e)q + [h(i [i + xe42 + wgi2
is true if and only if -2o1x + xyi - 2ep is true.
+ fly < 0
e)jp
518
CIMENT AND SWEET
Since (Y> 0, for 4 < 0 < 1, the above inequality is satisfied for all h > 0. When 0 < 0 < 4, the above inequality is true if
An upper bound on A, which insures our inequality, may be determined by finding the minimum of
where 01+ ip lies in the union of the Gershgorin discs of D-% 9j has the form
A typical disc
sj = {z : / z - rj j < rj , where rj = (P’J)~~}.
(Ol
-
rj)2
+
p”
<
rj2.
Hence, 1 d% PI d a2 + rj2 - (a - rJ2 = 2rj ’ 01
so
Finally, if we take h < 2(1 - %) SM2 ’ then Ipi/ < 1.
3. SOLUTION OF EQUATIONS
For 8 > 0, the solution of (3) at each time step involves solving a linear system of equations of the form
where A is almost a tridiagonal matrix and B is a vector formed from known quantities. Solving this system is best accomplished by initially factoring A as the
519
MESH REFINEMENTS
product of a lower unit triangular matrix L, and an upper triangular matrix R (Wilkinson [S]), and then, at each time step, solving two triangular systems of equations. Since A is nearly tridiagonal, one would expect L and R to have nearly the same form-namely, bidiagonal-as when A is tridiagonal. We show now that very few nonzero elements are added in L and R, due to the form of A. Let A = (Q) be a N x N matrix where aij = 0, if 1i - j 1 > 1, except aB,D+m# 0 and u*,~-~ # 0. Letting L = (tij) and R = (rii), one easily verifies that the nonzero elements of L are eiii = 1
(i = 1, 2,..., N)
cz,i-l= -.Lli-l,i-1
(i = 2,..., N-
G.i-1
1; i # q)
1 e&,-M
=
aqsq-M
ra-M,q-M
G9,q-i
=
&
(aq.q-i
-
e9p9-i-lr9-i-l,9-i)
(i = M - 1, M - 2,..., l),
and the nonzero elements of R are r I,1 = %l rii = Uii - ti,i-,ri-, ri-1.i = ai-l,i r n,v+M = rp+i,p,+M
=
(i = 2,..., N) (i = 2, 3,..., N;i#p+m)
i
aB,efM avii.D+M
-
~~+i.P+i-lrp+i-l,p+iu
(i = 1, 2,..., M - 1).
An operation count (one operation is either a multiplication or division) shows that a total of 2(N - 1) + 3(M - 1) operations are required to compute L and R, initially. The term 3(M - 1) is due to the “extra” elements of A. To solve the two tirangular systemsat each time step, 3N - 2 + 2(M - 1) operations are required, where the factor 2(M - 1) is, again, from the ,,extra” elements of A. Clearly, M will be small compared to N, and so, very little additional work is required to solve the implicit system of equations. 4. HIGHER DIMENSIONS Consider the model initial-boundary value problem Ut = u,, + u,, + f(x, Y, 0;
dx, Y, 0) = @,(x,Y>,
0
& 1
t>o (6)
520
CIMENT
I
‘+ 6+
II+
AND SWEET
2+ 7
3 +
+
4
30
31
8
32
33
9
34
.
.
.
.
.
40
41
1:
.
.
.
.
47
12 i
.
.
.
.
17 .
16
46
.
.
.
.
.
5.q
55
.
.
.
.
.
61
+
5
P
+t9
20+ 1
‘;rj-+26
2:
AxI
‘Ay,.;
+2e
Ax*=
Ay2=~
+ 29
1
FIGURE 1
with the solution prescribed on the sides of the square. We treat the situation of refining the mesh in a region of interest (seeFig. 1) while maintaining a uniformly coarser grid elsewhere. The difficulty usually encountered is a sharp increase in the complexity of the equations to be solved. Our method avoids this by adapting the alternating direction (ADI) scheme of Peaceman and Rachford [4] to our general mesh refinement pattern. The ADI equations are a”,
, k)n+1’2
+
+
+J;(;k
At ur”$’ - v;;1/2 = ; [(83Jk) %+1/Z + (Sy2vj,k)fi+l]+ -#l/Z, where h = At/Ax2, and = vxl& - 2v;f, + &&. With ADI, it will be seen that our one-dimensional method becomesthe basic component for higher dimensional mesh refinement problems. We were unable to prove the matrix stability of this procedure, but numerical experiments presented in Section 5 indicate that the method is accurate. (Sx2vj,k)n
MESH REFINEMENTS
521
For simplicity, we consider a mesh pattern where the region
D, = ((x, Y) : f f x < 3,i < Y d $1
(8)
has different (finer) mesh spacings than the remaining region. Namely, Ax,=Ay,=$inD,;
Ax, = Ay, = i outside D, .
The way one numbers (orders) the mesh points (unknowns) can significantly alter the easewith which one can solve the resulting equations. Figure 1 exhibits the resulting pattern, containing 65 points and their ordering. We now describe the approximation for the mesh pattern of Fig. 1. The generalizations to more complicated patterns will be straightforward and, thus, will be omitted from our discussion. Initially, say t = IZAt, we solve Eq. (7) for the grid values at the points l-29. Note that to obtain these values on any one line, we solve the associated system of equations along that line as in a one-dimensional implicit method of Section 2. For this example, on the lines (l-5}, {6-lo}, (20-24}, (2529}, we obtain the standard tridiagonal systems. In this simple pattern only one “coarse” line cuts through D, , namely, points (11-19). The system for this line can be solved by the one-dimensional methods of Section 2. Among values at (l-29}, we have also obtained the values, at time t = (n + 4) At, along all the “coarse” points of D, . The data at the points (7, 12, 21, 9, 18, 23) allows us to use any desired degree of interpolation to obtain the values at the points (34, 41, 48, 55, 40, 47, 54, 61). Now the x sweep inside D, is carried out, using uniform alternating direction on (34-61). Completing this, we now interpolate in the data on the top and bottom of D, , at points (30-33) and (62-65). Note that we describe this interpolation at the end to emphasize that in doing the x sweep in D, , one uses the data from t = n At on the top and bottom of D, in the explicit part of (7). This completes the x sweep, and all our grid values are at t = (n + 4.)At. The y sweep is carried out analogously, and completes the full time step. In the previous example, the method on each line is the uniform grid method, except on the lines cutting through D, . To investigate which condition is best to use at the interface points, we used all interface conditions listed in Section 2 in several computational experiments. V. COMPUTATIONAL RESULTS We solved Eq. (la) on the interval [0,4], using three different grid regions. The true solution U(X,t) = exp[x( -(37/6) + (13/6) x) + t] grows very rapidly near the right side. Starting with a Ax, = .20, we refined, at x, = 3.0, by a factor of 10,
522
FIG. A
CIMENT
2A.
Crank-Nicolson
AND
SWEET
Percent Error, IO:1 Mesh Ratio at X = 3.0, 20 time steps.
4PT; 03PT; DINT.
To.12
FIG.
2B. Crank-Nicolson Percent Error, 4:l Mesh Ratio at X = 3.5. A 4PT;
0
3PT; 0 INT.
523
MESH REFINEMENTS
so that dx, = -02, and again refined at about X, = 3.5, down to dx, = XXX, and dt = .Ol. The computing times for all methods were within a 2 y0 difference of each other. The results shown in Figs 2A and 2B are characteristic for all our calculations. Namely, the three-point method causes an abrupt oscillation of the approximate solution, while the other two methods are smoothly varying across the refinement interface. For smooth problems, even though this overshoot creates a larger error for the three-point method, it will eventually be damped out, and can actually causean improvement of the error profile after many time steps. However, in most situations the refined area is usually changed after several time steps, and in such casesone would expect that the three-point method might introduce many oscillations. TABLE I Ax, = Ay, = l/15; Ax, = Ay, = l/90; PI = P, = 8; At = 2 x 10-k Alternating direction Mesh refinement (a) Integer-fit” (b) 4-point (c) 3-point
t, time
IOAr 1OAt
1OAt
EC
Eint
3.7 x 30-2 3.1 x 10-Z 5.5 x 10-z
7.6 x lO-s 7.8 x 10-B 1.2 x 10-Z
a For the entries in (a), ECis a 5 % relative error while Eint is 0.12 %.
In Table I, we list the results of solving (6) by using our method, with the three interface conditions of Section 2. We took the analytical solution as
4~ Y,t>= (t +! a)p . ew i (x -
1/2y
+ (u -
-4(t + a)
w
1’
with j3 = 40 and OL= l/432. To describe the accuracy of the refined scheme away from the interface, we compare the error in the interior of D, at points at a distance of at least, a coarse mesh-width away from the boundary of D, . The maximum deviation of our approximation from the solution over all such points we denote by Eint . The maximum error over the remaining points in D, and over all the coarse points, we denote by E, . Note that, for all our calculations, we always used third-order Lagrange interpolation along the interface. For the time interval indicated in Table I, this function exhibits large gradients near D, . In fact, at t = 0, u(x., y, 0) has a maximum value of 10.8, and U, , U,
524
CIMENT AND SWEET
take on values greater than 80 in absolute value in D, . However, the region of large gradients begins to spread out, and our method remains an effective computational procedure only if we have a dynamically changing grid. In this case, the relative errors reflect the effectivenessof the refinement technique. TABLE II Integer-fit: 10 time steps; At = 0.006 -5 1. No refinements Ax = Ay = l/l5 2. No refinement Ax = Ay = l/60 3. Refinement Ax, = Ay, = l/l5 Pi = P2 = 4
&t
Work estimate
0.004608
0.010235
“(15)?
0.000277
0.000671
h16(15)2
0.00506
0.00201
~2.77(15)~
In Table II, we show the results for a case where the region of sharp gradients is more static. We solved for U(X,y, t) = 1 + exp[--2Ot[(x - i)>” + (y - $))“I] a solution of (6) on D, as in (8), using the integer-# method. With the mesh refinement, we achieve an increase of accuracy by a factor of, at least, five over the interior square (much more away from the interface), by using only 2.77 as many operations. As shown in Table II, the error, with no refinement, is inversely proportional to the number of grid points. Thus, the only other way to have obtained this increased accuracy in the interior region would have been to use a uniform grid
Bx=By=l--.
15 2/5
This would have required nearly twice as much storage and calculation as compared to the mesh refinement scheme. REFERENCES 1. M. CIMENT, Mesh Refinement Techniques-Fundamental Considerations, Internal Report No. 69-5, Marathon Oil Company, Denver Research Center, January, 1969. 2. M. CIMENT, Stable difference schemes with uneven mesh spacings, A&h. Comp. 25 (1971), 219-227. 3. S. OSHER, Mesh refinements for the heat equation, SIAM J. Numer. Anal. 7 (1970), 199-205. 4. D. W. PEACEMAN AND H. H. RACHFORD, JR., The numerical solution of parabolic and elliptic differential equations, SIAM J. 3 (1955), 28-41.
MESH REFINEMENTS
525
5. R. D. RICHTM~ERAND K. W. MORTON, “Difference Methods for Initial-Value Problems,” 2nd Edition, Interscience, New York, 1967. 6. J. VARAH, Stability of difference approximations to the mixed initial boundary value problem for parabolic systems, SIAM J. Namer. Anal. 8 (1971), 598-615. 7. R. S. VARGA, “Matrix Iterative Analysis,” Prentice-Hall, New York, 1962. 8. 3. H. WILRINSON,“The Algebraic Eigenvalue Problem,” Oxford University Press, London/ New York, 1965.