A method of solving stiff boundary-value problems based on operator splitting

139 9. BURKOV V.N., Elements of the Mathematical

Theory of Active Systems,

Nauka, Moscow,

1977.

Translated by Z.L.

U.S.S.R. Comput.gaths.Muth.Phys.,Vol.29,No.6,pp.139

0041-5553/89 $i0.00+0.00

-146,1989

01991

Printed in Great Britain

Pergamon Press plc

A METHOD OF SOLVING STIFF BOUNDARY-VALUE PROBLEMS BASED ON OPERATOR SPLITTING~ A.A. ABRAMOV and L.M. ALVAREZ

An iterative method of solving some stiff linear boundary-value problems for equations of fourth and second order is constructed. The rate of convergence is shown to be not slower than for a geometric progression whose ratio decreases as the stiffness parameter of the original equation increases. Results of some numerical experiments are presented.

Introduction. In this paper, we propose and investigate a method of solving some stiff linear boundary-value problems for ordinary differential equations. The main focus is on the Orr-Sommerfeld fourth-order equation, which plays a central role in hydrodynamic stability (see, e.g., /I/). This equation contains a parameter (Reynolds number), and the original boundary-value problem becomes stiff for large values of this parameter. When this bounary-value problem is solved by some standard version of the privotal-condensation method (see, e.g., /2/), the Cauchy problems for the pivotal-condensation equations turn out to be stiff, which essentially complicates their solution. The paper proposes an iterative method of solving the original problem. Several Cauchy problems are solved in each iteration. Only some of the Cauchy problems for first-order linear equations are stiff, and their solutions can be written in explicit form. The construction of such quadrature formulas that also ensure good accuracy for large p a r a m e t e r values is a much simpler problem than the solution of stiff pivotal-condensation equation. The convergence of the iterative method is accelerated as the parameter value increases. As few as 2-3 iterations are needed for the large parameter values that are relevant in applications. The entire method can be generalized to a broader class of equations. some numerical examples are presented for fourth-order boundary-value problems and for second-order problems with rapidly oscillating solutions.

I. Statement of the problem. we consider the equation

0
y'V--A(x, R)V'+B(x, R)y=l(x),

(t.i)

with the boundary conditions

y (0) =Y,o~

y~(0) =y.,,

R is a large parameter, B~i, x, a , and In what follows, we assume that

R

y (a) =y,.,

are real, all other quantities may be complex.

A(z, a)=RAo(z, R), and A,(z, R), Bo(x, B), yoo, yo,, ylo, y,,, a, ](x) existence of

B(x, R)=RBo(z, ~)

area all bounded;

A|

more precisely,

We assume that

x'-A (z, R)x'+B(z, R)=0 ~=x'.

we assume the

B.(z)=IimB.(x.B),

where A.(x) and B| are smooth functions in [0, a]. take non-positive real values. For each z, construct the characteristic equation

and put

(i.2)

y' (a) =y,,.

Then by our assumptions

the quadratic

