- Email: [email protected]
An unsaturated method for the discretization of systems of linear differential equations unsolved with respect to derivatives
100
5. 6. 7. 8. 9.
10. 11.
12.
13.
Nauka, Moscow, 1980. LAX P.D., Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves. Sot. Industr. Appl. Math., Philadelphia, 1972. LAX P. and WENDROFF B., Systems of Conservation Laws. Communs. Pure Appl. Math., 13, 217-237, 1960. HARTEN A., HYMAN J.M. and LAX P.D., On finite-differenceapproximations and entropy conditions for shocks. Communs. Pure Appl. Math., 29, 297-322, 1976. OSTAPBNKO V.V., On equivalent definitions of the concept of conservativenessfor finite-differenceschemes. Zh. Vychisl. Mat. Mat. Fiz., 29, 8, 1114-1128, 1989. OSTAPENKO V.V., A method for the theoretical estimation of imbalances of nonconservative difference schemes on shock waves. Dokl. Akad. Nauk SSSR, 295, 2, 292-297, 1987. HARTEN A., High resolution schemes for hyperbolic conservation laws. J. Comput. 1983. Phys., 49, 357-393, OSTAPENKO V.V., On the standard approximation of differential operators. Dinamika Sploshnoi Sredy, Inst. Gidrodinamiki Sibirsk. Otdel. Akad. Nauk SSSR, Novosibirsk, 72, 67-83, 1985. KOPCHENOV V.I. and KRAIKO A.N., A second-order monotone difference scheme for hyperbolic systems with two independent variables. Zh. Vychisl. Mat. Mat. Fiz., 23, 4, 848-859, 1983. HARTEN A., On a class of high resolution total-variation-stablefinite-difference schemes. SIAM J. Numer. Anal., 21, 1, l-23, 1984. Translated by D.L.
U.S.S.R. Comput.Maths. MatkPhys., Printed in Great Britain
Vo1.30,No.5,pp.100-106,199O
0041-5553/90 $10.00+0.00 01991 Pergamon Press plc
SHORT COMMUNICATIONS AN UNSATURATED METHOD FOR THE DISCRETIZATION OF SYSTEMS OF LINEAR DIFFERENTIAL EQUATIONS UNSOLVED WITH RESPECT TO DERIVATIVESs P.N. DENISENKO
1.
The method is based on the replacement of a system of differential equations by an equivalent system of integral Volterra equations of the third kind and the approximation in the system obtained of the operator of integration with variable upper limit by the interpolationoperator (with Chebyshev interpolationnodes of the form (cos(ti/n))on in the interval of integration) of numerical integration. After approximation a system of algebraic linear equations is obtained. The approximation, stability, and convergence of the method are investigated.The program realization of the method is described and its computational complexity is analysed. The results of the analysis of the computational experiment are presented where the program realizes the method on model problems of various types illustrating the effectiveness of the method. Consider the Cauchy problem for the system of m differential equations A(z)Y'(I)=B(z)Y(z)+D(z), Y(O)=Yo,
(1)
where I=[O,H],A(Z)and E(m) are matrix functions, and Y(s) and D(z) are, respectively, unknown and given vector functions. Suppose that in the interval of integration the matrix A(X) is continuously differentiable and the matrix B(x) and the function D(x) are continuous. Then this system is equivalent to the system of integral Volterra equations of the third kind
qyw, where
~l-w+,
(2)
101
R[Y(z),+]=A(s)Y(t)-
j
C(t)-d?(r)+A'(z),
c(t)Y(t)dt,
GM)== U(L)dt+A(O)Yp. P
In system (2) we approximate the integration operator i * operator 1% of numerical integration (see [ll):
in the interval [O, HI by the
(3) where L
-
H
z
R
%f(q):
=LIf(~),+tl,
=,=(I-Bf)H/2,
zi=eos(l~/n),
n
i-0,i ,...,
n,
1-o
are Chebyshev nodes of the second kind in the interval [O, HI, ” I, a~=“&?+a,f . . . + a._l+aJ2. c f-e = of the operator In are defined as the integrals The coefficients of the matrix {arj)o of fundamental interpolation polynomials on Chebyshev nodes ZI “‘j==--
-
“,s
4 (0 6
1,t-0
* 2*....4
(4)
*
where, according to [21, p.94,
and Tk(t) are Chebyshev polynomials. Explicit formulas for the coefficients are obtained in 131. After approximating the integration operator in (2) by the operator In we obtain a system of linear algebraic equations of order n-m = N of the following form:
(5) where
D,~n(rr)-(aI(zi),...,
d&c)); &-A(q), C&'(+,) are square matrices of order rrt, i = O,l,. . ..n. *
m
)I
The unknown variable in this system is the net vector function
t-i, 2,...,n, Y<=~Ylt,...,!flnc~, approximating the required solution Y(x) at the nodes of the net (XS)P.In this way we obtain a method of discretizing 14, Chapter 71 the Cauchy problem (1). 2. The approximation of the original system of linear differential equations (1) by the system of linear algebraic equations (5) can be estimated as the residue in the approximation of the integral operator P[U(zf,
zl:=K[U(zf,
SW(Z)
(6)
(here the operator K and the function G are defined in (2)) by the operator on net functions with Chebyshev nodes (3)
102
Fd (U,ll”, z,]=&[
(uj)~“,
(7)
a-G,
are defined where Uj-= U(zj), j = 1,2,...,n and the operator K, and the net function {G,),” The residue of approximation of the operator F by the operator Fn on the function in (5). 'J(a)@(F) is taken to be equal to the net function
i=l, 2,..., R. rin(~(z)l=PfU(~),rcl-Fn~(U(~~)}t", al, Lem
I.
The residue of approximation of the operator F by the operator F,
(8) on the
function V(x) is equal to the residue of approximation of the integration operator i 0 the numerical integration operator In on the function
by
R(z):=R[U(z)]=(A’(.r)+B(z))U(z)fD(z),
namely, we have the equality
7
rr"(U(s)l=-P,,[ ,R(r)].
i=i,Z,...,n.
0
Proof.
By the definition of the operators F in (6) and Fn in (7) we have
CorotZary 2. If in the Cauchy problem (11 the matrix A(x) is continuously differentiable, the matrix B(x) and function D(x) are continuous, then the condition of approximation (4, p.4451, of the operator F by the operators Fn is satisfied in the space of continuous functions C[O, H] and the method of discretization (5) is unsaturated. The residue of the Chebyshev integration process on a continuous function Proof. f(x) does not exceed
where EnIf] is a quantity of the best uniform approximation of the function f(x) by algebraic polynomials of order n, and the uniform norm of the integration operator in uniformly bounded in n (see (11): IlMicp, x]=H[1+0fn-‘Inn)].
Consequently, the order of decay of the residue of the process of discretization (5) increases together as the regularity of the solution Y(r) and functions A(r), B(r), D(X) increases, i.e., this process is unsaturated. 3. The stability of the discretizationmethod (5), unlike that of finite-difference methods, can be investigatedmore naturally in terms of general theory of abstract approximation schemes (4, p.2971. An investigationof the stability of the method in the case of the Cauchy problem for ordinary differential equations solved with respect to derivatives is carried out in (51 where the approximation of the integration operator by the polynomial analogue (11 of the operator of numerical integration (3) is considered. The results obtained in that paper can be extended in a natural way to the Cauchy problem for the system of differential equations (1) solved for the derivatives (A(Z)aE={8ij}tlli). The solvability condition for system (5) in this particular case (by the method of simple iteration) takes the form
where
(((fi),m((, is the length of the vector in the space Em. When this condition is satisfied it is obvious that the following conclusions hold: 1) the matrix K, of the system (1) is continuously invertible in the space of net vector-functions{F(z~)},* with the norm
103
2)
the norm of the inverse matrix is uniformly bounded in n,nz=no: ilK~-'iLc(l-q)-',
and, consequently, the stability condition is satisfied [4, p.2971. Suppose that the following conditions are satisfied for system (1) in the Theorem 1. interval [O, HI : 1) the matrix A(X) is continuously invertible, 2) the matrices A-~(r),A’(+) and B(z) are continuous. Suppose that the length of the integration interval does not exceed the quantity q
Proof. We consider the scalar case (m = 1). By the condition of the theorem, the function kl(s) exists and is continuous. Therefore, in this case system (5) is equivalent to the following: Yp.4
(91
-‘(M
The solvability of system (9) by the method of successive approximations obviously occurs in the case when the following inequality is satisfied: IIA-‘(~)licllz”llcllE(~)llc=q,.;;q<1. By the condition of the theorem this inequality occurs if
According to the estimate of the norm of the numerical integration operator [l] the last inequality holds for asno( Moreover, the sequence of norms of matrices inverse to the matrices of system (8) is uniformly bounded:
which means that the sequence of norms of the matrices inverse to the matrices of system (5) is also bounded:
The theorem in the vector case (m t 2) can be proved analogously. 4. The convergence of the sequence of solutions {(Y&n)~zi of the discrete problem (5) to the solution Y(Z) of the initial Cauchy problem (1) follows from the approximation and stability [4, p.2981. Moreover, the order of the error of the method
is equal to the order of approximation of the integral equation (2) by system (51, i.e., to the order of the residue (~P[Y])~,, of the approximation (8) of the operator F from (6) by the operator Fn from (7) on the exact solution Y(X) of Cauchy problem (1). Theorem 2. If the solution of Cauchy problem (1) exists and is unique and the operator X, in system (5) is continuously invertible, then the following equalities hold: (e~")l"=Kn-'(r,nIY(z)l),",
(10)
Ile”ll=C,llr”ll,
ill)
(IIK*II)-‘c;C,r:IIR*.-‘II.
Proof. We consider the integral equation (2) in the nodes of the net {z,),". Subtracting from each row of the system of equations obtained the corresponding row of system (5) we
104
obtain the following system of equations i31,2,.... K[Y(r), zr]-K,[(Yj)I", XiI==G(~I)--Gil Adding and subtracting from the left-hand side of the equations of this system the zi]and transferring the difference MY(r),+~I-&I(Y(~J)}I”, 4 over to the quantities K,[{Y(zj))P, right-hand side we obtain the equivalent system i! x~[{e"(r,)),,,t,l=r"[~ ,R[Y(z)l), i=1,2,... .
(10')
0 By the condition of the theorem, the operator Kn (the matrix of system (5)) is continuously invertible. Consequently, system (10') can be solved and its solution can be written in the form of Eq. (10). The estimate of the operator parts of systems (10) and (10') in the norm /I*[ I,, of the space RNcompatible with C[O, HI is a uniform norm of m-dimensional space of functions and gives Eq. (11). CorotZary. If the sequence of norms of the inverse matrices of system (5)
where E,[l?[Y]Jis the quantity of the best approximation of the function R[Y(xll by algebraic polynomials of order n, then the Chebyshev method (5) of the discretization of the Cauchy problem (1) is convergent. 5. The program realization of the discretizationmethod (5), (4) and (2) was executed as a consecutive realization of the following algorithms: 1) computation of the coefficients of the matrix of the numerical integration operator (3); 2) computation of the coefficients of system (5); 3) the solution of system (5) by Gauss's method without choosing the principal element. The computations by these algorithms requires that ~(nr)+O(nnz)+O((mn)s)=O((nan)g)is satisfied for arithmetic operations and o(n)+O(mQ~)=O(m*n) for the operation of computing functions. We will estimate the number of operations required by the method for the class of Cauchy problems of the form (1) satisfying the conditions of Theorem 2 with matrices A(x) and B(m), the function F(x), and the solutions belonging to the class of analytic functions in the ellipse of the complex plane with foci in the ends of the integration interval. According to Theorem 1 the number of net nodes depends on the required accuracy in the following way: n=O(lne) Hence, the dependence of the number of operations on the accuracy of the solution is O((rnhl8)~)for arithmetic operations and O(mslne)for the operations of computing functions. In the case of classical numerical methods of order p the dependence of the number of net nodes (and the number of arithmetic operations and operations of computing functions) on the accuracy has higher order o(mt+'p). 6. A numerical experiment to verify the effectiveness obtained in the assertions of Sections 2-5 was performed on a computer. The problems presented in [S-l03 and a whole series of equationsof the form(l) were used as model ones. An analysis of the results of the numerical experiment completely confirmed the effectiveness of the assertions obtained in Sections 2-5. The results of the analysis for each of the classes of differential equations is presented below. 7. The analysis of the results of a numerical experiment for classical nonstiff Cauchy problems with a unique solution belonging to the class of analytic functions 1 7 3 5 ii enabled us to draw the following conclusions. FIG. 1. The error of the method for the Cauchy problem is en and the residue of the approximation of integration operator (2) on the SXSCt SOlUtiOtl Of this problem is Pn and decreases as o($") when the number of nodes
105
increases. The decrease continues up to some value n,=no(e ) determined by the camp characteristicsof the computer (the length of computer words). Any further increase in the number of nodes n>no does not in practice lead to any change in the error of the method and the residue of the approximation. does not in practice depend on n. 2. For a number of nodes n
Y
(0) = I, y(O)=l,
ZE
(0,IOm],
.z=[O,
lam], m=O, i.
(12) (13)
Graph 1 shows the error of the method for problem (12) when H = 10; graph 2 shows the residue rn of integration of the right-hand side of Eq. (12) for H = 10; graph 3 shows the error of the method in problem (13) for H = 10; graph 4 shows the residue rn of integration of the right-hand side of (13) for H = 10;
graph 5 shows the error of the method and
the residue of integration of the right-hand side of problem (12) for H = 1 (in this case graph 6 shows the error of the method and the residue of the integration of the lle"ll/llr"ll~i); right-hand side of problem (13) and the model problem of [91 for H = 1. In the same figure (graph 7) we show the results of the solution of the Cauchy problem of type (13) in the interval LO, 11 by the Adams-Bashforthmethod presented in 191. A comparison of graph 7 with graph 6 well illustrates the above conclusions on the higher numerical stability of the proposed method. Using graphs 6 and 7 it is easy to compare the economy of the methods. Thus, to find the solution of the problem [91 with error Il~,li
the decrease of the ratio
lie,ll/llr,ll
as the stiffness of the system increases
up to a certain value lIBIIHza-l; b) the error of the met~~~pis independent of the length of the integration interval for II>H,,where ~~ze-1 camp /llBll. These conclusions are well illustrated by the results presented below on the solution of the following model problems: from [7, 101, y’=-l@y+f(s),
y(O)=&
(14)
where f(m) = (lOk - l)exp(-m), k = 0,...,12; from [a], I-2exp(-2) .zexp(-2) Y'(z)=lO‘ Y(s)+F(z), Y(O)= I II -sesp(-s) I-2exp(-s) II1' II
II
(15)
where F(x) is chosen so that Y(S) = (exp(z), exp(-z)), and also the one-dimensionala alogue of this system. These model problems are solved in the interval [O, HI, where H = 10't , 1 = -2,...,12. The results show that if for problem (14) we have H < 100, then Ilenll/llr,ll=O(lOLH)-l,n~lO, and if H 2 100, then lle.ll=lOW. For problem (15) and its one-dimensionalanalogue if H 2 100, then IleJ=l. 9. An analysis of the results of a numerical experiment on well-posed Cauchy problems for systems of the form (l), unsolved for the derivative, in particular on the problem
106
showed that such systems can be solved as well by the Chebyshev method the systems solved for the derivatives. 10. An analysis of the results of numerical experiment on well-posed Cauchy problems for systems of the form (1) with singularitiesof nodal type, in particular on the problem 2zy’=ky-f(r),
=P,
11,
v(O)=k-‘f(O),
(16)
where k = 0, 1. 1.5, 2-10-p, 109, p = 1,...,12, q = 1 ,...,4, fir) = (2% + kfexp(-zz),and the initial point is singular, enabled us to draw the following conclusions: 1) the solution by the Chebyshev method is unique; 2) this solution converges, as a rule, to the most regular of exact solutions of the Cauchy problem; 3) the ratio ~~~~/~r~~ essentially depends on which type the singular point belongs to. In particular, the error of the solution of problem (16) by the Chebyshev method depends on the parameter k in the following way: if k increases from 0 to 2, then the errors lieloll increase from 10-la to 10, and if k increases from iO*to IO’, the errors lleloll decrease from IO-' to fO+O. REFERENCES DENISENKO P.N., Chebyshev integration processes, in: Proc. Int. Conf. on Approximation Theory, 131-135, Nauka, Moscow, 1987. 2. PASHKOVSKII S., Numerical Applications of Chebyshev Polynomials and Series, Nauka, Moscow, 1983. DEYAD~ V.K., An approximation-iterativemethod of approximating solutions of 3. Cauchy's problem for ordinary differential equations, Preprint 84-27, Inst. Matem. Akad. Nauk UkrSSR, Kiev, 1984. 4. TRENOGIN V.A., Functional Analysis, Nauka, Moscow, 1980. 5, DENISENKQ P.N., To the solution of the Cauchy problem by interpolationon quasiChebyshev nodes, Differents.Uravn., 5, 5, 908-913, 1979. HAMMING R.V., Numerical Methods, Nauka, Moscow, 1968. 6. 7. RAKITSKII YU.V., USTINOV S.M. and CHERNORUTSKII I.G., Numerical Methods of Solving Stiff Systems, Nauka, Moscow, 1979. 8. HALL G. and WATT J.M. (Eds.)., Modern Numerical Methods for Ordinary Differential Equations, Clarendon Press, Oxford, 1976. BABUSKA I., PRAGER M. and VITkEK E., Numerical Processes in Differential Equations, 9. Wiley, London, 1966. 10. GEAR C.W., Numerical Initial Value Problems in Ordinary Differential Equations, Prentice Hall, Englewood Cliffs, 1971. 1.
Translated by M.C.
5. 6. 7. 8. 9.
10. 11.
12.
13.
Nauka, Moscow, 1980. LAX P.D., Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves. Sot. Industr. Appl. Math., Philadelphia, 1972. LAX P. and WENDROFF B., Systems of Conservation Laws. Communs. Pure Appl. Math., 13, 217-237, 1960. HARTEN A., HYMAN J.M. and LAX P.D., On finite-differenceapproximations and entropy conditions for shocks. Communs. Pure Appl. Math., 29, 297-322, 1976. OSTAPBNKO V.V., On equivalent definitions of the concept of conservativenessfor finite-differenceschemes. Zh. Vychisl. Mat. Mat. Fiz., 29, 8, 1114-1128, 1989. OSTAPENKO V.V., A method for the theoretical estimation of imbalances of nonconservative difference schemes on shock waves. Dokl. Akad. Nauk SSSR, 295, 2, 292-297, 1987. HARTEN A., High resolution schemes for hyperbolic conservation laws. J. Comput. 1983. Phys., 49, 357-393, OSTAPENKO V.V., On the standard approximation of differential operators. Dinamika Sploshnoi Sredy, Inst. Gidrodinamiki Sibirsk. Otdel. Akad. Nauk SSSR, Novosibirsk, 72, 67-83, 1985. KOPCHENOV V.I. and KRAIKO A.N., A second-order monotone difference scheme for hyperbolic systems with two independent variables. Zh. Vychisl. Mat. Mat. Fiz., 23, 4, 848-859, 1983. HARTEN A., On a class of high resolution total-variation-stablefinite-difference schemes. SIAM J. Numer. Anal., 21, 1, l-23, 1984. Translated by D.L.
U.S.S.R. Comput.Maths. MatkPhys., Printed in Great Britain
Vo1.30,No.5,pp.100-106,199O
0041-5553/90 $10.00+0.00 01991 Pergamon Press plc
SHORT COMMUNICATIONS AN UNSATURATED METHOD FOR THE DISCRETIZATION OF SYSTEMS OF LINEAR DIFFERENTIAL EQUATIONS UNSOLVED WITH RESPECT TO DERIVATIVESs P.N. DENISENKO
1.
The method is based on the replacement of a system of differential equations by an equivalent system of integral Volterra equations of the third kind and the approximation in the system obtained of the operator of integration with variable upper limit by the interpolationoperator (with Chebyshev interpolationnodes of the form (cos(ti/n))on in the interval of integration) of numerical integration. After approximation a system of algebraic linear equations is obtained. The approximation, stability, and convergence of the method are investigated.The program realization of the method is described and its computational complexity is analysed. The results of the analysis of the computational experiment are presented where the program realizes the method on model problems of various types illustrating the effectiveness of the method. Consider the Cauchy problem for the system of m differential equations A(z)Y'(I)=B(z)Y(z)+D(z), Y(O)=Yo,
(1)
where I=[O,H],A(Z)and E(m) are matrix functions, and Y(s) and D(z) are, respectively, unknown and given vector functions. Suppose that in the interval of integration the matrix A(X) is continuously differentiable and the matrix B(x) and the function D(x) are continuous. Then this system is equivalent to the system of integral Volterra equations of the third kind
qyw, where
~l-w+,
(2)
101
R[Y(z),+]=A(s)Y(t)-
j
C(t)-d?(r)+A'(z),
c(t)Y(t)dt,
GM)== U(L)dt+A(O)Yp. P
In system (2) we approximate the integration operator i * operator 1% of numerical integration (see [ll):
in the interval [O, HI by the
(3) where L
-
H
z
R
%f(q):
=LIf(~),+tl,
=,=(I-Bf)H/2,
zi=eos(l~/n),
n
i-0,i ,...,
n,
1-o
are Chebyshev nodes of the second kind in the interval [O, HI, ” I, a~=“&?+a,f . . . + a._l+aJ2. c f-e = of the operator In are defined as the integrals The coefficients of the matrix {arj)o of fundamental interpolation polynomials on Chebyshev nodes ZI “‘j==--
-
“,s
4 (0 6
1,t-0
* 2*....4
(4)
*
where, according to [21, p.94,
and Tk(t) are Chebyshev polynomials. Explicit formulas for the coefficients are obtained in 131. After approximating the integration operator in (2) by the operator In we obtain a system of linear algebraic equations of order n-m = N of the following form:
(5) where
D,~n(rr)-(aI(zi),...,
d&c)); &-A(q), C&'(+,) are square matrices of order rrt, i = O,l,. . ..n. *
m
)I
The unknown variable in this system is the net vector function
t-i, 2,...,n, Y<=~Ylt,...,!flnc~, approximating the required solution Y(x) at the nodes of the net (XS)P.In this way we obtain a method of discretizing 14, Chapter 71 the Cauchy problem (1). 2. The approximation of the original system of linear differential equations (1) by the system of linear algebraic equations (5) can be estimated as the residue in the approximation of the integral operator P[U(zf,
zl:=K[U(zf,
SW(Z)
(6)
(here the operator K and the function G are defined in (2)) by the operator on net functions with Chebyshev nodes (3)
102
Fd (U,ll”, z,]=&[
(uj)~“,
(7)
a-G,
are defined where Uj-= U(zj), j = 1,2,...,n and the operator K, and the net function {G,),” The residue of approximation of the operator F by the operator Fn on the function in (5). 'J(a)@(F) is taken to be equal to the net function
i=l, 2,..., R. rin(~(z)l=PfU(~),rcl-Fn~(U(~~)}t", al, Lem
I.
The residue of approximation of the operator F by the operator F,
(8) on the
function V(x) is equal to the residue of approximation of the integration operator i 0 the numerical integration operator In on the function
by
R(z):=R[U(z)]=(A’(.r)+B(z))U(z)fD(z),
namely, we have the equality
7
rr"(U(s)l=-P,,[ ,R(r)].
i=i,Z,...,n.
0
Proof.
By the definition of the operators F in (6) and Fn in (7) we have
CorotZary 2. If in the Cauchy problem (11 the matrix A(x) is continuously differentiable, the matrix B(x) and function D(x) are continuous, then the condition of approximation (4, p.4451, of the operator F by the operators Fn is satisfied in the space of continuous functions C[O, H] and the method of discretization (5) is unsaturated. The residue of the Chebyshev integration process on a continuous function Proof. f(x) does not exceed
where EnIf] is a quantity of the best uniform approximation of the function f(x) by algebraic polynomials of order n, and the uniform norm of the integration operator in uniformly bounded in n (see (11): IlMicp, x]=H[1+0fn-‘Inn)].
Consequently, the order of decay of the residue of the process of discretization (5) increases together as the regularity of the solution Y(r) and functions A(r), B(r), D(X) increases, i.e., this process is unsaturated. 3. The stability of the discretizationmethod (5), unlike that of finite-difference methods, can be investigatedmore naturally in terms of general theory of abstract approximation schemes (4, p.2971. An investigationof the stability of the method in the case of the Cauchy problem for ordinary differential equations solved with respect to derivatives is carried out in (51 where the approximation of the integration operator by the polynomial analogue (11 of the operator of numerical integration (3) is considered. The results obtained in that paper can be extended in a natural way to the Cauchy problem for the system of differential equations (1) solved for the derivatives (A(Z)aE={8ij}tlli). The solvability condition for system (5) in this particular case (by the method of simple iteration) takes the form
where
(((fi),m((, is the length of the vector in the space Em. When this condition is satisfied it is obvious that the following conclusions hold: 1) the matrix K, of the system (1) is continuously invertible in the space of net vector-functions{F(z~)},* with the norm
103
2)
the norm of the inverse matrix is uniformly bounded in n,nz=no: ilK~-'iLc(l-q)-',
and, consequently, the stability condition is satisfied [4, p.2971. Suppose that the following conditions are satisfied for system (1) in the Theorem 1. interval [O, HI : 1) the matrix A(X) is continuously invertible, 2) the matrices A-~(r),A’(+) and B(z) are continuous. Suppose that the length of the integration interval does not exceed the quantity q
Proof. We consider the scalar case (m = 1). By the condition of the theorem, the function kl(s) exists and is continuous. Therefore, in this case system (5) is equivalent to the following: Yp.4
(91
-‘(M
The solvability of system (9) by the method of successive approximations obviously occurs in the case when the following inequality is satisfied: IIA-‘(~)licllz”llcllE(~)llc=q,.;;q<1. By the condition of the theorem this inequality occurs if
According to the estimate of the norm of the numerical integration operator [l] the last inequality holds for asno( Moreover, the sequence of norms of matrices inverse to the matrices of system (8) is uniformly bounded:
which means that the sequence of norms of the matrices inverse to the matrices of system (5) is also bounded:
The theorem in the vector case (m t 2) can be proved analogously. 4. The convergence of the sequence of solutions {(Y&n)~zi of the discrete problem (5) to the solution Y(Z) of the initial Cauchy problem (1) follows from the approximation and stability [4, p.2981. Moreover, the order of the error of the method
is equal to the order of approximation of the integral equation (2) by system (51, i.e., to the order of the residue (~P[Y])~,, of the approximation (8) of the operator F from (6) by the operator Fn from (7) on the exact solution Y(X) of Cauchy problem (1). Theorem 2. If the solution of Cauchy problem (1) exists and is unique and the operator X, in system (5) is continuously invertible, then the following equalities hold: (e~")l"=Kn-'(r,nIY(z)l),",
(10)
Ile”ll=C,llr”ll,
ill)
(IIK*II)-‘c;C,r:IIR*.-‘II.
Proof. We consider the integral equation (2) in the nodes of the net {z,),". Subtracting from each row of the system of equations obtained the corresponding row of system (5) we
104
obtain the following system of equations i31,2,.... K[Y(r), zr]-K,[(Yj)I", XiI==G(~I)--Gil Adding and subtracting from the left-hand side of the equations of this system the zi]and transferring the difference MY(r),+~I-&I(Y(~J)}I”, 4 over to the quantities K,[{Y(zj))P, right-hand side we obtain the equivalent system i! x~[{e"(r,)),,,t,l=r"[~ ,R[Y(z)l), i=1,2,... .
(10')
0 By the condition of the theorem, the operator Kn (the matrix of system (5)) is continuously invertible. Consequently, system (10') can be solved and its solution can be written in the form of Eq. (10). The estimate of the operator parts of systems (10) and (10') in the norm /I*[ I,, of the space RNcompatible with C[O, HI is a uniform norm of m-dimensional space of functions and gives Eq. (11). CorotZary. If the sequence of norms of the inverse matrices of system (5)
where E,[l?[Y]Jis the quantity of the best approximation of the function R[Y(xll by algebraic polynomials of order n, then the Chebyshev method (5) of the discretization of the Cauchy problem (1) is convergent. 5. The program realization of the discretizationmethod (5), (4) and (2) was executed as a consecutive realization of the following algorithms: 1) computation of the coefficients of the matrix of the numerical integration operator (3); 2) computation of the coefficients of system (5); 3) the solution of system (5) by Gauss's method without choosing the principal element. The computations by these algorithms requires that ~(nr)+O(nnz)+O((mn)s)=O((nan)g)is satisfied for arithmetic operations and o(n)+O(mQ~)=O(m*n) for the operation of computing functions. We will estimate the number of operations required by the method for the class of Cauchy problems of the form (1) satisfying the conditions of Theorem 2 with matrices A(x) and B(m), the function F(x), and the solutions belonging to the class of analytic functions in the ellipse of the complex plane with foci in the ends of the integration interval. According to Theorem 1 the number of net nodes depends on the required accuracy in the following way: n=O(lne) Hence, the dependence of the number of operations on the accuracy of the solution is O((rnhl8)~)for arithmetic operations and O(mslne)for the operations of computing functions. In the case of classical numerical methods of order p the dependence of the number of net nodes (and the number of arithmetic operations and operations of computing functions) on the accuracy has higher order o(mt+'p). 6. A numerical experiment to verify the effectiveness obtained in the assertions of Sections 2-5 was performed on a computer. The problems presented in [S-l03 and a whole series of equationsof the form(l) were used as model ones. An analysis of the results of the numerical experiment completely confirmed the effectiveness of the assertions obtained in Sections 2-5. The results of the analysis for each of the classes of differential equations is presented below. 7. The analysis of the results of a numerical experiment for classical nonstiff Cauchy problems with a unique solution belonging to the class of analytic functions 1 7 3 5 ii enabled us to draw the following conclusions. FIG. 1. The error of the method for the Cauchy problem is en and the residue of the approximation of integration operator (2) on the SXSCt SOlUtiOtl Of this problem is Pn and decreases as o($") when the number of nodes
105
increases. The decrease continues up to some value n,=no(e ) determined by the camp characteristicsof the computer (the length of computer words). Any further increase in the number of nodes n>no does not in practice lead to any change in the error of the method and the residue of the approximation. does not in practice depend on n. 2. For a number of nodes n
Y
(0) = I, y(O)=l,
ZE
(0,IOm],
.z=[O,
lam], m=O, i.
(12) (13)
Graph 1 shows the error of the method for problem (12) when H = 10; graph 2 shows the residue rn of integration of the right-hand side of Eq. (12) for H = 10; graph 3 shows the error of the method in problem (13) for H = 10; graph 4 shows the residue rn of integration of the right-hand side of (13) for H = 10;
graph 5 shows the error of the method and
the residue of integration of the right-hand side of problem (12) for H = 1 (in this case graph 6 shows the error of the method and the residue of the integration of the lle"ll/llr"ll~i); right-hand side of problem (13) and the model problem of [91 for H = 1. In the same figure (graph 7) we show the results of the solution of the Cauchy problem of type (13) in the interval LO, 11 by the Adams-Bashforthmethod presented in 191. A comparison of graph 7 with graph 6 well illustrates the above conclusions on the higher numerical stability of the proposed method. Using graphs 6 and 7 it is easy to compare the economy of the methods. Thus, to find the solution of the problem [91 with error Il~,li
the decrease of the ratio
lie,ll/llr,ll
as the stiffness of the system increases
up to a certain value lIBIIHza-l; b) the error of the met~~~pis independent of the length of the integration interval for II>H,,where ~~ze-1 camp /llBll. These conclusions are well illustrated by the results presented below on the solution of the following model problems: from [7, 101, y’=-l@y+f(s),
y(O)=&
(14)
where f(m) = (lOk - l)exp(-m), k = 0,...,12; from [a], I-2exp(-2) .zexp(-2) Y'(z)=lO‘ Y(s)+F(z), Y(O)= I II -sesp(-s) I-2exp(-s) II1' II
II
(15)
where F(x) is chosen so that Y(S) = (exp(z), exp(-z)), and also the one-dimensionala alogue of this system. These model problems are solved in the interval [O, HI, where H = 10't , 1 = -2,...,12. The results show that if for problem (14) we have H < 100, then Ilenll/llr,ll=O(lOLH)-l,n~lO, and if H 2 100, then lle.ll=lOW. For problem (15) and its one-dimensionalanalogue if H 2 100, then IleJ=l. 9. An analysis of the results of a numerical experiment on well-posed Cauchy problems for systems of the form (l), unsolved for the derivative, in particular on the problem
106
showed that such systems can be solved as well by the Chebyshev method the systems solved for the derivatives. 10. An analysis of the results of numerical experiment on well-posed Cauchy problems for systems of the form (1) with singularitiesof nodal type, in particular on the problem 2zy’=ky-f(r),
=P,
11,
v(O)=k-‘f(O),
(16)
where k = 0, 1. 1.5, 2-10-p, 109, p = 1,...,12, q = 1 ,...,4, fir) = (2% + kfexp(-zz),and the initial point is singular, enabled us to draw the following conclusions: 1) the solution by the Chebyshev method is unique; 2) this solution converges, as a rule, to the most regular of exact solutions of the Cauchy problem; 3) the ratio ~~~~/~r~~ essentially depends on which type the singular point belongs to. In particular, the error of the solution of problem (16) by the Chebyshev method depends on the parameter k in the following way: if k increases from 0 to 2, then the errors lieloll increase from 10-la to 10, and if k increases from iO*to IO’, the errors lleloll decrease from IO-' to fO+O. REFERENCES DENISENKO P.N., Chebyshev integration processes, in: Proc. Int. Conf. on Approximation Theory, 131-135, Nauka, Moscow, 1987. 2. PASHKOVSKII S., Numerical Applications of Chebyshev Polynomials and Series, Nauka, Moscow, 1983. DEYAD~ V.K., An approximation-iterativemethod of approximating solutions of 3. Cauchy's problem for ordinary differential equations, Preprint 84-27, Inst. Matem. Akad. Nauk UkrSSR, Kiev, 1984. 4. TRENOGIN V.A., Functional Analysis, Nauka, Moscow, 1980. 5, DENISENKQ P.N., To the solution of the Cauchy problem by interpolationon quasiChebyshev nodes, Differents.Uravn., 5, 5, 908-913, 1979. HAMMING R.V., Numerical Methods, Nauka, Moscow, 1968. 6. 7. RAKITSKII YU.V., USTINOV S.M. and CHERNORUTSKII I.G., Numerical Methods of Solving Stiff Systems, Nauka, Moscow, 1979. 8. HALL G. and WATT J.M. (Eds.)., Modern Numerical Methods for Ordinary Differential Equations, Clarendon Press, Oxford, 1976. BABUSKA I., PRAGER M. and VITkEK E., Numerical Processes in Differential Equations, 9. Wiley, London, 1966. 10. GEAR C.W., Numerical Initial Value Problems in Ordinary Differential Equations, Prentice Hall, Englewood Cliffs, 1971. 1.
Translated by M.C.