The DI methods for directly solving a system ofa general higher order ODEs are discussed. The

convergence of the constant stepsize and constant order formulation of the DI methods is proven first

before the convergence for the variable order and stepsize case

Suleiman, Mohamed

Universiti Putra Malaysia Institutional Repository

Pertanika 8(1), 59 - 66 (1985)Convergence of the Variable Orderand Variable Stepsize Direct Integration Methodsfor the Solution of the HigherOrder Ordinary Differential EquationsMOHAMED SULEIMANMathematics Department,Faculty of Science and Environmental Studies,Un:iversiti Pertanian Malaysia,Serdang, Selangor, Malaysia.Key words: Convergence; higher order ordinary differential equations (ODEs); Direct Integra-tion (DI) methods; predicted and corrected values; integration coefficients; variableorder and stepsize.ABSTRAKKaedah Pengkamilan Terus bagi penyelesaian terus satu sistem persamaan pembezaan biasaperingkat tinggi dibincangkan. Penumpuan bagi kaedah Pengkamilan Terus dengan langkah danperingkat yang tetap dibuktikan dahulu, sebelum penumpuan bagi kaedah tersebut dengan langkahdan peringkat yang berubah dibuktikan.ABSTRACTThe DI methods for directly solving a system ofa general higher order ODEs are discussed. Theconvergence of the constant stepsize and constant order formulation of the DI methods is proven firstbefore the convergencefor the variable order and stepsize case.1. INTRODUCTIONThe general system of higher order initialvalues ODEs is given by(d.) .Yj 1 = f j (x, Y), 1 = 1, 2, ... , S (1)with initial conditionsYea) = 1)in the interval a.,;;: x .,;;: b, where the i-th equa-tion is of order dj and(dl-l) (ds - 1)yT(x) =(Yl' ... , Yl , ;.. , Ys ' ... , Ys ),Key to author'. name: M.B. Suleiman.Many problems in engineering and sciencecan be formulated in terms of such a system.Some of the higher order ODEs problems foundin'the literature are on the symmetrical bendingof a laterally loaded circular plate, the bendingof a thin beam clamped at both ends. Both prob-lems are given in Russell and Shampine (1972).Others are on the steady flow of a viscoelasticfluid parallel to an infinite plane surface withuniform sunction; given in Serth (1975), controltheory in Enright et al. (1974), coplanar circulartwo body motion in Shampine and Gordon(1975), gravitational nobody problem in Dol etal. (1972).The current technique C)f solving the prob-lem in (1) is by reducing it to a system of firstM.B. SULEIMANorder equations. LetZll = yl' Z12 =Y~, ..., ZI d = Y11Integrating (3) d - r times, r == 0,1, ... , d1 leads to the identityd-l-r .hIy(r) (xn+ l ) = L - y(r+i) (xn ) + (4)i =0 i!Jxn + I JX ... fxx f(x, Y(x» dx ." dx .X n xn nand+- d - r times "'"*ZT = [Zll' Z12' .•• , ZI d , ..• , ZSI' Zs2"",Zsd ]I Sandf is then approximated with an interpolationpolynomial Pk n(x) of degree k - 1 which inter-polates fat the'k preceding points (fn, xn), (fn _ l ,xn _ I ), ... , (fn _ k + l , Xn _ k + 1 ). Thus expressedin terms of divided differences,then (1) reduces to the equivalent first orderproblem, Pk (x) = f[ ] + (x - x ),n n nZ' = F(x, z),with(2) f[ n, n- I J + ... + (x - x n ) ...z(a) = 77 (x - X n - k + 2 ) fIn, n-I, ... , n-k+1 J2 THE DI METHODS IN THEPREDICTOR-CORRECTOR MODE(5)hi . (r+i) JXn+IJx IX-Yn + '"i! x x xn . n nd-I-rLi = 0Hence a discrete replacement for (4) is(x - Xn_i+1 ) dx dx ... dx . (6)+- t times "'"*Define, the integration coefficients git tobe the t fold integralwhere p n~~ are predicted values associated withtotally know back values.f[n,n_I, ... , n-i](r)Pn+1wheref -f, .[n,n-I,.... , rl.-i+I J n-I, n-2, ... , n-I](3)y(d) = f(x, Y), Yea) = Awhere y T = (y, y', ... , Y (d-I»)and AT =(77,1]', ... , ned-I) ).The system in (2) is then solved using eitherthe Adams or the Runge-Kutta class of methods.But the work of Krogh (1969)' and Suleiman(1979) have shown that for many problems, it ismore efficient to solve (1) by integrating directly,rather than reducing it to (2). This DI method isin fact a generalization of the Adams-type ofmethods to higher order ODEs. Hall and Sulei-man (1r:,81) and Gear (1978) have studied thestability of the method for solving second orderODEs. Gear'(1971) proves the convergence forcertain methods which solved (1), those whichcan be transformed into the Nordsieck's arrays,while Shampine and Gordon (1975) prove thecase for the first order problem.For simplicity of description and withoutloss of ger>erality, the discussion which followswill be for the single equation60 PERTANIKA VOL. 8 NO. I, 1985CONVERGENCE OF THE Dl METHODSThen the predicted values in (5) can be 0 1 2 k + d + 1written asd-I-r hi Yn(r+i) 0per) L - + goo gOI g02 gO,k+d+1n+1 .,1- 1i = 0 glo gll gl2 gl,k+dk-I 2 g20 g21 g22 g2,k+d-1L gi,d-r f[n, n-l, "', n-i] , (7a)i = 0For the case i = 0,From the definition of gil in (6) it can beshown by integrating by parts thatg = (xn+1 - Xn-ifl)gi-I, t - tg i _ l , t+I' (8)i, tFormulae (7a, b) represent the k-step pre-dictor-corrector mode of the DI methods. In theva;;able order implementation of the method, k,if changed, is usually either decreasing or in-creasing by one from the value at the previousstep.hn + 1gl t = gl, t - 1 -- , t = 1, ... , k + d,t+13. CONVERGENCE OF THEPREDICTOR-CORRECTOR MODE OFTHE DI METHODSk + 1 gk+ I, PDespite the simplicity and advantages inspeed, storage and accuracy of the DI methods,they are yet to be popularised as a general pur-pose code for the solution of the general orderODEs.and hence the computation of the first row isavoided. Constant stepsize will further accountfor the saving of computation, since if the step-size remains constant for R steps, and the co-efficients of the R-th row depend only on thecurrent stepsize, therefore the coefficients thefirst R rows need not be recomputed and com-putation need only start at row R + 1. A changein order corresponds to adding or dropping aninterpolation point. If at the current point to becomputed, the step is not changed for R + 1steps and the order is reduced, no special actionis required. If however, the order is raised, oneadditional term needs to be added to the endtable for each of the first R rows of the inte:-gra-tion coefficients, before commencing with com-putation of the row R + 1.From (8) and (9) it is observed that gOt =gl I-I giving rise to the relation,(9)(7b)n-k+II'dx dx...dx = ht n+ It!r = 1, 2, 00" d(d- r) f*- w+- Pn+1 gk,r [n+l, n,(d-r)Yn+1The corrector formula is obtained by ap-proximating fin (3) by the k degree polynomialPk+l , n+1 (x) interpolating the (k + 1) points(f*0+1' x n+I), (fn, x n ), .. " (fn-k+.l' Xn_ k+I ),where f*n+1 == f(xn+l , P*n+I)' P*n =(d-I)(Pn, P~, ... , Pn ),In terms of Pk,n(x),Pk+l,n+l(x) =Pk,n(x) + (x - x n)(x - xn_ l ) ... (x - X n _ k + l )f*[n +I ,n '00" n_ k+ I I ' where the last divided diffe-rence is obtained using the f evaluated by the pre-dicted values p*n+I' Replacing Pk+l , n+1 (x) in(4) and integrating, we have the corrected valuesThese integration coefficients are generatedin a triangular array of the form shown belowstarting with the first row and from left to rightwith k + d + 2 coefficients along the first rowand (k + 2) along the first column,Before the proof of convergence of thevariable stepsize and order method in (7a, b) isgiven, the proof for the constant stepsize handorder k formulation of the method is first given.PERTANlKA VOL. 8 NO. I, 1985 61M,B, SULEIMANDefinitionA linear k-step method is said to be conver-gent if for any function f(x, Y), for which thesolution of the problem in (3) exists, and for allconstants 77' (,In),om = 0, 1, ... , d· - 1, lim Y~s) = y(s)(x),h--+Onh=x-as = 0, 1, ... , d - 1;polynomial interpolating the set of values{ (xn+l _ r, Yn+l-r)' ..., (Xn_ k+l , Yn-k+I)} .This form of the interpolation polynomial is usedhere for convenience of proving the convergenceof our theorem. The constant stepsize and orderformulation of the predictor in (7a) is(d-r) (d-r) (d-r+l) hr- I (d-I)'Pn+l =Yn +hYn + ... + -- Yn• (r-l)!k+hr ~ ri, d-r fn+l - i, r = 1, ... , d (10)i= Ip. = 0, 1, ...., k-l;where n = (x - a)/h for any x €. [a, b], withstarting valuesyCm) = 77(m )(h) satisfying the con-p. p.di . (m) (m)tlon lim 77p. = 77 0 'h-+Om = 0, 1, ... , d - 1.where the coefficients+- r times -+Li (x, 1) dx ... dxare constants independent of h. To see this, putx = sh + x n thenLet en be the global error-between the truesolution Y(xn) and the computed solution Yn'i.e.,:;; (y(xn) - Yn, Y (xn ) - Y~, ... ,jd-I) (d-I)Ylxn) - Yn )ri,d-r =J:J: ..'J: ciand ci are constants.k-IIT (s + j) ds...dsj= 0j:J:i-lWe now quote the following theorem. The constant stepsize, Lagrangian form ofthe corrector polynomial is given byTheorem I If f(x, Y) is a function such that thesolution to (3) exists and the starting errors aresuch that E =max II e· II .,.;;; ch, c > 0, then theoilconstant stepsize h and constant order k formu-lation of 7 (a, b) converges.Pk+l,n+,1 (x) =k~ Li (x, 0) fn + l - i + Lo (x, 0) t'ri+1i=·1(11)i='.li!khiyn(d-r+i) + hr ~ r* fi,d-r n+l-ir-I+ hr r~,d',.-r fri+1> r = 1,2, ... ,d.i = 0(d-r)Yn+1and the correctea values satisfy(Xn+l_i - x n+l _ j )is the LagrangianPERTANI,KA VOL. 8 NO, I, 1985kLi(x, r) = ITj=rrhwhere62kPk,n(X) = ~ Li(x, l)fn +l _ ii=1Proof The Lagrangian form of the polynomialPk,tJx) which interpolates the points, { (xn + I-i'f ) I· = 1· k} associated with the pre-ri+l-i ' ... ,dictor isCONVERGENCE OF THE Dr METHODSSubtracting (15) from (13), the global errore . satisfiedn+.}where+- r times -)-L j (xn + sh, 0)ds ... ds.ken+1 = Aen + ~i= 1a:H~(~n+l-)(y(X .)r ay n+l-1Let HT = (hd, hd- 1, ... , h), a i a d x ddiagonal matrix, = diagonal (y. 0, -y. l' ... ,1, 1,-Yi,d':"'I)' (X7 a d x d diagonal matrix, atdiagonal (-y:" -Y:"I' .•• , -Y:"d -1)'1,0' 1, 1,A a d x d matrix, A =h h 2 hd - 121 (d - 1)10 h h d - 2(d - 2)1l.e.ken+1 =Aen + ~ a7HG~+I_i en+1- ii=lk+ cx~HGTn+l (Aen t ~ a i HG~+l_i en+1 _ i)i=lkYn+1 = AYn + ~ at Hfn+1_ i + a~ Hf\+1i= 1 (13)Then using the above notation, (10) and(11) can be written in the matrix form askPn+1 = AYn + ~ a i Hfn+1_ i (12)i= 1o o8 - * T *where n+l - ( a o HG n+1 Tn + 1 + Tn + 1 )is the local truncation error for -the combinedpredictor-corrector method andarwhere n+l-i. =-aY-~n+1_ i -is a vector between{Y(xn+l_)·and Yn+1- i ,i=I,2, ...kY(xn + 1) and Pn +l' i= O.Therefore,klIen+ll1~IIAllllenll+h*L 2: lIatlillen+l_illi=1Since Y(x n+l) is the true solution at x n+ l'then(14)kAY(xn) + ~ a i Hf(Y(xn+1 _)) + T~+li= 1kY(xn + 1 ) = AY{xn) + ~ a7 Hf(Y(xn+ 1 _))i=lwhere T *n+1 and Tn +1 are vectors for the localtruncation errors of the explicit and implicitmethods respectively.+h*La* IIA [I lien lI+h*2 L2 a*k(~ Ilen+1_ illllllai ll)+8i=lwhere h*= IIH II,L=max IIGn+1_ i II,k8 = max 118 m II, a* = ~ lIai * II,m i=Oka= ~ lIai II.i= 1(16)PERTAN1KA VOL. 8 NO. I, 1985 63M.B. SULEIMAN(18)Here, we introduce {E.} a sequence suchthat II e· II ~E. for J' = 0, I, .'.. , i and for all i =J 1. •0, I, ... , n (n;;"k - 1) and Eo~s given, E o= maxII y (x r) - Yr II, r = 0, 1, ... , k - 1. Notingthat II A II ..,;;;; f + clh and, h* ..,;;;; c2 h, for someconstants cI > 0, c2 > 0 and sufficiently small h,then from (17), we havel1en+1 II..,;;;; U+clh+2hc2w~+c3h2}En+o,(17)for some constant c3 >.0, where c3 > (c I c2w*-tc2 ~L2a*a) ).DefineEn + 1= REn +0where R = 1 + clh + 2hc2w* + c3h2then II en+1 l1"';;;;En+1 and II ej l1"';;;;En <En+1 ,j =0,1, ... , n. (19)4, CONVERGENCE FOR VARIABLE•STEP, VARIABLE ORDER METHODThe proof is very similar to that of constantstepsize and order. In this case the variable step-size is first bounded and then a lemma is stated.Let the variable stepsize hj and order kj, beused on the j - th step, Xj_1 to xj, j = I, 2', ... ,NN, such that b - a = ~ hj.,Let H* = max hj ,j= I . jD*, K are positive constants such that NH* ..,;;;; D*and k j ~ K for all j. Let the ratio of successivestepsize be bounded,h. IJ-11";;; ..,;;;; 71 for all j and kj > 1. Fork j 1,hjthe above restriction is not necessary.The equation in (18), for Rand 0 positivegives rise to the relationshipRn+1 - 1En+1 =Rn+1 Eo + 0 ,n=O, 1, ... (20)R-lFurther, since (1 + X)"'''';;;; exp (nx), x ;;,. 0,using (18), (19) and (20) then (17) becomesRn + 1II en+1 I1";;;Rn+1 Eo'+ -- 0, since (Rn+l ;;,.R-1oRn+1 _ 1) < (Eo +--) exp «n + 1) x),R-1Lemma'Yi,d-r = f: f: ... f: Li (xn + S~+I' 1) ds ... dsis bounded, where 0 < s'< 1.Proof 'Yi,d-r I ..,;;;;fl f S •••JS ~ 18 ji (s)ldso 0 0 j= 1j=FiXn + shn +I - X n+1-jwhere 8ji (s) = --------kand Li (Xn + s~+ I' 1) = n 8ji (S);j= Ij=Fi[Xn+l_i-Xn+l_j I~;t'.=J ~+IandthatWe have only to show that each term.Ieji (s) I is bounded in order for ''Yi,d-r to bebounded. This is done by showingIXn + S~+I -xn+l _ j~. =J ~+Iwhere x = (clh + 2hc2 w*' + c3h2 ). Notingthat (n + 1) h = X p.+1 - a,. the above inequalitybecomesSince 0 = O(h k +'2:) and since Eo < ch,therefore II en+I 11-+ 0 as h -+ O.Therefore the constant stepsize and orderformulation of (7a, b) converges. If .Eo isO(h k +I ), then II en+I II converges with orderk + 1.oIlen + l l1< {Eo+ .}h(cI + 2c2w* + c3 h)exp «xn+1 - a) (c I + 2c2w* + c3 h))64 PERTANIKA VOL. 8 NO. I. 1985CONVERGENCE OF THE DI METHODSare bounded.xn+1 - xn+1_ j = hn+1 + ... +hn+2_ j I{3j < ~+l hn +1this leads to8lIen +1 11< tEo'+ --------}H*(c1 + 2c2L a* + c3 H*< 1 + 71 + 712 + ... + ~-I =(~ - 1)/(71- 1)for allj, since for aRm = 0, 1, ....hn_m ~-m ~+1-m h n---=~+I ~+l-m ~+2-m ~+I< 71m+1 .Hence the convergence of the variable step-size, variable order method (7a, b).REFERENCESIn a similar manner, the denominator ofI () .. (s)1 is bounded below, the bound depend-ing bn the relative sizes of i and j. For j > i,X n+ l - i - xn + l _ j ~+I-i + ~ ....j-+ ... + ~+2-jhn + 1 hn + 1pi _ pi;;;. pi + pi+l + ... + pi-1 = .p-lSimilar bounds exist for j < i.Hence Iai,r I is bounded for all r. There-fore II a i II is bounded. Similarly II at II ISbounded.Theorem 2 If f(x, y) is a function such that thesolution to (3) exists' and the starting errors aresuch that Eo=max l1ei 11< c*H*,i = 0, 1, ... ,DOLD, A. and ECKMAN, B. (1972): Proceedingsof the Conference on the Numerical Solu-tz'on of Ordinary Differential Equations.19 - 20 October 1972, The University ofTexas at Austin. Lecture Notes No. 362.Springer-Varlag.ENRIGHT, W.H., HULL, T.E. and LINDBERG, B.(1974): Comparing Numerical Methods forSystems of Ordinary Differential Equations.Technical Report No. 69, Department ofComputer Science, University of Toronto.GEAR, C.W. (1971): Numerical Initial ValueProblems in Ordinary Differential Equa-tions. Prentice Hall. .GEAR, C.W. (1978): The Stability of NumericalMethods for Second Order Ordinary Diffe-rential Equations. J. SIAM Numer. Anal.15(1): 188-197.+8,k - 1 then the variable stepsize, variable order1 'method (7a, b) converges.then we find thatII en+t II < (1 + c i H* + 2H*c2L* + c3H*2 ) EnProof The remaining argument in proving con-vergence for the variable stepsize, variable ordercase is the same as in Theorem 1. Equation (17)remains true except that H is replaced by H n+Iwhere H~+l =(~+I,h~~1 ""'~+l)'HALL, G and SULEIMAN, M.B. (1981): Stabilityof Adams-Type Formulae for Second OrderOrdinary Differential Equations. IMAJ. ofNumer. Anal. 1: 427 -428.KROGH, F.T. (1969): A Var~able Step VariableOrder Multistep Method for the NumericalSolution of Ordinary Differential Equation.Information Processing 68 - North Hol-land Publishing Company Amsterdam,194-199.RUSSELL, R.D. and SHAMPINE, L.F. (1972): ACollocation Method for Boundary ValueProblems. Num. Math. 19: 1- 28.SHAMPINE, L.F. and GORDON, M.K. (1975):Computer Solution of Ordinary DifferentialEquations. Freeman.ka*;;;' L II ajil,i=Oka;;;' L II ai II andj= 1LetPERTANlKA VOL. 8 NO. I, 1985 65M.B. SULEIMANSERTH, R.W. (1975): Solution of Stiff BoundaryValue Problems by Orthogonal CollocationInt.J.foT Num. Meth. in Eng. 9: 691-199.SULEIMAN, M.B. (1979): Generalised MultistepAdams and Backward DifferentiationMethods for the Solution of Stiff and Non-Stiff ODEs. Ph. D. Thesis, University ofManchester.Received 24 August, 1984)66 PERTANlKA VOL. 8 NO. I, 1985

Convergence of the Variable Order

and Variable Stepsize Direct Integration Methods

for the Solution of the Higher

Order Ordinary Differential Equations

http://psasir.upm.edu.my/2295/1/Convergence_of_the_Variable_Order.pdf

Convergence of the Variable Order and Variable Stepsize Direct Integration Methods for the Solution of the Higher Order Ordinary Differential Equations

Abstract

Similar works

Full text

Available Versions

Universiti Putra Malaysia Institutional Repository