~Zh.vych{sZ.Hat.mut.Eiz.,29,12,1800-1810,1989

equation

A.(x) does

not

140

p'-A (z, B) ~+B(z, n) =0 for large H has one large root and one bounded root. respectively. It is easy to verify that

^

B,

-~-=A,-

A,R

We uenote these two roots by

B,"

B,

A,sR 2+ ....

~ and p,

B,"

P='~,+'A-.','R § ....

By the assumption for A., the square roots of ~/R for large R do not lie on the imaginary axis. D e n o t e by ~ the square root of ~/R that lies in the right halfplane. Then there exists a positive number s0 such that for sufficiently large R we have Red--g0. Choose such s0 and fix its value in all that follows. For brevity, we will write ~(x) and ~(z) instead of ~(x,R) and X(x,R).

2. Description of the operator splitting method. We can directly verify that Eq.(l.l) is equivalent to the system

y"-p(z)y=n, n'R-'I'+;L(x) u=v, v'R-V'--~(z) v=g(x),

(2.t) (2.2) (2.3)

g(x) =l(x)R-'- (p"y+2p'y')R-'+I'nR -~.

(2.4)

where

Eq.(l.l) is reduced to the form (2.1)-|2.4) for the following reason. If we take any initial approximations O(x) and ~(x) and substitute them into the right-hand side of (2.4), the next approximations y(x) and u(zJ are obtained by solving Eqs.(2.1)-(2.3). Let us first make more precise the algorithm that produces y(x) and u(z) from 9(z) and ~(x).

Step 1. Step 2.

Use (2.4) to construct

g(x)

from

Obtain a general expression for Cauchy problems.

~(z)

v(x)

v,'a-~-X(z) v,=g(x), v,'R-~--~(z)v,=O,

and

~(z)

from Eq.(2.3).

To this end solve the

v,(a) =0,

v,(a)=l.

Then the general solution of Eq.(2.3) has the form

v(z)=Vo(Z)+C,v,(z) (here and in what follows, Cs and

Step 3.

Cz are arbitrary constants).

From Eq.(2.2) we obtain a general expression for

where the functions

u(z)==,(x)+C,u,(z)+C~u~(z), no(X),u,(x), un(x) are solutions of the

n,'a-"+~.(z)u,=Vo(Z), u,'a-"+~.(z)u,=v,(z), .,'a-"+~. (z) u.=O,

u(x)

in the form

following Cauchy problems:

u,(O)=O, u,(O)=O, u, (0) = 1.

Note that v0, pt, u,, ut and u~ were obtained by solving stable Cauchy problems: and ~t are computed from right to left, and u,, u,, u2 from left to right. The sign of the coefficient l(x) is precisely such that it ensures stability of integration of the equations in the required direction (we recall that Ro~(x)~0>0). Therefore, in particular, for large R all the functions remain bounded uniformly in R.

Step 4.

Finally, we solve Eq.(2.1). Take two of the four boundary conditions in (1.2): y(a) =yt,. These two conditions are the most useful for practical calculations and investigation. Solve the following boundary-value problems:

y(O) =yo,

and

y,"-~(x)y,=u,(x), y,"--9(x)yt=ut(x), y,"--p(x)y,=ui(z), Then the general expression for the solution ditions from (1.2) has the form

y,(O) =you, y,(a)=y,,, y,(O)=y,(a)=O,

y,(O) =y,(a)=O. y(x)

satisfying (2.1)-(2.3) and the two con-

y(z) =y.(~) +c,y, (z)+c,~,(~) Here

(2.5)

~(x) is assumed to be such that the homogeneous problem

y"--p(x)y=O, does not have non-trivial solutions.

y(O)=y(a)=O

If, for instance,

Re~>10,

then this assumption holds

141 and the boundary-value problems given above are numerically stable. The two as yet unknown constants in (2.5) will be determined by using the two remaining conditions from 11.2): y'(0)=y~ and y'(a)=y,,. For C, and C, we obtain a system of linear algebraic equations

c,:,," (o) +c,~,,' (o) =y.-y,' (o),

'(2.6)

c,v,'(a)+c,u,'(,,)=u,,-u.'(,O.

If

I v, (a) y,'Ca) [ #o, we obtain

u(=)

Ci and

C,

by solving this system and thus find the next approximation

y(x) and

