EL-SAFTY, A.

English

Periodica Polytechnica (Budapest University of Technology and Economics)

ON THE NUMERICAL SOLUTION OF RETARDED ORDINARY DIFFERENTIAL EQUATIONS BY A CLASS OF FORMULAE BASED ON SPLINE FUNCTION APPROXIMATION A. EL-SAFTY Department of Mathematics of the Faculty of Electrical Engineering. Technical University H-1521 Budapest Received November 15. 1982 Presented by Dr. P. RozsA Summary This paper presents a class of spline approximation method for the numerical solution of retarded ordinary differential equations' (RODE). The proposed method is explicit, for convergence we need several restrictions and assumptions on the right side of the RODE. For other application of spline functions see [I], [2J, [3]. A scheme for constructing mt;thods of arbitrary ordel and numerical examples is given. Introduction In recent years there has been a growing interest for the numerical treatment of delay differential equations. This is due to the versatility of such equations in the mathematical modeling of processes in various application fields, where they provide the best and sometimes the only realistic simulation of observed phenomena. In difference-differential equations (DDF), or more generally, in functional differential equations (FDE) the past exerts its influence in a significant manner upon the future. For general treatment of FDE and for their history we could refer to the book of Hale [4]. Nevertheless the numerical methods which deal with this kind of equation were not developed so far. Cryer and Tavernini [5J studied the Cauchy problem for Volterra functional differential equations: y'(t) = F (y, t), t E [a, bJ (1.1) y(t) = g(t), t E [ex, aJ Here, F:C([e<,bJ-En)X[a,bJ-En is a Volterra functional, that is, F(y,t) depends on t and on y(s) for s> t; and the function 9 E C([ ex, a] - En) is specified initial function. The problem (1.1) includes as special case the initial value problems for ordinary differential equations (ODE), retarded ordinary differential equations (RODE), and Volterra integro-differential equations (VIDE). 124 SAfTY, A, EL Numerical method for solving special cases of (1.1) has been considered by several authors, Feldstein and Goodman [6], Bleyer [7], Bleyer and Preuss [8], Tavernini [9], Kemper [10], Oberle and Pecsh [11]. Description of the method We consider a RODE in the form y'(x) = J(x, ytx), y(a(x))), x E [a, b] with y(a)= Ya and a(a)=a (2.1 ) (2.2) where a(x) ~ x, y(x) is an unknown function. The function a(x) is usually called the "retardation" or "lag" function. Generalizations to multiple lags are carried on in the obvious way: Let us denote djy/dx j by yU); if yU) is continuous at Xi(Xi E [a, b]), we denote it by yV). The notation yV) will be for the left-hand limit of yU) at Xi, when y(j) is discontinuous at Xi and sometimes, for convenience, even when y(j) is continuous there. The solution of th~ problem (2.1) with (2.2) is denoted by v(x). J(j) stands for the higher derivatives of y'(x}. Choosing h* > 0, define a~ = a + 11* = xo, then we define the algorithm on [a, a~] by using Taylor's expansion as follows: Let us start the algorithm by VIOl = v(a) = V .. a " ~ a "Ul = JU - 1 )(a VIOl ,;(0»). 1 ::;_J'::;_ m )la ' ... a ')':XQ, where m is fixed integer. For X E (a, a~], y(x) is given by formula m 1 y(x) = I -k' y~k) [x-aY k=O . we define for O~j~p m h*k - j y~) = I, (k- _ ')' y~k) k =} J . where O~p~m-l, and for p+ 1 ~j~m VU) = [U - 1) (X VIOl v(O») _0. O"O'.:xo where y~~) can be calculated from (2.5) m 1 y(a(xo)) = kf:o k! y~k) [a(xo)-aY (2.3) (2.4) (2.5) (2.6) (2.7) (2.8) NUMERICAL SOLUTION Of' ORDINARY DIFFERENTIAL EQUATIONS 125 and the higher derivatives yU>(Ct:(xo)) = f, (k~ ')' y~k) [:x(xo)-ar- j k =) } • (2.9) we continue the algorithm in term of "mesh points" Xn and "non mesh points" Ct:(xn) defined as follows. Let us Pick h > 0 and integer N so that N h = b - ad, denote xn = ad + nh for n=O, 1, ... , N. Let us define a'! = max {xj/xj > a'!-I; :x(xd;£ a'! I for k;£j;£ N} (2.10) It may happen that the procedure defining the numbers a,!, which are always mesh points, terminates before a'! = x N = b. In this case let the last defined a'! be denoted by x j and by definition let a'! + I = X j + I . Since Ct:(x j + I) ~ a'! we could say that in the interval (x j' x j + I) we have no delay and :x(x) ::::: X can be assumed. By assumption :x(x);£x on [a, bJ, the sequence at is finite, the possible number of a'! -s is at most N. The algorithm now will be defined first on (ad, ai], then on (at, an and so on. If Xn E (ad, ai] then by (2.10) x: = :X(Xn);£ ad. Let Xi E (ad, ai] and define for O;£j;£p m hk - j I,(j) _ '\' __ \!(k) Ji - L., (k- ')' J(i 1) k=) } • where 0;£ p;£ m-I is fixed. If p + 1 ;£j;£ m then by (2.5) and (2.7) v(jJ = j(j- I)(X" v(O) v(O») ..- z Z,-,l '", ,:);1 is well - defined. For x E (ad, ai], y(x) is given by formula lI(x) - f 1 I,(k) [x - X Jk -'. - kf-O k! J(i-I) • i-I where X i - I <X;£Xi' (2.11 ) (2.12) (2.13) In the next step we move to the interval (af, a:rJ, a!;£ b. If x/ E (af, an then a(x/) E [a, ai], if :x(x/) E (ad, ai], then we need the approximation of v(Ct:(x/)). y<~) will be counted from (2.13) :x (2.14) whenever Xi- 1 < a(x/);£ Xi' In the case a(x/) E (a, ad] then y~?) is given by (1.5) at x = Ct:(x/). Higher derivatives can be get similarly fromj-th derivative of (2.13) at x = Ct:(xtJ or respectively from the derivatives of (1.5) evaluated at x = Ct:(x/). vU)(a(x))=vU) = f _l_y<k)[Ct:(x)_aJ k - j • -/ ':Xl k=o(k-j)! a / (2.15) 126 SAFTY, A. EL. and (j)(cx(x »_yUl - f _1_ (j) [() Jk- j Y I - '" - ~ (k-j)! YU-1) CXX, -Xi-1 (2.16) for Xi _ 1 < cx(x/) ~ Xi and for 0 ~j ~ m. Hence similarly to (2.11) and (2.12) for O~j~p m hk - j (j) _ " (k) YI - Lo. (-k ')1 Y,-1 k=) -j. where for the smallest 1, X , - 1 =af, y(j)= jU-1)(Xh y!O), y~~») (2.17) (2.18) when p+ 1 ~j~m. By the assumption on cx(x) the algorithm terminates in at most N steps, if we repeat the above method for (a!, an, .... Obviously (2.11) is thej-th derivative of(2.13) evaluated at Xi' Therefore y(x) E CP[a, bJ while the (p+ l)-st through m-th derivatives of y(x) generally have finite jumps at the mesh points, the values given by (2.12) (or for the next interval (aT, a!], by (2.18) being left-hand limits of these derivatives of y(x). . Convergence We begin this section with five lemmas which will prove useful in the following convergence (for to proves see [12]). Lemma 1 Hj, p and r are non-negative integers and hand ri , are real numbers, then f f (rh)k- j h' - k ril = f «r+ 1)/~)k-j CXk. k=jl=k (k-j)! (l-k)! k=j (k-J)! Lemma 2 Let Bj,i' l1j,i and f3 j ,i be real numbers which satisfy ~ hk - j m hk - j h m - j + 1 Bj • i = /;;j (k-j)! Bk.i-1 + k=~+l (k-j)! l1k,i-1 - (m-j+ I)! f3 j ,i-I for O~~p, 1 ~i~N. Then for any integer q such that 1 ~q~i ~(qh)k-j ~[~ m «r-l)h)l<-j hI-I< ] Bj.i = 4:-. (k-j)1 Bk.i-q + ~ 4:-. _L (k-j)1 (1_k)I 111•i - r -k-J • r-1 k-J I-p+l • • ~ [ ~ «r-l)ht- j hm-k+ 1 ] - r~l /:;!j (k-j)! (m-k+ 1)! f3k,i-r for O~j~p, 1 ~i~n NUMERICAL SOLUTION Of'ORDINARY DIFFERENTIAL EQUATIONS 127 Lemma 3 Let j, p, I and r be non-negative integers, and let h ~ 0 be real. If 1-;£ p-;£j and r~ 1, then Lemma 4 P «r-l)h)k- j h1- k ,. (2h)l- j _' L --s--rP J k=j (k-j)! (l-k)! - (l-j)! Let q ~ 0 and r ~ 1 be integers. Then i L r q -;£ iq + 1 r= 1 Lemma 5 Let (Xi be real numbers and C and D be non-negative real numbers satisfying the recurrent inequality i-I I ctil -;£ D L I (Xk I + c k=O Then the following bound holds l(Xd-;£C eiD, i~O Let us consider the errors on the interval (a, a6J (a6 =xo). G - y(k) _ V(k) k.O- 0 0 , n - V(k) V(k) 'lk.O-.O - 0' O-;£k-;£p p+l-;£k~m Between (a, a~) we obtain by Taylor expansion m h*k - j I *m - j + 1 (j)_" (k) 1 (m+1)(V) Vo - /;;j(k-j)! Da + (m-j+l)! V t,j,a where a<~j.a<ao. From (2.6) and (3.1) it follows that f h*k - j m kk - j h*m - j + 1 _ + " (m+l)(" ) G}.o - k=j (k-j)! Gk.a k=7;+1 (k-j)! '1k.a(m_j+ I)! V C;j.a for O~j~p. Since Gk,a=O and '1k.a=O, we get h*m- j +1 _ (m+ 1)(;: ) Gj.o- - (m-j+ I)! v ~j.a Assuming (3.1) (3.2) (3.3) (3.4) (3.5) 128 then SAFfY. A. EL. sup Iv(m+ 1)1 ~ Vm+ 1 [a.b] h*m- j+ 1 IG'ol< V.·1 ). = (m-j+ I)! m-r (3.6) (3.7) For the continuous error bounds, by Taylor expansion of y(x) and v(x) on the right-side end of (a, a6'J we can obtain for 0 ~j ~ m h*m- j +1 sup lyU) - vU) I < V. [a,ali] (x) (x) = (m _ j + I)! m+ 1 Assuming 0 ~j ~ p and set hP - j+ 1 bP1o) = (m~j+ I)! Vm+1 sup lyU) -lP) 1< b rh )h*m- P (x) (x) = )1. 0 [a.aO] by taking h* ~ho. For the case p+ 1 ~j~m, we have b rJ Vm+ 1 )1.10) = (m-j+ I)! Then for p+ 1 ~j~m sup I VU) - vU) I < b .(h )h*m- j+ 1 _ (x) (x) = ) 0 [a.aol Now we have given the error bound for 1Jj.o as follows Y/j.o = fU - 1 )(xo, Vo + Go. 0' V:>o + C:>. obo(Jlo)h*m - P) - fU - 1)(XO' Vo, v:>o) for p + 1 ~j ~ m, where 0 ~ Ic:>. 0 I ~ 1. Therefore by the virtue of the mean value theorem: C . 1Jj,O = Co, 0 ay pr 1)(XO, Pj, 0 Vo + (1 - Pj, o)y~o>, y;~)+ (3.8) (3.9) (3.10) (3.11) + c b (h )h*m-p~fU-l)(X v· " ,. +(1 " )y(O» :>0.00 0 62 0' O,rj,o":>o -'j,O:%o where 0 < Pj, 0 and r'j.o < 1. NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS Using the following notation: rn = a jU-l)(x R ° v +(1- R ° )y(O) y(O) }.O ay 0' p}.O 0 p).o 0, 20) F<2>- G jU-l)(X v y. v +(1 ') )11 ) j.O - GZ 0' 0' j.O 20 - {j.O Jao we can write - FO) ;;: (h )h*m- P F<2) t7j.O-CO.O j.O +C:ZOUO 0 j.O Assuming then by (3.13) I ;lj.o I ~ Ico. 0 I Mj~>1 + 60 (h o)h*m- P Mj:>1 set Mj=max {Mf>, M?>}. Then for p+ 1 ~j~m lt7j.ol ~ [I co.ol + 60(h o)h*m - PJ M j- 1 From (3.7) and (3.9) we have O~j~p and p+ 1 ~j~m Now let us consider the errors on the mesh points of (a6, aiJ O~k~p p+l~k;£m Between mesh points (Xi-I' Xi) we obtain by Taylor expansion m hk - j hm - j + 1 .(j)_ \' ,(k) .m+l(- ) Vi - 4:!o(k_ O ), Vi-l+( _'+l)I V <;j.i k-) J. m J . where Xi-l < ~j.i-l <Xi' From (2.11) and (3.19) it follows that hk - j m hk - j Cj.i= fo(k_ ')1 ck.i-l + I (k- ')1 t7k,i-l-k=) J. k=p+l J. hm- j+ 1 m + 1( y ) (m_j+l)!V <;j,i-l 129 (3.12) (3.13) (3.14) (3.15) (3.16) (3.17) (3.18) (3.19) (3.20) (3.21) (3.22) 130 SAFTY, A, EL. By applying Lemma 2 and by using (3.5) one can obtain _ f (iht- j hm-k+l (m+l)(Y) Gj,j- - /~j(k-j)! (m-k+ I)! V C;k,a + i [P m «r-l)h)k- j Iz' -k ] + rf:t k~jl=~+l (k-j)! (l-k)! ~"i-r -" " .(m + t)( y ) i [P «r- t)h)k- j hm - k+ t ] - r~l l~j (k-j)! (m-k+ I)! r; C;k.i-r (3,23) for O;£j;£ p, 1;£ i;£ N, To obtain some estimation for Gj.i we must have some bounds for 1Jj,i j j-1( )_ a jU-1)( P . (I-P ) (0) (0» - Xi' Vi' V". - Go i -::;- Xi' J' iVi + J' i Yi • y". + I , ay , , . + ;: (Iz )Izm- p °jU-1)( , +(1-') (0) G".uo 0 '" Xi. Vi' '1j iVa.' '1j i y". I CZ . I .,where 0 <pj,i' }'j,i < 1, 0;£ 1 GO!; 1;£ 1 with notation similar to (3.12) and (3.14), and by applying Lemma 3, we obtain for O;£j;£ p and (3.25) i [ m (2h)' - j ] ' 0;' *m-p (2) __ p-j + L L (lGO.i-rI M ,-l+(jo(lzo)1z M ,- 1) (1_')1 r r=l I=p+l j, Denoting 11 = {i/l ;£ i;£ N. Xi E (a~, aTJ} then by applying Lemma 4 we get (2h)m- j+ 1 , (2h)m- j+ 1 , ! I<V. (' l)P-J V. 'p-J+l Gj.i = m+l (m-j+ 1)! l+ + m+l (m-j+ I)! I' + m (2h)l- j 1 ! 'P- j+ 1 " MO) + max GO,i ! L... -(I ')1 1-1 + 1\ I=p+l -j. m (2h)l- j +b (h )h*m- PjP-j+l '" __ M<2~ o 0 L... (I ')1 1-J • I=p+l -j. (3.26) From (3.25) takingj=O, since rP;£iP, we obtain NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS 131 I 1< ~ 'p ~ MO) (2h)1 eO.i = L... eO.i-rl L... 1-1 -t'- + r=l I=p+l • V. (2h)m+1 2P 'p V. 'p+l (2hr+1 + + m+l(m+l)! 1+ m+l' (m+l)! + (jo(ho)h*m-PiP+ 1 [ f Mi:~ (2t7l ]. (3.27) I=p+l • In order to write (3.26) and (3.27) in a denser form we introduce the notation m 21- jhl-p-l A -(b )P-j ~ M<l) 0 j- -a L... 1-1 (t ')' I=p+l -J. 2m - j + 1 . 2m - j + 1 BJ. = (b - alP - j + 1 V. + (b a \P- ) V. 2P + m+l m-j+ 1 -} m+l (m-j+ I)! +(b-a)P-j+l(j (h) ~ M<2) 21-j hl - p - 1 . 0 0 _L... I-I (/-J')' 0 I-p+ 1 . Now taking h;£ho and Ji=max {h*,h}, '<N- b-ao ,= --h-i-I leO.il;£AoJi L leo.rl+BoJim- p r=O by applying Lemma 5, (3.30) implies le .I~B· Jim-PeiiiAo~B e(b-ao)AoJim- p 0,1 - 0 _ 0 therefore (3.28) (3.29) (3.30) (3.31) lej,il;£ [Boe{b-aoAo Aj(b -a) + Bj]Jim- p =O(Jim- P) (3.32) for O;£j;£p, i Ell' and then l1Jj,d;£leo.dMJ~>l + (jo(ho)Jim-PMJ=>1 . set Mj=max {MJI>, MJ2>} , Now we give the continuous error bounds on (ao, aT], By Taylor expansion of y(x) and v(x) on the right end of each subinterval (Xi' Xi + 1] of (a~, aT] we can obtain for O;£j;£ m 132 S.4FTY. A. EL. max {SUp {lyW)-vl~)I: [Xi' Xi + I]}: 1 ~ i~ 11 -1} ~ ~ It - j m It - j Jim - j + 1 ~ L.(k_')If!1aX 1ek •d+ L. (k_')If!1ax l'1k.iI+ ( -'+1)1 Vm + 1 k:J J. ,el! k:p-r 1 J, ,ell m ) ' Assuming that O~j~p and set hk - j Mj(ho)= f'(k~ ')1 [Bor(b-ao)AOAk(b-a)+Bk]+ one gets for 0 ~j ~ p k:J } • m hk - j + L ('~')IMk_1[Boe(b-ao)AO+bo(ho)]+ k:p+ 1 X J, h'O- j+ 1 + (m-j+ 1) Vm + 1 sup lyUJ(x) - v(j)(X) I ~ M /ho)Jim - p [ao.ali] (3.34) (3,35) (3,36) For p + 1 ~j ~ m, M j(ho) is different from (3.35) since in this case the last term of (3,34) gives the essential part to determine the error bounds: m hk - p + 1 V (Miho)= k"5;;j (~-j)! Mk_l[Boe(b-a)Ao+bo(ho)+ (m~~ 1)! (3,37) Then for p+ 1 ~j~m max {sup {Iy(j)(x)- v(j)(x)I : [Xi' x i+ I]}: 1 i 11 :..-1} ~ Miho)Jim - j+ 1 (3,38) Now we are able to give error bounds whenever our procedure is continued to next steps (a l' a2), , " and it is obvious that the same error bounds could be obtained as (3,32), (3.33) (3.36) and (3,38), For given El' E~ > 0, ho can be chosen such that bo(ho)h'o- p and M o(ho)h'O-P<E(E=maxE1 ,E2), This implies that y(x)-tv(x) in Ghebyshev norm on [a, b]. Since Aj and Bj are independent of Ji, bo(ho) and Mj(ho) are bounded, Hence bo(ho)Ji'O:8 and Mo(ho)h'O-P-tO as ho-tO, We can summarize our results, Theorem Let the function J(x, y, z) be k-times continuously differentiable with respect to X, y, z in same domain D of (x, y, z)-space (k ~ 1), Let us assume that tX(x) ECk [a, b]. The restriction of D into (x, y) - space will be denoted by D1 , Let us assume that Dl contains the exact solution v(x) of the RODE on [a, bJ, so that v(x) E Ck + l[a, b]. Let 1 ~m~k and O~p~m-l, Then there is ho >0 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS 133 such that the algorithm producing the spline function y(x) E PCm • p [a, b], satisfies the discrete error bounds max lij)(Xi)-VU)I~LjJim-p O;iii;iiN Moreover, y(x) also satisfies the Chebyshev error bounds sup lyU)(x)-vU)(x)1 ~Liho)Jim-p XE[a.a('H sup ly(j)(x)-v(j)(x)I~Mj(ho)Jim-p, [a~.bl O~j~p sup lij)(x)- v(j)(x) I ~ Liho)Jim- j+ 1 (a.a~l max sup lyU)(x)-vU)(x)I~Mj(ho)Jim-j+l, 1 ;iii;iiN (Xi-I.Xil p+ 1 ~j~m where L j , 6j (ho), Mj(ho) are given constants. Numerical Examples (3.39) (3.40) (3.41 ) To illustrate convergence, we applied the method to the following examples. Example 1 y'(x) = 2y(.jx) y(l)= 1 1 ~x~2 which has the exact solution Example 2 y'(x)= 1 2~' ) y(O)=O -y smx which has the exact solution y=arcsin x. re O<x<-= = 2 The errors for h=0.05, m =3, p=2 are tabulated in Table 1, for example 1 and in Table 2, for example 2. 4 Per. Pol. El. 27/2 x 1.1 1.2 1..3 L4 1.5 x 0.1 0.3 0,5 0.7 0.9 1.21 1.4·4 L69 1.96 2.25 SAFTY, A. EL. Table 1 h=0,05 m=3 p=2 exact solution 2.1 2.4 1.6 2,8 3.0 2,0 0,0 2,0 0,0 2,0 0.0 2.0 0.0 2,0 0.0 Table 2 h=0,05 m=3 p=2 1.21 '2.2 1.44 2.4 1.69 2.6 1.96 2.8 2.25 3,0 appr. soL error -~- -----~~.- .~- -~ - ---_._-0,10016800 0.30479035 0.52425301 0.77820489 1.11976955 0.00000058 - 0,000 097 00 - 0.000 654 23 -0.00287387 0.00568089 -0,00042433 -0.00119856 -0.00517953 -0.01769241 0.29033550 References -0,00683370 - 0.009351 49 -0,03243105 -0,04298330 7.03998204 2,0 -0,00003 2.0 -0,00002 2.0 0,0 2,0 0,00001 2,0 0.0 -0,06348575 -0.17004752 -0,67401591 -0,94177484 159.625309 1. KEMPER, G. S,: Spline function approximation, Siam, J. Num. Anal., 12, 73 (1975), 2, MDLTHEI, H, N.: Spline approximationen von beliebigen Defekt zur numerischen Lasung gewohnlicher Differentialgleichungen. Teil I, H, HI, Numer. Math" 32,147 (1979) 32,343 (1979) 34, 143 (1980). 3. LOSCALZO. F. R.-T ALBOT, T D.: "Spline function approximations for solution of ordinary differential equations." Bull. Amer. Math. Soc." 73, 438 (1967). 4. HALL, J.: Functional differential equations, New York, Springer Verlag 1971. 5. CRYER. C. W.-TA.VER"I"'I. L: "The numerical solution of Volterra functional differential equation by Euler"s method." Siam J. Numer. Anal. 9, No. 1, March (1972). 6. FELDSTEI!', A.-GOODMA"', R.: "Numerical solution of ordinary and retarded differential equations with discontinuous derivatives." Numer. Math. 21, 1 (1973). 7. BLEYER, A.: "On global methods to solve difference-differential equations with retarded arguments," Periodica Polytechnica. 22, 141 (1978). 8. BLEYER, A.-PREUSS, W.: "Spline function approximations for solution of differential equation with retarded arguments." Acta Math. Acad. Sci. Hung. 39 (1982). 9. TWER"INI, L.: "One-step methods for the numerical solution of Volterra functional differential equations. Siam J. Numer. Anal. 8, No. 4, December (1971). 10. KEMPER, G. A.: "Linear multistep methods for a class of functional differential equations." Numer. Math. 19, 361 (1972). I!. OBERLE, H. T-PESCH, H. T: "Numerical treatment of delay differential equations by Hermit interpolation." Numer. Math. 37, 235 (1981). 12. HULME, B. L.: "Piecewise polynomial Taylor methods for initial value problems." Numer. Math. 17, 367 (1971). A. EL-SAFTY H-1521 Budapest 

ON THE NUMERICAL SOLUTION OF RETARDED ORDINARY DIFFERENTIAL EQUATIONS BY A CLASS OF FORMULAE BASED ON SPLINE FUNCTION APPROXIMATION

https://pp.bme.hu/ee/article/download/4712/3817

ON THE NUMERICAL SOLUTION OF RETARDED ORDINARY DIFFERENTIAL EQUATIONS BY A CLASS OF FORMULAE BASED ON SPLINE FUNCTION APPROXIMATION

Abstract

Similar works

Full text

Available Versions

Periodica Polytechnica (Budapest University of Technology and Economics)