- Email: [email protected]
The construction of stable difference schemes for second-order alternating linear differential operators
86
Yu. A. Eremin and A. S. Leotwv
against the integrationstep for differentvaluesof the parameterS2= w L/c. The continuouscurves givethe errorobtainedwhen the scheme (2.4), (2.5) (u = %) is used, and the broken curves the error when the standard Runge-Kutta method with a constant step is used. It can be seen that the scheme (2.4), (2.5) gives high accuracy with a step over 10 times the maximum permissible step for the Runge-Kutta method. The authors thank E. F. Sabaev for su~esting the problem and for valuable co~ents. 2hnslated by D. E. Brown REFERENCES 1.
KATSKOVA, 0. N., and KRAIKO, A, N., Computation of plane and axisymmetric supersonic flows in the presence ofirreversible processes, PrtkLMot. tekhn, Fiz., No. 4,116-l 18,1963.
2.
KATSKOVA, 0. N., and KRAIKO, A. N., Computation of plune and axisymmetric supersonic flows in the presence of irreversibleprocesses, VTs Akad, Nauk SSSR, Moscow, 1964.
3.
GALYUN, N. S., and KRAIKO, A. N., Computation Mekhan ~Mashinos~, No. 6,41-47,1964.
4.
KAMZOLOV, V. N., and PIRUMOV, U. G., Computation of non-equIIibrIum flows in nozzles, Zzv. Akad Nauk SSSR, Mekhan Zhidkosti i Gaza, No. 6,25-33,1966.
5.
of non-equilibrium
flows, Izv. Akud Nauk SSSR,
GROMOV, V. G., Chemically non~qu~b~um boundary layer in dissociated air, Izv. Akad Nauk SSSR, Mekhan, Zhidkosti i Gaza, No. 2,3-g, 1966.
6.
STULOV, V. P., and TURCHAK, L. I., Supersonic flow past blunt bodies in the presence of rapid nonequilibrium processes, Izv. Akad. Nauk SSSR, Mekhan. Zhidkosti i Gaza, No. 5,88-93,1967.
7.
KRAIKO, A. N., Numerical integration of equations In which the derivative contains a smaII parameter, Zh. vj%hisLMat. mat. Fiz.., 9, No. 2,438-442,1969.
8.
MOROZOV, I. I., and GERLIGA, V. A., Stability of boilers Wstoichivost’ kipyashchikh apparatov), Atomizdat, Moscow, 1969.
9.
GORYACHENKO, V. D., Methods of stability theory in the dynamics of nuclear reactors (Metody teorii ustoichivosti v dmamike yademykh reaktorov), Atomizdat, MOSCOW, 1971.
10. S~ARSKII,
A. A., rn~od~ction to the theory of difference schemes (Vvedenie v teoriyu raznostnykh
skhem),Nauka,MOSCOW, 1971. 11. BORISENKO,A. I., Gusdynamicsof engines (Gazovayadinamikadvigateleil,Oborongia,Moscow,1962.
THE CONSTRUCTION OF STABLE DIFFERENCE SCHEMES FOR SECOND-ORDER ALTERNATING LINEAR DIFFERENTIAL OPERATORS* Y-U.A. EREMIN and A. S. LEONOV Moscow (Received
17 July 1973 ; revised 15 October 1973)
AH AWORITHM for constructing stable difference schemes for second-order alternating linear differential operators is described and proved; it is based on the method of Tikhonov realization and the generalized disparity principle. An estimate is obtained for the accuracy of the approximate solution. The method is illustrated by a numerical example. The problem of constructing stable difference schemes for alternating linear differential *Zh. $ihisL Mat. mat. Fiz., 15,3,635-643,
1975.
87
Second-order alterhating linear differential operators
operators has still not been solved. At the same time, a considerable number of problems of practical interest, such as those concerning electromagnetic wave propagation in a waveguide, or diffraction at a local body of arbitrary filling, etc., lead to the requirement for constructing stable difference schemes for weakly elliptic operators, which are not in general of f=ed sign. In this connection, the present paper describes and proves a new and fairly general approach to this problem. The underlying idea is to use an arbitrary (possibly unstable) difference scheme, approximating the initial boundary value problem, in order to construct a stable algorithm for the approximate solution of the problem. The construction of the algorithm is based on the method of Tikhonov regularization, with a constructive choice of the regularization parameter based on the generalized disparity principle introduced in [I] . Our algorithm can be efficiently realized numerically; this is illustrated by the example of solving a boundary value problem for the Helmholtz equation in a circle. Estimation of the accuracy of the approximate solution obtained by our method is also discussed.
I
Suppose we are given the linear, in general alternating, operators
with domains of definition D(S?‘,) =D(5?,,)=D(_P) =Wzl(IIIn), 122, n<21, where IL= {s= (Si, . . . , Sn): S,
(1) where ZED (P), u,=R($? 1) , and u2 is the trace of a function uz*=R (9”) on the boundary of the domain II,. We shall write the problem (1) for brevity in the operator form L?z=u,
(2)
where z belongs as before to I%‘,‘[ II,], 9= (Zi, LZo1seann),and the vector rz= (u,, ~2) belongs to a new space U=L2 (II,) XL2(a&),the norm in which is given by the relation
Let the problem (2) have a unique solution TED (9)) functions Li= (U’,, ii*) EU. Further, let ZNi and U,,, i=l, Euclidean Ni-dimensional vector spaces.
corresponding to the collection of 2, . . . , n,-- be finite-dimensional
We consider spaces .Z,n and UN”, N= (N,, . . . , NN) , which can be written as the direct products Z,n=ZN,XZN,X . . . XZN, and U,n=U,,XU,,X.. . XUNn. Henceforth we shall denote elements of spaces Z# and UNn by the symbols 2, and fi,, and their norms by lbN1l~N and IIdavIIwn. We introduce the operators (PN: we' operators, and the interpolation operators QN
-+ZN",
$N
: ZNn*Wz'(II,),
: U+UN",
which will be called drift It will be assumed
ijjN : UNn+U.
88
Yu. A. Eremin and A. S. Leonov
1) the drift operators are linear and bounded, 2) for any .& =ZNn and ii,~lJ,n (3) (the constants C, and C, are independent of iN and &), 3) on the elements rir and Uz’6Lz (II,) we have the limit relation lim il~Nqh~--uIIv=o; N-PC.3 4) for any zEvz’ (IL)
Il~~Zllz~~~llzIIw,~~~,,,
m=const.
Conditions 2) and 4) are conditions for matching of the norms in I$‘: (II,) and ZNn. Notice that we can choose e.g., simple drift operators [2] as (PN and $N. The construction of the interpolation operators was described in [3]. We consider next the finite-dimensional linear operator AN with domain of defmition
D(A,)=ZNn and range OfVahleSR(A,)GJ,n.We shall call AN the approximating operator for the operator 9 of the initial problem on the element Z*ED (9) , provided that the measure of the approximation
4, (z*) ~ll9~*-$,vAN(~Nz*llv
(4)
tends to zero asN+ m, i.e., as N1, . . . , N,+m. We shall assume that the approximation is such that, on elements z&(P)
where XN,is independent of z and xN+O
asN+
00.
Instead of problem (2), let us consider the problem &^Z,=&,
b
with the operator approximating 9 of problem (2) in the sense of (4), (5), and with right-hand side satisfying the condition
Then, instead of the collection of exact input data (II 1, .!I*)we are given approximate input data of the problem (2): p= {N, &fn,
uNn,
qN,
$N,
TN,
$N,
AN,
UN,
br 6,
XN),
of which it is required to construct the framework of th_”approximate Solution [2] zIy pEZ,” of the problem (2), in such a way that the element zr=(pNzN. PEW*’ (II,,) is close to ZmthemetricofCasN+-and6+0. on
the
basis
89
Second-order alternating linear differential operators
2
A stable difference scheme for solving the above problem can be constructed by using a regularizing algorithm, based on minimization of the Tikhonov functional
followed by a choice of the regularization parameter a on the basis of the generalized disparity principle: if Z$ p is the extremal minimizing the functional (7) in Z# for the given (Y> 0 and fixed p, then the required o(p, e) is the solution of the equation
where J+,e is the following quantity, found to an accuracy E: (9) i.e., O~y,"--pC1,~~. The basis for this method, in the case of incorrectly posed (unstable) problems of a more general hind, was provided in [3]. The main results of [3] remain unchanged in the case of problem (2) and are stated in the form of theorems below. Theorem
1
For any (Y> 0 and a futed collection p, a unique extremal ;L”N,,, exists, minimizing the functional (7) in Z# . Theorem 2 For all 6~ (0, &I, EE (0, ~1,
N~No let (10)
Then Eq. (8) has a unique solution a (p, E) >O. Theorem
3
Let condition (10) hold for the sequences{&}, tend to zero as k -+ 00.Then,
,wk,ek)
where s N k,pk of Eq. (8).
.
IS the
{ek}, {UN:},
i=l,
2, . . . , n,which
extremal of the function (7), corresponding to the solutiona=
a (pk, --.sk)
Notice also that, if the solution of problem (2) is not unique, with u = ii the sequence of -a(P,+,~k) functions ($NkzNk,pk } converges in the metric of C to the normal solution of problem (2) (see 141).
90
Yu.A. Eremin and A. S. Leonov
The results stated in Theorems l-3 enable us to state the following stable algorithm for solving problem (2): 1) there exists to an accuracy E a number pP, which is the solution of the extremal problem (9); methods for solving such extremal problems are well developed [S] ; 2) for every fured collection p and every (Y> 0, there exists a unique extremal of the functional (7), convergent to the solution of Euler’s equation
(uTN+A,+As)^ZN=AN+^UN,~,
(11)
where TN is a finite-dimensional linear operator, acting from Z,vn into ZNn, the concrete form of which is determined by the normalization of the space Z# ; 3) the value of the regularization parameter cu(p.e) is found from Eq. (8); the function zr (s) =q~%,E’~)is taken as the approximate solution of problem (2). It must be noted that, if, during the computations, it turns out that ~~“<6~, then, see [I] , condition (8) for choosing (11 may be modified as follows: p(a) -IlnN~i&-&,,ll;~
-
(s+xNll&~pllz~J”=o.
To perform the algorithm effectively, we require upper and lower bounds for the parameter c&, E). As the lower bound we can take the value (II* = 0. The problem of obtaining the upper bound (Y**is discussed in [6,7]. In short, our algorithm for approximately solving problem (2) can be effectively realized numerically, and Eq. (1 I), with the parameter a! chosen in accordance with (8), represents the required stable difference scheme.
3 We shah now discuss the estimation of the accuracy of the approximate solution, in a given metric, when the algorithm described above is applied to problem (2). Let the operator 9,
act from Wz’ (IL)
into LZ(IL),
132,n<24 while
This relation holds e.g., if the spectral set g, corresponding to the operator pi and the boundary =O, is such that 1p 1>p'>O for arbitrary ~E%X (see [8]). condition L?,,z 1eeon, i.e., the We now construct the modulus of continuity of the inverse to p’I, in space Lz(II,), function o(r), continuous, non-negative, and monotonically non-decreasing for r > 0, such that w(O) = 0 and
Ilz’llL*,n,,-
(7))
provided that
lI%Z*lILbwI”)~~,
z-D(P).
91
Second-order alternating linear differential operators
Lemma
Let relation (12) hold, the number y* > 0 being known. Then. 0
(t) =7/p*.
(13)
The proof follows immediately from the defmition of the modulus of continuity. The next theorem estimates the accuracy of the approximate solution of problem (2). Theorem 4 Let
rp,^zyy’ be the approximate solution of problem (2) obtained by the method of Section
2, corresponding to the collection p of initial data. Then,
(14)
where w(r) is given by Eq. (13).
It is easily seen, from conditions (3) and (8) that
(16) Finally, by assumption (6), (17)
Using the extremal properties of the element i,?:*e’ we obtain
(18)
92
Yu. A. Eremin and A. S. Leonov
Adding (15)-( 17) in t’ e light of (1 g), we finally get
The estimate (14) now follows, on using the defmition of the modulus of continuity o(r). Notice that, for many practical problems, the number P* is either known? or can be found approximately.
4 As an example of the application of our algorithm, we take the boundary value problem of the 3rd kind for the well-known Hehnholtz equation. Given the rectangular domain IL= { (p, rp) : O&cp
k=const,
(19)
and on the boundary of the domain, the conditions
($- B”)1P=n 12(0, $4I <+w2
-
=u2,
p=const>o,
z (p. 0) =z
(P, 2n)
(20) ’
FIG. 2 FIG. 1 ) exact solution; o) approximate solution (mesh N = 36); x) mesh N = 144
93
Second-order alternating linear differential operators
Assuming that the constants k and /I are such that the resonant case is excluded, it can easily be seen that problem (19), (20) satisfies all the requirements imposed on problem (2) in Section 1. We can thus use the algorithm of Section 2 to obtain an approximate solution. The algorithm was programmed for the BESM-4 computer. The axisymmetric case with p1 c 0 and ~2 z 1 was considered. The one-dimensional boundary value problem then obtained was approximated by a homogeneous difference scheme. The quantities x and 6, required for numerical realization of the algorithm, were written in the explicit form. For instance, for the difference scheme
AN=
a,
a2 0 0
... 0
5 0
r2
r,
0
... 0
0
0
rp
rg r6
... 0
0
0
0
. . .
.
.
0
0
. . . . . . . . . . . . . . . . . . 0
0
r3N_g
0 0 0 0 ... 0 (where di=k2/4-l/h’, (1-0.5/i)hm2, i=l,
dz=lIh2,
rs+l=k2-2/h2,
2, . . . ,
0
&=-l/h,
r3N-7
:rgN_B
d3
dp *
&=1/h-_B
rsi=(1+0.5/i)hw2,
&,a=(O,
N-2, h=al (N-l))they are equal to
xN=h/24,
(1+0.5hla)
-0.5h12,
psi-z=
0,. ..,0,If 0.5hla), 6=0.5h.
The computations were performed for different values of k with a = 1 and with different numbers of mesh base-points. The exact solution of problem (19) (20) is known to be
Z(p,cF)= -
Jo(kp) kJ, (k) +BJo(k) ’
The computational results are illustrated by the curves of Figs. l-3. In Figs. 1 and 2 we give the results for k = 12 and 18 respectively. While on the coarse mesh the error of the approximate solution is 8%, on the mesh with 144 points it is under 1%. For comparison, in Fig. 3 we give the results for k = 12 and 288 mesh points, computed by the usual method, i.e., with a! = 0. The curve of the exact solution is shown for comparison on all the figures.
FIG. 3. -)
exact solution; x) approximate solution; 9. . ) solution with (Y= 0
94
Yu. A. Eremih and A. S. Leonov
Some general comments may be made regarding the numerical realization of our algorithm. The average number of iterations with respect to (IIfor finding the root of Eq. (8) was not greater than 10. The approximate solution was assumed to have been reached here if no further refinement of the generalized disparity occurred. For the 288-point mesh, the time for one CIiteration was 10 sec. Notice also that, due to the presence of the estimate (14), the required accuracy for the approximate solution of the problem can be achieved by choosing a finer mesh. In conclusion, let us comment on the scope for using our method. The greatest difficulty, in the sense of computer time, is presented by the realization of Para. 2) of the algorithm of Section 2, when solving problems with n > 2 dimensions. For solving Eq. (1 l), therefore, the most economic methods need to be used. It would seem that these are iterative methods [9]. The authors thank A. N. Tikhonov, A. G. Sveshnikov, and V. B. Glasko for guidance, and also A. V. Goncharskii, A. G. Yagola and Yu. M. Anisimov for valuable discussions. Danslated by D. E. Brown REFERENCES 1.
GONCHARSKII, A. V., LEONOV, A. S., and YAGOLA, A. G., The generalized disparity principle, Zh. vjchisl. Mat. mat. Fiz., 13, No. 2, 294-302,1913.
2.
GAVURIN, M. K., Lectures on computational methods (Lektsii po metodam vychislenii), Nauka, Moscow, 1971.
3.
GONCHARSKII, A. V., LEONOV, A. S., and YAGOLA, A. G., Finite-difference approximation incorrectly posed problems, Z/r. vjchisl. Mat. mat. Fiz., 14, No. 1, 15-24, 1974.
4.
TIKHONOV, A. N., On incorrectly posed problems of linear algebra and a stable method for solving them, Dokl. Akud. Nauk SSSR, 163, No. 3,591-594, 1965.
5.
BUDAK, B. hi., and VASIL’EV, F. P., Approximate methods for solvingoptimal control problems (Priblizhennye metody resheniya zadach optimal’nogo upravleniya), No. 2, 152-160, Izd-vo MGU, Moscow, 1969.
6.
GONCHARSKII, A. V., LEONOV, A. S., and YAGOLA, A. G., The solution of two-dimensional Fredholm integral equations of the 1st kind with a kernel dependent on the difference between the arguments, Z/J. vj%hisl. Mat. mat. Fiz., 11, No. 5, 1296-1301, 1971.
7.
VINOKUROV, V. A., Two remarks on the choice of regularization parameter, Zh. vjchisl. Mat. mat Fiz., 12, No. 2,481-483, 1972.
8.
AKHIEZER, N. I., and GLAZMAN, I. M., Theory of linear operators in Hilbert space (Teoriya lineinykh operatorov v gil’bertovom prostranstve), Nauka, Moscow, 1966.
9.
FEDORENKO, R. P., Iterative methods for solving elliptic difference equations, Us. mat. Nauk, 28, No. 2, 120-182. 1973.
of linear
Yu. A. Eremin and A. S. Leotwv
against the integrationstep for differentvaluesof the parameterS2= w L/c. The continuouscurves givethe errorobtainedwhen the scheme (2.4), (2.5) (u = %) is used, and the broken curves the error when the standard Runge-Kutta method with a constant step is used. It can be seen that the scheme (2.4), (2.5) gives high accuracy with a step over 10 times the maximum permissible step for the Runge-Kutta method. The authors thank E. F. Sabaev for su~esting the problem and for valuable co~ents. 2hnslated by D. E. Brown REFERENCES 1.
KATSKOVA, 0. N., and KRAIKO, A, N., Computation of plane and axisymmetric supersonic flows in the presence ofirreversible processes, PrtkLMot. tekhn, Fiz., No. 4,116-l 18,1963.
2.
KATSKOVA, 0. N., and KRAIKO, A. N., Computation of plune and axisymmetric supersonic flows in the presence of irreversibleprocesses, VTs Akad, Nauk SSSR, Moscow, 1964.
3.
GALYUN, N. S., and KRAIKO, A. N., Computation Mekhan ~Mashinos~, No. 6,41-47,1964.
4.
KAMZOLOV, V. N., and PIRUMOV, U. G., Computation of non-equIIibrIum flows in nozzles, Zzv. Akad Nauk SSSR, Mekhan Zhidkosti i Gaza, No. 6,25-33,1966.
5.
of non-equilibrium
flows, Izv. Akud Nauk SSSR,
GROMOV, V. G., Chemically non~qu~b~um boundary layer in dissociated air, Izv. Akad Nauk SSSR, Mekhan, Zhidkosti i Gaza, No. 2,3-g, 1966.
6.
STULOV, V. P., and TURCHAK, L. I., Supersonic flow past blunt bodies in the presence of rapid nonequilibrium processes, Izv. Akad. Nauk SSSR, Mekhan. Zhidkosti i Gaza, No. 5,88-93,1967.
7.
KRAIKO, A. N., Numerical integration of equations In which the derivative contains a smaII parameter, Zh. vj%hisLMat. mat. Fiz.., 9, No. 2,438-442,1969.
8.
MOROZOV, I. I., and GERLIGA, V. A., Stability of boilers Wstoichivost’ kipyashchikh apparatov), Atomizdat, Moscow, 1969.
9.
GORYACHENKO, V. D., Methods of stability theory in the dynamics of nuclear reactors (Metody teorii ustoichivosti v dmamike yademykh reaktorov), Atomizdat, MOSCOW, 1971.
10. S~ARSKII,
A. A., rn~od~ction to the theory of difference schemes (Vvedenie v teoriyu raznostnykh
skhem),Nauka,MOSCOW, 1971. 11. BORISENKO,A. I., Gusdynamicsof engines (Gazovayadinamikadvigateleil,Oborongia,Moscow,1962.
THE CONSTRUCTION OF STABLE DIFFERENCE SCHEMES FOR SECOND-ORDER ALTERNATING LINEAR DIFFERENTIAL OPERATORS* Y-U.A. EREMIN and A. S. LEONOV Moscow (Received
17 July 1973 ; revised 15 October 1973)
AH AWORITHM for constructing stable difference schemes for second-order alternating linear differential operators is described and proved; it is based on the method of Tikhonov realization and the generalized disparity principle. An estimate is obtained for the accuracy of the approximate solution. The method is illustrated by a numerical example. The problem of constructing stable difference schemes for alternating linear differential *Zh. $ihisL Mat. mat. Fiz., 15,3,635-643,
1975.
87
Second-order alterhating linear differential operators
operators has still not been solved. At the same time, a considerable number of problems of practical interest, such as those concerning electromagnetic wave propagation in a waveguide, or diffraction at a local body of arbitrary filling, etc., lead to the requirement for constructing stable difference schemes for weakly elliptic operators, which are not in general of f=ed sign. In this connection, the present paper describes and proves a new and fairly general approach to this problem. The underlying idea is to use an arbitrary (possibly unstable) difference scheme, approximating the initial boundary value problem, in order to construct a stable algorithm for the approximate solution of the problem. The construction of the algorithm is based on the method of Tikhonov regularization, with a constructive choice of the regularization parameter based on the generalized disparity principle introduced in [I] . Our algorithm can be efficiently realized numerically; this is illustrated by the example of solving a boundary value problem for the Helmholtz equation in a circle. Estimation of the accuracy of the approximate solution obtained by our method is also discussed.
I
Suppose we are given the linear, in general alternating, operators
with domains of definition D(S?‘,) =D(5?,,)=D(_P) =Wzl(IIIn), 122, n<21, where IL= {s= (Si, . . . , Sn): S,
(1) where ZED (P), u,=R($? 1) , and u2 is the trace of a function uz*=R (9”) on the boundary of the domain II,. We shall write the problem (1) for brevity in the operator form L?z=u,
(2)
where z belongs as before to I%‘,‘[ II,], 9= (Zi, LZo1seann),and the vector rz= (u,, ~2) belongs to a new space U=L2 (II,) XL2(a&),the norm in which is given by the relation
Let the problem (2) have a unique solution TED (9)) functions Li= (U’,, ii*) EU. Further, let ZNi and U,,, i=l, Euclidean Ni-dimensional vector spaces.
corresponding to the collection of 2, . . . , n,-- be finite-dimensional
We consider spaces .Z,n and UN”, N= (N,, . . . , NN) , which can be written as the direct products Z,n=ZN,XZN,X . . . XZN, and U,n=U,,XU,,X.. . XUNn. Henceforth we shall denote elements of spaces Z# and UNn by the symbols 2, and fi,, and their norms by lbN1l~N and IIdavIIwn. We introduce the operators (PN: we' operators, and the interpolation operators QN
-+ZN",
$N
: ZNn*Wz'(II,),
: U+UN",
which will be called drift It will be assumed
ijjN : UNn+U.
88
Yu. A. Eremin and A. S. Leonov
1) the drift operators are linear and bounded, 2) for any .& =ZNn and ii,~lJ,n (3) (the constants C, and C, are independent of iN and &), 3) on the elements rir and Uz’6Lz (II,) we have the limit relation lim il~Nqh~--uIIv=o; N-PC.3 4) for any zEvz’ (IL)
Il~~Zllz~~~llzIIw,~~~,,,
m=const.
Conditions 2) and 4) are conditions for matching of the norms in I$‘: (II,) and ZNn. Notice that we can choose e.g., simple drift operators [2] as (PN and $N. The construction of the interpolation operators was described in [3]. We consider next the finite-dimensional linear operator AN with domain of defmition
D(A,)=ZNn and range OfVahleSR(A,)GJ,n.We shall call AN the approximating operator for the operator 9 of the initial problem on the element Z*ED (9) , provided that the measure of the approximation
4, (z*) ~ll9~*-$,vAN(~Nz*llv
(4)
tends to zero asN+ m, i.e., as N1, . . . , N,+m. We shall assume that the approximation is such that, on elements z&(P)
where XN,is independent of z and xN+O
asN+
00.
Instead of problem (2), let us consider the problem &^Z,=&,
b
with the operator approximating 9 of problem (2) in the sense of (4), (5), and with right-hand side satisfying the condition
Then, instead of the collection of exact input data (II 1, .!I*)we are given approximate input data of the problem (2): p= {N, &fn,
uNn,
qN,
$N,
TN,
$N,
AN,
UN,
br 6,
XN),
of which it is required to construct the framework of th_”approximate Solution [2] zIy pEZ,” of the problem (2), in such a way that the element zr=(pNzN. PEW*’ (II,,) is close to ZmthemetricofCasN+-and6+0. on
the
basis
89
Second-order alternating linear differential operators
2
A stable difference scheme for solving the above problem can be constructed by using a regularizing algorithm, based on minimization of the Tikhonov functional
followed by a choice of the regularization parameter a on the basis of the generalized disparity principle: if Z$ p is the extremal minimizing the functional (7) in Z# for the given (Y> 0 and fixed p, then the required o(p, e) is the solution of the equation
where J+,e is the following quantity, found to an accuracy E: (9) i.e., O~y,"--pC1,~~. The basis for this method, in the case of incorrectly posed (unstable) problems of a more general hind, was provided in [3]. The main results of [3] remain unchanged in the case of problem (2) and are stated in the form of theorems below. Theorem
1
For any (Y> 0 and a futed collection p, a unique extremal ;L”N,,, exists, minimizing the functional (7) in Z# . Theorem 2 For all 6~ (0, &I, EE (0, ~1,
N~No let (10)
Then Eq. (8) has a unique solution a (p, E) >O. Theorem
3
Let condition (10) hold for the sequences{&}, tend to zero as k -+ 00.Then,
,wk,ek)
where s N k,pk of Eq. (8).
.
IS the
{ek}, {UN:},
i=l,
2, . . . , n,which
extremal of the function (7), corresponding to the solutiona=
a (pk, --.sk)
Notice also that, if the solution of problem (2) is not unique, with u = ii the sequence of -a(P,+,~k) functions ($NkzNk,pk } converges in the metric of C to the normal solution of problem (2) (see 141).
90
Yu.A. Eremin and A. S. Leonov
The results stated in Theorems l-3 enable us to state the following stable algorithm for solving problem (2): 1) there exists to an accuracy E a number pP, which is the solution of the extremal problem (9); methods for solving such extremal problems are well developed [S] ; 2) for every fured collection p and every (Y> 0, there exists a unique extremal of the functional (7), convergent to the solution of Euler’s equation
(uTN+A,+As)^ZN=AN+^UN,~,
(11)
where TN is a finite-dimensional linear operator, acting from Z,vn into ZNn, the concrete form of which is determined by the normalization of the space Z# ; 3) the value of the regularization parameter cu(p.e) is found from Eq. (8); the function zr (s) =q~%,E’~)is taken as the approximate solution of problem (2). It must be noted that, if, during the computations, it turns out that ~~“<6~, then, see [I] , condition (8) for choosing (11 may be modified as follows: p(a) -IlnN~i&-&,,ll;~
-
(s+xNll&~pllz~J”=o.
To perform the algorithm effectively, we require upper and lower bounds for the parameter c&, E). As the lower bound we can take the value (II* = 0. The problem of obtaining the upper bound (Y**is discussed in [6,7]. In short, our algorithm for approximately solving problem (2) can be effectively realized numerically, and Eq. (1 I), with the parameter a! chosen in accordance with (8), represents the required stable difference scheme.
3 We shah now discuss the estimation of the accuracy of the approximate solution, in a given metric, when the algorithm described above is applied to problem (2). Let the operator 9,
act from Wz’ (IL)
into LZ(IL),
132,n<24 while
This relation holds e.g., if the spectral set g, corresponding to the operator pi and the boundary =O, is such that 1p 1>p'>O for arbitrary ~E%X (see [8]). condition L?,,z 1eeon, i.e., the We now construct the modulus of continuity of the inverse to p’I, in space Lz(II,), function o(r), continuous, non-negative, and monotonically non-decreasing for r > 0, such that w(O) = 0 and
Ilz’llL*,n,,-
(7))
provided that
lI%Z*lILbwI”)~~,
z-D(P).
91
Second-order alternating linear differential operators
Lemma
Let relation (12) hold, the number y* > 0 being known. Then. 0
(t) =7/p*.
(13)
The proof follows immediately from the defmition of the modulus of continuity. The next theorem estimates the accuracy of the approximate solution of problem (2). Theorem 4 Let
rp,^zyy’ be the approximate solution of problem (2) obtained by the method of Section
2, corresponding to the collection p of initial data. Then,
(14)
where w(r) is given by Eq. (13).
It is easily seen, from conditions (3) and (8) that
(16) Finally, by assumption (6), (17)
Using the extremal properties of the element i,?:*e’ we obtain
(18)
92
Yu. A. Eremin and A. S. Leonov
Adding (15)-( 17) in t’ e light of (1 g), we finally get
The estimate (14) now follows, on using the defmition of the modulus of continuity o(r). Notice that, for many practical problems, the number P* is either known? or can be found approximately.
4 As an example of the application of our algorithm, we take the boundary value problem of the 3rd kind for the well-known Hehnholtz equation. Given the rectangular domain IL= { (p, rp) : O&cp
k=const,
(19)
and on the boundary of the domain, the conditions
($- B”)1P=n 12(0, $4I <+w2
-
=u2,
p=const>o,
z (p. 0) =z
(P, 2n)
(20) ’
FIG. 2 FIG. 1 ) exact solution; o) approximate solution (mesh N = 36); x) mesh N = 144
93
Second-order alternating linear differential operators
Assuming that the constants k and /I are such that the resonant case is excluded, it can easily be seen that problem (19), (20) satisfies all the requirements imposed on problem (2) in Section 1. We can thus use the algorithm of Section 2 to obtain an approximate solution. The algorithm was programmed for the BESM-4 computer. The axisymmetric case with p1 c 0 and ~2 z 1 was considered. The one-dimensional boundary value problem then obtained was approximated by a homogeneous difference scheme. The quantities x and 6, required for numerical realization of the algorithm, were written in the explicit form. For instance, for the difference scheme
AN=
a,
a2 0 0
... 0
5 0
r2
r,
0
... 0
0
0
rp
rg r6
... 0
0
0
0
. . .
.
.
0
0
. . . . . . . . . . . . . . . . . . 0
0
r3N_g
0 0 0 0 ... 0 (where di=k2/4-l/h’, (1-0.5/i)hm2, i=l,
dz=lIh2,
rs+l=k2-2/h2,
2, . . . ,
0
&=-l/h,
r3N-7
:rgN_B
d3
dp *
&=1/h-_B
rsi=(1+0.5/i)hw2,
&,a=(O,
N-2, h=al (N-l))they are equal to
xN=h/24,
(1+0.5hla)
-0.5h12,
psi-z=
0,. ..,0,If 0.5hla), 6=0.5h.
The computations were performed for different values of k with a = 1 and with different numbers of mesh base-points. The exact solution of problem (19) (20) is known to be
Z(p,cF)= -
Jo(kp) kJ, (k) +BJo(k) ’
The computational results are illustrated by the curves of Figs. l-3. In Figs. 1 and 2 we give the results for k = 12 and 18 respectively. While on the coarse mesh the error of the approximate solution is 8%, on the mesh with 144 points it is under 1%. For comparison, in Fig. 3 we give the results for k = 12 and 288 mesh points, computed by the usual method, i.e., with a! = 0. The curve of the exact solution is shown for comparison on all the figures.
FIG. 3. -)
exact solution; x) approximate solution; 9. . ) solution with (Y= 0
94
Yu. A. Eremih and A. S. Leonov
Some general comments may be made regarding the numerical realization of our algorithm. The average number of iterations with respect to (IIfor finding the root of Eq. (8) was not greater than 10. The approximate solution was assumed to have been reached here if no further refinement of the generalized disparity occurred. For the 288-point mesh, the time for one CIiteration was 10 sec. Notice also that, due to the presence of the estimate (14), the required accuracy for the approximate solution of the problem can be achieved by choosing a finer mesh. In conclusion, let us comment on the scope for using our method. The greatest difficulty, in the sense of computer time, is presented by the realization of Para. 2) of the algorithm of Section 2, when solving problems with n > 2 dimensions. For solving Eq. (1 l), therefore, the most economic methods need to be used. It would seem that these are iterative methods [9]. The authors thank A. N. Tikhonov, A. G. Sveshnikov, and V. B. Glasko for guidance, and also A. V. Goncharskii, A. G. Yagola and Yu. M. Anisimov for valuable discussions. Danslated by D. E. Brown REFERENCES 1.
GONCHARSKII, A. V., LEONOV, A. S., and YAGOLA, A. G., The generalized disparity principle, Zh. vjchisl. Mat. mat. Fiz., 13, No. 2, 294-302,1913.
2.
GAVURIN, M. K., Lectures on computational methods (Lektsii po metodam vychislenii), Nauka, Moscow, 1971.
3.
GONCHARSKII, A. V., LEONOV, A. S., and YAGOLA, A. G., Finite-difference approximation incorrectly posed problems, Z/r. vjchisl. Mat. mat. Fiz., 14, No. 1, 15-24, 1974.
4.
TIKHONOV, A. N., On incorrectly posed problems of linear algebra and a stable method for solving them, Dokl. Akud. Nauk SSSR, 163, No. 3,591-594, 1965.
5.
BUDAK, B. hi., and VASIL’EV, F. P., Approximate methods for solvingoptimal control problems (Priblizhennye metody resheniya zadach optimal’nogo upravleniya), No. 2, 152-160, Izd-vo MGU, Moscow, 1969.
6.
GONCHARSKII, A. V., LEONOV, A. S., and YAGOLA, A. G., The solution of two-dimensional Fredholm integral equations of the 1st kind with a kernel dependent on the difference between the arguments, Z/J. vj%hisl. Mat. mat. Fiz., 11, No. 5, 1296-1301, 1971.
7.
VINOKUROV, V. A., Two remarks on the choice of regularization parameter, Zh. vjchisl. Mat. mat Fiz., 12, No. 2,481-483, 1972.
8.
AKHIEZER, N. I., and GLAZMAN, I. M., Theory of linear operators in Hilbert space (Teoriya lineinykh operatorov v gil’bertovom prostranstve), Nauka, Moscow, 1966.
9.
FEDORENKO, R. P., Iterative methods for solving elliptic difference equations, Us. mat. Nauk, 28, No. 2, 120-182. 1973.
of linear