Remarks I. Note that in these Cauchy problems, for large R, the functions v0, v,, u0, u,, u2 vary rapidly only near the ends of the interval [0, a], while throughout the rest of the interval their variation is "normal". This enables us to solve the Cauchy problems by methods that produce the prescribed accuracy in a limited number of operations (as a+-). 2. It is easy to note that the functions v,,u,,u,,y,,y, are independent of the initial approximation and do not change during the iterative process; Thus, they may be computed only once. In each iteration, we only recalculate Vo, U0, y0, i.e., solve two Cauchy problems for first-order linear equations and one boundary-value problem for a second-order linear equation. 3. The convergence of the method. The right-hand side of (2.4) contains the functions y and u with the multipliers R -i and R-% respectively, which suggests that the iterative process converges rapidly (not slower than a geometrical progression with a common ratio of order R-~). We will give a rigorous proof of this statement. Denote by ~.(x) and ~.(=) the limiting values of ~(z) and ~(x) as R-+~; X.'(z)=A.(z), BeX.(x)>~0>0, ~.(x)=B.(z)/A.(z). In what follows, we assume that the problem

/'-~.(=)y=o,

y(o)=~(.)=o

(3.t)

does not have non-trivial solutions. The functions v~, u,, uz~ y~, y,, as we have noted above, are computed only once. values determine the matrix

Their

M=II u,'(~) u,:(o)II (a)ll " By assumption, problem 13.1) does not have non-trivial solutions, and therefore for large R the problem

y"-~(z)y=u(z),

y(O)=y(a)=O

(3.2)

has a (unique) solution, which is given by a

(3.3)

y (z) = S K(z, z) u (z) dz, g

where

K(x,z) is Green's function of this problem. We will first establish the principal terms of the representation elements of the matrix M. In order to simplify the subsequent calculations, we introduce

A (z) = i ~"(z) dz,

A~

(0)

(A (a) =0).

z

Then

v,(=)=exp[--a'~A(z)],

u,(z)=exPl--a"[Ao--A(z)]),

u, (x)= R 'sexp[-R'I'A(x) ] i exp{--2RV'[A(z)--A(z) ]}dz. 0

Together with 13.3), we have the equality

y' (x) = I Kz' (x, z) u (z) d.. 0

(for large

R) of the

142 Hence

y'(O)= K,,'(O,z)a(z)dz,

y'(a)=~ K,/(a,z)u(z)dz.

o

o

The function K='(0,z) is continuous in z in (0, a); its limiting value is -I as z--0 and 0 as z~a. The function K.'(a,z) is also continuous in (0, a). Its limiting value is 0 as z ~ 0 and 1 as z~a. Hence in the usual way we obtain a

y,'(o) = i ~p{-m'[Ao-A(~)]}d==- ~ ~p[-m'z~ (0)]d~, o

o

because for large R only a small part of the integration interval near the point z=0 is relevant; outside this part, ReB~[A0 -A(z)] is noticeably greater and exp{-R~[A0-A(z)]} is exceedingly small. Finally we obtain

y,' (0) = -- [ a ~ (0) ] - ' + o (a -~). We similarly have 'I y, ' (a)=o(n-,),

', v, ' (o)=o(n-,),

y,' (a) = [2R't'l ' (a) ] - ' + o (R-V'). Therefore

l-,+o(R-',.)

]]

and

