- Email: [email protected]
Efficient Two-Step Numerical Methods for Parabolic Differential Equations
EFFICIENT TWO-STEP NUMERICAL METHODS FOR PARABOLIC DIFFERENTIAL EQUATIONS Arieh Iserles Department of Applied Mathematics and Theoretical Physics University of Cambridge Cambridge CB3 9EW, England
We develop two-step A-stable methods of maximal order for the numerical solution of ordinary differential systems. When these methods are applied to the stiff, large systems which originate from parabolic differential equations they yield a set of algebraic equations of special form. This set is easier to solve than the algebraic equations which are obtained when using one-step methods of the same order. When the spatial variables of a parabolic differential equation are discretized either by finite differences or by finite elements, one obtains a stiff ordinary differential system of the form
u'
(1)
=
+ f(t,u), u(to)
Mu
=
uo
,
where M is a large, sparse, constant, stable matrix, p ( M ) > > O , and the function f has a Lipschitz constant much smaller than the spectral radius of M.
To illustrate the approach of our paper, let us compare two onestep methods for the system (1): (i) diagonal Obrechkoff method [7]: U1
where
u1
=
'
- 'jhul
u(to
+
+
1 2 " 12 h u 1 = u 0
-
+
%hue'
+
1 2 " 12 h u0
h);
(ii) Ndrsett-Makinson restricted Obrechkoff method [ 4 ] ,
.
[6]:
Both methods employ the first two derivatives of u They are A-stable. The first is of order 4 , while the second is of order 3. While solving (2), it is u s u a l to evaluate an algebraic linear system with a matrix 1 2 2 I - S h M + - h M . 12 The term M 2 causes large fill-in, which makes the application of either iteratives or direct methods for the solution of linear algebraic systems expensive. However, method (3) was proposed because the matrix of its linear system has the form 319
320
(I
A . ISERLES
-
'.
(Ji+fi/G)hM)
T h e r e f o r e t h e s y s t e m c a n b e s o l v e d i n two s t a g e s ,
where t h e m a t r i x o f e a c h s t a g e i s I advantage of avoiding f i l l - i n .
- (f +
J 3 / 6 ) h l l , which h a s t h e
Hence, t h e method ( 2 ) ( a n d t h e d i a g o n a l O b r e c h k o f f m e t h o d s i n general [81) is superior a s f a r a s t h e order i s concerned, while t h e method ( 3 ) ( a n d t h e r e s t r i c t e d O b r e c h k o f f m e t h o d s i n g e n e r a l ) y i e l d s a more c o n v e n i e n t a l g e b r a i c s y s t e m . W e p r o p o s e m e t h o d s which s h a r e t h e m e r i t s o f b o t h d i a g o n a l and r e s t r i c t e d O b r e c h k o f f m e t h o d s . The p r i c e w e pay i s moving from o n e s t e p t o t w o - s t e p m e t h o d s . L e t u - ~= u ( t o - h ) . An e x a m p l e o f a method w h i c h c o r r e s p o n d s t o ( 2 ) and ( 3 ) i s
u1
-
(I+J?/3)hu;
+
( 1 / 3 + f i / 6 ) h 7 u i = -2fi/3huA
+
(l+6/3hull
+
+
(1/3+fi/6)h2u11
u - ~+
.
I t i s A - s t a b l e and o f o r d e r 4 and y i e l d s a n a l g e b r a i c s y s t e m w i t h
t h e matrix (I
- (f +
0/6)hM)'.
I n g e n e r a l , w e c o n s i d e r m e t h o d s o f o r d e r 2n o f t h e form
where a > 0 , i n o r d e r t h a t t h e m a t r i x o f t h e l i n e a r s y s t e m h a s t h e form ( I
-
I t i s c o n v e n i e n t t o a n a l y s e m e t h o d s ( 4 ) i n terms o f a p p r o x i m a t i o n s t o t h e exponential function. I n d e e d , f o r ut;C" ( 4 ) i s e q u i v a l e n t , i n operational notation, t o
n hD {(l-ahD) e where D =
(5)
d
.
ez(l-az)"
-
n
1 b k ( h D )k - n1 c k ( h D ) k e - h D } u ( t o ) = 0 k=o k=o
,
C o n s e q u e n t l y ( 4 ) i s o f o r d e r 2n i f and o n l y i f n
- 1 bkzk k=o
e-'
n
1 ckzk
k=o
= O(z
2n+l
)
.
Equation (5) g i v e s a u s e f u l t o o l f o r t h e d e r i v a t i o n of t h e e x p l i c i t form o f t h e c o e f f i c i e n t s bk and c k . Theorem 1:
-
L e t co = 0 ,
-
bo = 1
and
321
E F F I C I E N T TWO-STEP N U M E R I C A L METHODS
k = l ,
. .,
k = l , . Then the method (4) with bk = bk, ck order 2n.
=
.,
n;
n
.
ck, k = 0,
...,
n,
is of
In the proof of the theorem we show that the given coefficients satisfy (5). This involves long and tedious algebra and numerous applications of the theory of hypergeometric functions. For any fixed N (5) reduces to 2n + 1 linear equations in 2n + 2 unknowns. It is an easy exercise in linear algebra to show that the matrix of this system is of rank exactly 2n + 1. Therefore we have an extra degree of freedom, in addition to a . Theorem 2: The general form of the coefficients of the method (4) of order 2n is
(2n-k)! C k = C k+ - (2n)!
(:)a
, k
=
0,
...,
n;
5
where bk and ck
are given by (6) and
is an arbitrary constant.
The proof of this theorem rests upon the fact that we can perturb the second and third terms in (5) by -BPn(z) and BQn(z) respectively, where Pn/Qn is the n-th diagonal Pad6 approximation to exp(z) and is an arbitrary constant, without changing the order. Using the theory of weakly nonlinear methods [5] we are able to show that method (4) is zero-stable [3] if and only if -1 < B 51. It is evident that the classical definition of A-stability [ 3 ] is equivalent to /R+(z)I, IR-(z)I < 1 for every z, Rez < 0, where R+ and R- are the solutions of the quadratic (7)
A R ~- BR A(z) = (
method (4) is
- c = o ,
322
A. ISERLES
Theorem 3: The method (4) is A-stable if and only if the three following conditions are satisfied: (i)
(ii)
> 0,
-1 < p c
[A(it)1 ’
2 IC(it)
ci
1 ; 7
.
{/A(it)I2-1C(it)I for every real t ;
It is L-stable if and only if in addition
deg8, degC
5 n-1 .
As usual with A-stability analysis, the proof is based upon the application of the maximal modulus principle. Hence, for ci > 0 it is sufficient to investigate the solutions of (7) along the imaginary axis. It is easy to see that the A-stability is equivalent to the satisfaction of the root condition by the quadratic (7) for z = it, --m < t < m Using a suitable modification of Schur-Cohn criterion [ 2 ] , we are able to establish (ii) and (iii). The condition for L-stability follows from an investigation of the asymptotic behaviour of R+ and R- when I z l + m .
.
Figure 1 gives the A-stability regions in the ( a , B ) plane for 151153. Note that the regions are closed, with the exception of the portion of the boundary where
a=
-1.
There are no L-stable methods for 1 & n 3, with the single exception of the second order formula of Gear 181, namely n = 1, a = -23 r P = - 17 . Furthermore, computer search showed that there are no A-stable methods (4) for n = 4. By noticing this and the steady contraction of the stability regions in Figure 1 for increasing n, we conjecture that there are no A-stable methods (4) for n 2 4
.
Finally we examine the effect of solving the algebraic system which is obtained by various methods from the diffusion equation
where u = u(t;xl, ..., xm) and A is the Laplace operator. We assume that the equation is solved in an equidistant square spatial grid and the second derivatives in respect to xi, i = 1, ..., m, are discretized by the usual three-point formulae. The comparison is between the usual methods (for example the diagonal Obrechkoff), which produce a system with a matrix which is a polynomial in M,
323
E F F I C I E N T TWO-STEP N U M E R I C A L METHODS
n=2
n=1 1
Figure 1
324
A . ISERLES
and the restricted methods (like the restricted Obrechkoff and method ( 4 ) ) , which produce a system with a factorizable matrix (I-ahM)n, where h is the time step and a > 0. If a direct method, like Hockney or CORF scheme [l], is used to solve a system with the matrix I-ahM, the computation consists mainly of a large number of solutions of (2m-l)-diagonal q x q systems, where q is the number of internal points in the spatial grid in any direction. However, if such a method is used to solve a system with the matrix P(hM), where P is a polynomial of degree n, one must solve instead (2n(m-l)+l)-diagonal q x q systems. Hence, if NR and No are the numbers of multiplications involved in the solution of systems with the matrices (I-ahM)" and P(hM) respectively, then
-
NR
=
C
+
No
=
C
+ (n(m-1)(2n(m-l)+1)+1)D
n(2m'
3m
+
2)D
;
where C and D are constants of equal order of magnitude. Table 1 displays the number of multiplications which are required The advantage of to solve the problem for m = 2 and m = 3 method ( 4 ) is obvious.
.
Now we examine the efficiency of iterative schemes [9] for sparse algebraic systems, in conjunction with the restricted and the usual methods. We restrict ourselves to two space variables only and to n = 2 . Let g be the size of the space grid and u = h/g2 . The spectral radius of the Jacobi method for the matrix I - clhM is
while the spectral radius of the SOR method for this matrix is 4aucos 2q+l
However, it is possible to show that the spectral radius of the Jacobi method for the matrix I - fhM + x1 h 2 M 2(which originates from the second diagonal Obrechkoff) satisfies
Computer calculations estimated the spectral radius of the SOR method for this matrix as l-2/max{qIu} + O(l/(max{q,a~)' ) . Hence,
325
E F F I C I E N T TWO-STEP NUMERICAL METHODS
Two space variables: diagonal Obrechkoff
+
restricted Obrechkoff
method (4)
order 4
C
11D
C
+ 12D
C
+
8D
order 6
C t 37D
C
+ 20D
C
+
12D
Three space variables: diagonal Obrechkoff
+
order 4
C
order 6
C t 127D
37D
restricted Obrechkoff
method (4)
C + 33D
C
C t 55D
C
+ +
22D 33D
Table I : The number of multiplications required by the various methods.
the SOR method for I - ahM converges asymptotically faster by a factor of %Jmax{q,o} than the corresponding method for 1 I - 4hM + l;?h2M2. References Golub, G., Direct methods for solving elliptic difference equations, in: Morris J.L. (ed.), Symposium on the Theory of Numerical Analysis (Lecture Notes in Math. 193, Springer Verlag, Berlin, 1971). Henrici, P., Applied and Computational Complex Analysis (Wiley, New York, 1976). Lambert, J.D., Computational Methods in Ordinary Differential Equations (Wiley, London, 1973). Makinson, G.J., Stable high order implicit methods for the numerical solution of systems of differential equations, The Computer Journal 3 (1968) 305-310. Makela, M., Nevanlinna, 0. and Sipila, A.H., On the concepts of convergence, consistency and stability in connection with some numerical methods, Numer. Math. 22 (1974) 261-274. N$rsett, S . P . , One-step methods of Hermite type for numerical integration of stiff systems, BIT 14 (1974) 63-17.
326
A . ISERLES
[7] Obrechkoff, N., Sur les quadraturcs mecaniques (Bulgarian, French summary), Spisanie Bulgar. Akad. Nauk 65 (1942) 191-289. [81
Wanner, G . , Hairer, E. and Ndrsett, S.P., Order stars and stability theorems, BIT 18 (1978) 475-489.
191
Young, D . M . , Iterative Solution of Large Linear Systems (Academic Press, New York, 1971).
We develop two-step A-stable methods of maximal order for the numerical solution of ordinary differential systems. When these methods are applied to the stiff, large systems which originate from parabolic differential equations they yield a set of algebraic equations of special form. This set is easier to solve than the algebraic equations which are obtained when using one-step methods of the same order. When the spatial variables of a parabolic differential equation are discretized either by finite differences or by finite elements, one obtains a stiff ordinary differential system of the form
u'
(1)
=
+ f(t,u), u(to)
Mu
=
uo
,
where M is a large, sparse, constant, stable matrix, p ( M ) > > O , and the function f has a Lipschitz constant much smaller than the spectral radius of M.
To illustrate the approach of our paper, let us compare two onestep methods for the system (1): (i) diagonal Obrechkoff method [7]: U1
where
u1
=
'
- 'jhul
u(to
+
+
1 2 " 12 h u 1 = u 0
-
+
%hue'
+
1 2 " 12 h u0
h);
(ii) Ndrsett-Makinson restricted Obrechkoff method [ 4 ] ,
.
[6]:
Both methods employ the first two derivatives of u They are A-stable. The first is of order 4 , while the second is of order 3. While solving (2), it is u s u a l to evaluate an algebraic linear system with a matrix 1 2 2 I - S h M + - h M . 12 The term M 2 causes large fill-in, which makes the application of either iteratives or direct methods for the solution of linear algebraic systems expensive. However, method (3) was proposed because the matrix of its linear system has the form 319
320
(I
A . ISERLES
-
'.
(Ji+fi/G)hM)
T h e r e f o r e t h e s y s t e m c a n b e s o l v e d i n two s t a g e s ,
where t h e m a t r i x o f e a c h s t a g e i s I advantage of avoiding f i l l - i n .
- (f +
J 3 / 6 ) h l l , which h a s t h e
Hence, t h e method ( 2 ) ( a n d t h e d i a g o n a l O b r e c h k o f f m e t h o d s i n general [81) is superior a s f a r a s t h e order i s concerned, while t h e method ( 3 ) ( a n d t h e r e s t r i c t e d O b r e c h k o f f m e t h o d s i n g e n e r a l ) y i e l d s a more c o n v e n i e n t a l g e b r a i c s y s t e m . W e p r o p o s e m e t h o d s which s h a r e t h e m e r i t s o f b o t h d i a g o n a l and r e s t r i c t e d O b r e c h k o f f m e t h o d s . The p r i c e w e pay i s moving from o n e s t e p t o t w o - s t e p m e t h o d s . L e t u - ~= u ( t o - h ) . An e x a m p l e o f a method w h i c h c o r r e s p o n d s t o ( 2 ) and ( 3 ) i s
u1
-
(I+J?/3)hu;
+
( 1 / 3 + f i / 6 ) h 7 u i = -2fi/3huA
+
(l+6/3hull
+
+
(1/3+fi/6)h2u11
u - ~+
.
I t i s A - s t a b l e and o f o r d e r 4 and y i e l d s a n a l g e b r a i c s y s t e m w i t h
t h e matrix (I
- (f +
0/6)hM)'.
I n g e n e r a l , w e c o n s i d e r m e t h o d s o f o r d e r 2n o f t h e form
where a > 0 , i n o r d e r t h a t t h e m a t r i x o f t h e l i n e a r s y s t e m h a s t h e form ( I
-
I t i s c o n v e n i e n t t o a n a l y s e m e t h o d s ( 4 ) i n terms o f a p p r o x i m a t i o n s t o t h e exponential function. I n d e e d , f o r ut;C" ( 4 ) i s e q u i v a l e n t , i n operational notation, t o
n hD {(l-ahD) e where D =
(5)
d
.
ez(l-az)"
-
n
1 b k ( h D )k - n1 c k ( h D ) k e - h D } u ( t o ) = 0 k=o k=o
,
C o n s e q u e n t l y ( 4 ) i s o f o r d e r 2n i f and o n l y i f n
- 1 bkzk k=o
e-'
n
1 ckzk
k=o
= O(z
2n+l
)
.
Equation (5) g i v e s a u s e f u l t o o l f o r t h e d e r i v a t i o n of t h e e x p l i c i t form o f t h e c o e f f i c i e n t s bk and c k . Theorem 1:
-
L e t co = 0 ,
-
bo = 1
and
321
E F F I C I E N T TWO-STEP N U M E R I C A L METHODS
k = l ,
. .,
k = l , . Then the method (4) with bk = bk, ck order 2n.
=
.,
n;
n
.
ck, k = 0,
...,
n,
is of
In the proof of the theorem we show that the given coefficients satisfy (5). This involves long and tedious algebra and numerous applications of the theory of hypergeometric functions. For any fixed N (5) reduces to 2n + 1 linear equations in 2n + 2 unknowns. It is an easy exercise in linear algebra to show that the matrix of this system is of rank exactly 2n + 1. Therefore we have an extra degree of freedom, in addition to a . Theorem 2: The general form of the coefficients of the method (4) of order 2n is
(2n-k)! C k = C k+ - (2n)!
(:)a
, k
=
0,
...,
n;
5
where bk and ck
are given by (6) and
is an arbitrary constant.
The proof of this theorem rests upon the fact that we can perturb the second and third terms in (5) by -BPn(z) and BQn(z) respectively, where Pn/Qn is the n-th diagonal Pad6 approximation to exp(z) and is an arbitrary constant, without changing the order. Using the theory of weakly nonlinear methods [5] we are able to show that method (4) is zero-stable [3] if and only if -1 < B 51. It is evident that the classical definition of A-stability [ 3 ] is equivalent to /R+(z)I, IR-(z)I < 1 for every z, Rez < 0, where R+ and R- are the solutions of the quadratic (7)
A R ~- BR A(z) = (
method (4) is
- c = o ,
322
A. ISERLES
Theorem 3: The method (4) is A-stable if and only if the three following conditions are satisfied: (i)
(ii)
> 0,
-1 < p c
[A(it)1 ’
2 IC(it)
ci
1 ; 7
.
{/A(it)I2-1C(it)I for every real t ;
It is L-stable if and only if in addition
deg8, degC
5 n-1 .
As usual with A-stability analysis, the proof is based upon the application of the maximal modulus principle. Hence, for ci > 0 it is sufficient to investigate the solutions of (7) along the imaginary axis. It is easy to see that the A-stability is equivalent to the satisfaction of the root condition by the quadratic (7) for z = it, --m < t < m Using a suitable modification of Schur-Cohn criterion [ 2 ] , we are able to establish (ii) and (iii). The condition for L-stability follows from an investigation of the asymptotic behaviour of R+ and R- when I z l + m .
.
Figure 1 gives the A-stability regions in the ( a , B ) plane for 151153. Note that the regions are closed, with the exception of the portion of the boundary where
a=
-1.
There are no L-stable methods for 1 & n 3, with the single exception of the second order formula of Gear 181, namely n = 1, a = -23 r P = - 17 . Furthermore, computer search showed that there are no A-stable methods (4) for n = 4. By noticing this and the steady contraction of the stability regions in Figure 1 for increasing n, we conjecture that there are no A-stable methods (4) for n 2 4
.
Finally we examine the effect of solving the algebraic system which is obtained by various methods from the diffusion equation
where u = u(t;xl, ..., xm) and A is the Laplace operator. We assume that the equation is solved in an equidistant square spatial grid and the second derivatives in respect to xi, i = 1, ..., m, are discretized by the usual three-point formulae. The comparison is between the usual methods (for example the diagonal Obrechkoff), which produce a system with a matrix which is a polynomial in M,
323
E F F I C I E N T TWO-STEP N U M E R I C A L METHODS
n=2
n=1 1
Figure 1
324
A . ISERLES
and the restricted methods (like the restricted Obrechkoff and method ( 4 ) ) , which produce a system with a factorizable matrix (I-ahM)n, where h is the time step and a > 0. If a direct method, like Hockney or CORF scheme [l], is used to solve a system with the matrix I-ahM, the computation consists mainly of a large number of solutions of (2m-l)-diagonal q x q systems, where q is the number of internal points in the spatial grid in any direction. However, if such a method is used to solve a system with the matrix P(hM), where P is a polynomial of degree n, one must solve instead (2n(m-l)+l)-diagonal q x q systems. Hence, if NR and No are the numbers of multiplications involved in the solution of systems with the matrices (I-ahM)" and P(hM) respectively, then
-
NR
=
C
+
No
=
C
+ (n(m-1)(2n(m-l)+1)+1)D
n(2m'
3m
+
2)D
;
where C and D are constants of equal order of magnitude. Table 1 displays the number of multiplications which are required The advantage of to solve the problem for m = 2 and m = 3 method ( 4 ) is obvious.
.
Now we examine the efficiency of iterative schemes [9] for sparse algebraic systems, in conjunction with the restricted and the usual methods. We restrict ourselves to two space variables only and to n = 2 . Let g be the size of the space grid and u = h/g2 . The spectral radius of the Jacobi method for the matrix I - clhM is
while the spectral radius of the SOR method for this matrix is 4aucos 2q+l
However, it is possible to show that the spectral radius of the Jacobi method for the matrix I - fhM + x1 h 2 M 2(which originates from the second diagonal Obrechkoff) satisfies
Computer calculations estimated the spectral radius of the SOR method for this matrix as l-2/max{qIu} + O(l/(max{q,a~)' ) . Hence,
325
E F F I C I E N T TWO-STEP NUMERICAL METHODS
Two space variables: diagonal Obrechkoff
+
restricted Obrechkoff
method (4)
order 4
C
11D
C
+ 12D
C
+
8D
order 6
C t 37D
C
+ 20D
C
+
12D
Three space variables: diagonal Obrechkoff
+
order 4
C
order 6
C t 127D
37D
restricted Obrechkoff
method (4)
C + 33D
C
C t 55D
C
+ +
22D 33D
Table I : The number of multiplications required by the various methods.
the SOR method for I - ahM converges asymptotically faster by a factor of %Jmax{q,o} than the corresponding method for 1 I - 4hM + l;?h2M2. References Golub, G., Direct methods for solving elliptic difference equations, in: Morris J.L. (ed.), Symposium on the Theory of Numerical Analysis (Lecture Notes in Math. 193, Springer Verlag, Berlin, 1971). Henrici, P., Applied and Computational Complex Analysis (Wiley, New York, 1976). Lambert, J.D., Computational Methods in Ordinary Differential Equations (Wiley, London, 1973). Makinson, G.J., Stable high order implicit methods for the numerical solution of systems of differential equations, The Computer Journal 3 (1968) 305-310. Makela, M., Nevanlinna, 0. and Sipila, A.H., On the concepts of convergence, consistency and stability in connection with some numerical methods, Numer. Math. 22 (1974) 261-274. N$rsett, S . P . , One-step methods of Hermite type for numerical integration of stiff systems, BIT 14 (1974) 63-17.
326
A . ISERLES
[7] Obrechkoff, N., Sur les quadraturcs mecaniques (Bulgarian, French summary), Spisanie Bulgar. Akad. Nauk 65 (1942) 191-289. [81
Wanner, G . , Hairer, E. and Ndrsett, S.P., Order stars and stability theorems, BIT 18 (1978) 475-489.
191
Young, D . M . , Iterative Solution of Large Linear Systems (Academic Press, New York, 1971).