~l_i=m, [[

o(0

2~'(=)+o(0

o(,)

-~.(o)+o(l)

(3.4)

II

We can now proceed with a direct proof of convergence and derive a bound on the rate of convergence. Since the iterative process is linear, it suffices to consider the case ](x)=O, yoo=yo,=y,~ and to establish the convergence (and rate of convergence) to zero of the functions vs u, y. Suppose that for the initial v a l u e s ~(x) and ~(x) we have

I~(x)l
I~(x)l
Ir

where h(x)=exp(--aoR~x)+exp[--aoR~(a-x)] and p is a number. We will show that y(z) and a(x) produced by one iteration satify bounds of the same type: only instead of p these bounds contain the multiplier pmR -%, where m is some bounded constant (as R - ~ ). This result obviously shows that the required convergence holds for large R. Let us derive the successive bounds. Under our assumptions, we have (in what follows, mr, m=,.., are some constants, bounded as B ~ ) '~ -'b Igl
Since

]K(x,z)l~m,

and

IK='(x, z)l~
pro, + pro, [ h(x)dz < pm. I~o1< n m' ; B and correspondingly

[Yo'[<~pm,R-";. C=- the solutions of system {2.6) ]C,l~pmi,R-~", IC,l~pm,,tl-" and t h e r e f o r e lyl
Using formula (3.4), we obtain for C, and

interval Y'(0)=y0,

[--a,a]. Then, in order to pass to the interval for x=O are replaced by the conditions y' (0) =y'" C0) =0.

[01 a] , the conditions

in the y(0)=y00, (4.1)

143 In what follows, A(x), B(x),](x) are even functions; then ~(x) and ~(z) are even, and we assume that ;t'(0)=0. These changes in the problem do not lead to any significant changes in the algorithm and in the analysis of its convergence. From Eq.(2.1) it follows that y'"--~'y--~y'=u'. By (4.1) and ~'(0)=0, this gives u'(0)=0. Then we obtain the condition

u(O)=v(O)ll(O).

(4.2)

The final algorithm in this case is the following: from before, v0, pl, u0, ut, u2; solve the problems

~(z)

y~176176 y,"-~(z)y,ffiu,(z),

y,'(o)=o, y~176 y,'(O)=y,(a)=O,

y,"-~(x)y,=u,(x),

y~'(o)=y,(a)ffiO

and

~(x)

find, as

and obtain

y(z)=yo(z)+C,y,(z)+C,y,(z), u(x)=Uo(X)+C,u,(x)+C,u,(x); these Ci and Ca are sought using the remaining two conditions: y'"(0)=0 (which is equivalent to (4.2)) and y'(a)=y,,. The convergence of the iterative process in this case is investigated along the same lines, introducing Green's function K(z, z) of the problem

y"-~(x)y=u(x),

y'(O)=y(a)=O.

4.8. 2he second-order equation. An important class of problems are boundary-value problems associated with a stiff second-order equation. We will describe the application of the operator splitting method to the solution of the boundary-value problem

y"-p' (x) y=l(z), O<.z<.a, y(O)=y,, y(a)=y,.

(4.3) (4.4)

If p'(x)=Rpo*(x,B), B ~ l , and po(x,l~) is a bounded function (of order of unity), then the boundary-value problem (4.3), (4.4) is stiff. Eq.(4.3) is equivalent to the system

y'B-~+po(x)y=z, z'a-~-po(x)z=l(z)R-'+po'(X)a-~y.

(4.5a) (4.5b)

The algorithm now takes the following form.

Step I. Rep0(x)~a0>0

Using ~(x) (the initial approximation) solve two Cauchy problems, which for have a stable numerical solution from right to left:

Zo'R-~-po(Z)Zo=]Cx)R-'+po'(X)R-~9(=), z/R-~-p~(x)z,=O, z,(a)=t. Then we obtain a general solution of (4.5b) with

y=~

zo(a)=O,

on the right-hand side:

z Cz) =Zo(X)+C,z,(x).

Step 2.

The last step of the algorithm solves three Cauchy problems from left to right, which also have numerically stable solutions:

yo'R-",+po(X)y~176

y~ =o,

yi'B-"+p~ (x) y,=z, (z),

y, (0) =0,

y,'R-'~'+po(x)y,=O,

y2(O)=l.

We obtain y(x)=yo(x)+C, yi(x)+C2y,(z). The two unknown constants Ct and C, are obtained from conditions (4.4) . If Rep0(x, R ) ~ a a > 0 for all z, then using the same technique as for the fourth-order equation, we prove that the iterative method converges not slower than a geometrical progression with the ratio NR-V'. Remarks. 3. Note that five Cauchy problems are solved in the first step of the iterative process and only two (for z0 and y0) in all other steps. 4. If the equation has the form

f'+f(z) y=/Cx), where

p(z) is real, it can be rewritten as

~'- tip(z)]'y=1(z).

144 The algorithm formally remains the same as before, but all calculations are carried out in the complex domain. The investigation of the rate of convergence conducted previously for R0;0(x,B)ma~>0 is obviously inapplicable in this case. 5 . If we consider a more general boundary-value problem

~"+~ (=)u'+q(=) y=l(=),

0
then the operator splitting method remains valid under certain assumptions; on the left-hand sides of the corresponding equalities now have the form y'-p,(z)y and d-p,(z):, where ~, and p, are the roots of the characteristic equation ~'+T(z)~+q(=)=0.

5. gumers

4ap~ementation of the method.

In this section, we consider numerical methods of solving the auxiliary initialboundary-value problems that arise in the process.

and

5.1. A numerical method fop soZuing the auziZ4ary Cauchy probZems, since all the auxiliary Cauchy problems involve first-order linear equations , they can be solved explicitly in quadratures. We will only consider the formulas for u(z); the problem for r(x) is solved similarly. We have

u'R-'%X(z)u=v(z), where

Re~(x)>a,>O

.(0)=0,

(s.1)

and N>t.

Define a grid O=x,<...
](x), we are

z

o

The computations are carried out from left to right. W e will give a recursive formula that expresses u,+, in terms of uA and the given quantities. From (5.1) for

k=0,| .....N--i

we have

u~+,=u, exp [ - R ~ (A~+,-A~)] +~,,

(5.2)

where fl.t

@,,=

lexp{_R,t,[A,H_A(t)] } ~(t) v!t! d{-lt'[A,+,-A(t)]}. =l

Reduce

~, to the form r..t

(~, Imm I e][p{R'h[Al+i--A(') ]1 zl

Put ~(t)=v(t)/l(t). Let

v-!t! dt__R,~(A,+,__A(t))].

~(t)

v(t) be a "good" function.

Represent

A~.,-A(t) [ (~,-~,+,)+O((t,,-t,)')

Then for the principal part of obtain

r ~

~(t) in the form

].

~(t) integration can be carried out explicitly, and

(5.3)

we

n+i {i - exp[-nh(A~,,-~) ]}+

vJl,- v,+,l~,+,{I- [i + R '~(A~+,-A,) ]exp [ - R '~(A,+t--A,)]}. a 'h(,%+,--A,) From (5.2) we now obtain the required recursive formula P& § i 'b u,+,~u, exp[-R ,& (/%+,--A,) ]+-7----{t-exp[--R (A,+,-A~) ])+

t,~+l

(5A)

145

v~/l~--v,+,ll~+, { t - [ t+Rv' (Ah+,-A~)]exp[--lr K" (&§

1},

which has the following properties: a) it is stable for any z-stepping and any R (provided, of course, b) it gives a global error O(h.mln(h,R-'h)) on [0, a] , where

ne ~.(z) ;~ ~0>0) ;

h=max (~+,-~). k This result was obtained by integrating the error in (5.3) over the interval [~,zA+,] and then summing over all intervals. For ~(t) we used an interpolation formula linear in A i 0 . If (5.3) is replaced with an expression which is a polynomial of higher degree in A(t), then we obtain an explicit formula with property a) that has an accuracy of higher order (an analogue of Adams formulas, see /2/).

5.2. Numer~caZ soZution of boundary-ua~ue problems. The auxiliary boundary-value problems (3.2) under our assumptions are numerically stable and can be solved easily by any version of the pivotal condensation method. The numerical experiments were carried out by the method proposed in /3/. The Cauohy problems for the pivotal condensation equations were solved by the RungeKutta-Felberg method using the subroutine RKF45 in /4/. For a prescribed accuracy, this subroutine automatically chooses the step and solves non-stiff and "almost non-stiff" Cauchy problems. This is precisely what we need, because the right-hand sides (uo,u,,u,) of the boundary-value problems have a boundary layer. Since we only have the numerical values of u0, u,, u, on a grid, a difficulty arises: these functions must be interpolated when the program changes the step. Cubic spline interpolation was used (the subroutines SPCOEF and SPLINE, see /5/). $.$. Selecting a non-un~foz~ gP~d. For the standard methods of solving the Cauchy problem that do not allow for the stiffness of the differential equations we cannot use large steps inside the interval (the solution "destabilizes"). But to integrate functions with a boundary layer we can use a large step in the middle part of the interval [0, a], and a fine grid is needed only near the ends. We select a symmetric grid that depends on R; the higher R, the greater is the condensation of the gird at the ends. After some experimentation, we obtain the following useful final formula for grid construction. Specify N - the number of intervals into which the interval .[0,a] is partitioned by the nodes 0=x0<...
,v A ,, B Y --'-~Y + " ~ v = O ,

y(t)=t,

t<~z~2,

#'(1)=4,

#(2)=16,

#'(2)=32,

where A = R , B=I2R--24, Rml. The corresponding solution is mation in the iterative process, we used #(z)=1, u(z)=l.

0

i ~l

1.3

1.999

Z z

0 ~ i ' I

I.!

Fig 91

1.3

#(z)=z ~. As the initial approxi-

1.g99

:r

Fig 92

"z(z) %(=)

O!

t

/.0~ /.7 1.9 Z z

Fig.3 USSR 29:6-J

01

t.OOi I 1.7

Fig.4

1.9 Z x

146 Figs.l-4 plot the results of calculations of the last approximations for R=I0' and N=51. The iterations were continued until the deviation of the current approximation from the previous approximation became uniformly less than 8 (in these examples, ~=|0-~). The deviation of the approximate solution from the exact solution was less than 10 -~. The number of iterations decreased significantly as R increased: for ~=I0',10~,|0 ' and 10' the number of iterations required was respectively 6, 4, 3, and 2. 2. Consider test examples for the second-order equations from /6/:

~"+~=i0'cos(80t), 0~t~, a. For

y(0)=1, ~(I)=0

~=10 t the solution is

#(t)=C, cos (tOOt)+C, sin (tOOt)+#(t), where

Ci and C2 a r e computed from t h e boundary c o n d i t i o n s and #=|O=(~-6400)-tXeos(80t). The calculations were carried out by the method described in Sect.4.2. We used a uniform grid with N = 2 0 0 and ~=I0 -~. b. For 7=--|0 ~ the solution has the form

y ( t ) =C,e'*~

Cze-tO~ f ( t ) ,

where

~(t) is calculated using the same formula as in paragraph a above. A uniform grid was used with N = 1 0 0 and 8=10 -~. Since the coefficients of the equation in cases a and b are constant, no iterations were required. The results reported in /6/ (they are identical with our results) were obtained by a fairly complicated method on a CDC-CYBERI76 computer. We used an IBM/XT personal computer and programs written in FORTRAN-77 (in double precision). All the results were generated in a few seconds. REFERENCES I. GOL'DSHTIK M.A. and SHTERN V.N., Hydrodynamic Stability and Turbulence, Nauka, Novosibirsk, 1977. 2. BAKHVALOV N.S., Numerical Methods, Nauka, Moscow, 1973. 3. ABRAMOV A.A., On transport of boundary conditions for systems of linear ordinary differential equations (a version of the pivotal condensation method), Zh. Vychisl. Matem~ i Mat. Fiz., i, 542-545, 1961. 4. SHAMPINE D.F. and WATTS H.A., Practical solution of ordinary differential equations by Runge-Kutta methods , SAND76-0586, 1976. 5. SHAMPINE D.F. and ALLEN R.C., Numerical Computing: An Introduction, Saunders, Philadelphia. 1973. 6. MEYER G.H., Continuous orthonormalization for boundary-value problems, J. Comput. Phys., 62, 248-262, 1986.

Translated by Z.L.