Numerical Solution of Ordinary Differential Equations PDF

UNIT V Numerical solution of Ordinary Differential equations Solution by Taylor’s series-Picard’s Method of successive Approximations-Euler’s Method-Runge-Kutta Methods. Numerical solutions of Laplace equation using finite difference approximation. 1. Ta...

UNIT V Numerical solution of Ordinary Differential equations Solution by Taylor’s series-Picard’s Method of successive Approximations-Euler’s Method-Runge-Kutta Methods. Numerical solutions of Laplace equation using finite difference approximation. 1. Taylor's Series method Consider the one dimensional initial value problem y' = f(x, y), y(x0) = y0 where f is a function of two variables x and y and (x0 , y0) is a known point on the solution curve. If the existence of all higher order partial derivatives is assumed for y at x = x0, then by Taylor series the value of y at any neighboring point x+h can be written as y(x0+h) = y(x0) + h y'(x0) + h2 /2 y''(x0) + h3/3! y'''(x0) +...... where ' represents the derivative with respect to x. Since at x0, y0 is known, y' at x0 can be found by computing f(x0,y0). Similarly higher derivatives of y at x0 also can be computed by making use of the relation y' = f(x,y)  y'' = fx + fyy' y''' = fxx + 2fxyy' + fyy y'2 + fyy'' and so on. Then y(x0+h) = y(x0) + h f + h2 ( fx + fyy' ) / 2! + h3 ( fxx + 2fxyy' + fyy y'2 + fyy'' ) / 3! + o(h4) Hence the value of y at any neighboring point x0+ h can be obtained by summing the above infinite series. However, in any practical computation, the summation has to be terminated after some finite number of terms. If the series has been terminated after the Pth derivative term then the approximated formula is called the Taylor series approximation to y of order P and the error is of order P+1. The same can be repeated to obtain y at other points of x in the interval [x0, xn] in a marching process. Algorithm 1. Specify x0, xn, y0, h ( (x0, y0) Initial point, xn point where the solution is required ,h the step length to be used in the marching process ) 2.Repeat compute f(xi, yi), f'(xi, yi), f''(xi, yi)... compute y(xi+h) = y(xi) + h f(xi, yi) + h2 /2 f'(xi, yi) + h3/3! f''(xi, yi) +... xi = xi + h until xi = xn However, in practice one usually computes two sets of approximations using step sizes h and h/2 and compares the solutions For p = 4, E4 = c * h4 and the same with step size h/2, E4 = c * (h/2)4, that is if the step size is halved the error is reduced by an order of 1/16. Solve the initial value problem y' = -2xy2, y(0) = 1 for y at x = 1 with step Example 1 length 0.2 using Taylor series method of order four. Solution: Given y' = f(x,y) = -2xy2 y'' = -2y2 - 4xyy' y''' = -8yy' - 4xy'2 - 4xyy'' yiv = -12y'2 - 12yy'' - 12xy'y'' - 4xyy''' yv = -48y'y'' - 16yy''' -12xy''2 - 16xy'y''' - 4xyyiv The forth order Taylor's formula is y(xi+h) = y(xi) + h y'(xi, yi) + h2 y''(xi, yi)/2! + h3 y'''(xi, yi)/3! + h4 yiv(xi, yi)/4! +... given at x=0, y=1 and with h =.2 we have y' = -2(0)(1)2 = 0.0 y'' = -2(1)2 - 4(0)(1)(0) = -2 y''' = -8(1)(0) - 4(0)(0)2 - 4(0)(1)(-2) = 0.0 yiv = -12(0)2 - 12(1)(-2) - 12(0)(0)(-2) - 4(0)(1)(0) = 24 y(0.2) = 1 +.2 (0) +.22 (-2)/2! + 0 +.24 (24)/4! =.9615 now at x =.2 we have y =.9615 y' = -0.3699, y'' = -1.5648, y''' = 3.9397 and yiv = 11.9953 then y(0.4) = 1 +.2 (-.3699) +.22 (-1.5648)/2! +.23(3.9397)/3! +.24 (11.9953)/4! = 0.8624 y(0.6) = 1 +.2 (-.5950) +.22 (-0.6665)/2! +.23 (4.4579)/3! +.24 (-5.4051)/4! = 0.7356 y(0.8) = 1 +.2 (-.6494) +.22 (-0.0642)/2! +.23 (2.6963)/3! +.24 (-10.0879)/4! = 0.6100 y(1.0) = 1 +.2 (-.5953) +.22 (-0.4178)/2! +.23 (0.9553)/3! +.24 (-6.7878)/4! = 0.5001 now at x = 1.0 we have y =.5001 y' = -0.5001, y'' = 0.5004, y''' = -.000525 and yiv = -3.0005 Error in the approximation E4 = h4 (y4(1.0) - y4(0.8))/ 5! = h4 (-3.0005-11.9953)/ 5! = -1.9994e-004 Analytical solution y(x) = 1/(1+x2) at x = 1, y =.5 Alternate method: By Taylor series method y(x0 + x - x0) = y0 + (x-x0)y'(x0) + (x-x0)2y''(x0)/2! + (x-x0)3y'''(x0)/3! + (x- x0)4yiv(x0)/4! +... Since x0 = y' = y''' = 0, y = 1, y'' = -2 and yiv = 24 at x = 0, y( x ) = 1. - x2 + x4 +...  Using Taylor series method of order four solve the initial value problem y' Example 2 = (x - y)/2, on [0, 3] with y(0) = 1. Compare solutions for h = 1, 1/2, 1/4 and S 1/8. Solution: Given y' = f(x,y) = (x - y)/2 y'' = (1 - y')/2, y''' = -y''/2, yiv = -y'''/2 The forth order Taylor's formula is y(xi+h) = y(xi) + h y'(xi, yi) + h2 y''(xi, yi)/2! + h3 y'''(xi, yi)/3! + h4 yiv(xi, yi)/4! +... = 1 + h (xi - yi)/2 - h2 (1 - y'i)/(2*2!) - h3 y''i/(2*3!) - h4 y'''i/(2*4!) + o(h5) =.5* ( 2 + h (xi - yi) -.5h2 (1 - y'i) -.1667h3 y''i -.0417h4 y'''i ) + o(h5) Comparison of the solution for various h yi Exact xi h=1 h =.5 h =.25 h =.125 yi 0.000 1.000000 1.000000 1.000000 1.000000 1.000000 0.125 0.943239 0.943240 0.250 0.897492 0.897491 0.897492 0.375 0.862087 0.862087 0.500 0.836426 0.836404 0.836402 0.836402 0.750 0.811870 0.811868 0.811868 1.000 0.820315 0.819629 0.819594 0.819592 0.819592 1.500 0.917142 0.917102 0.917100 0.917100 2.000 1.104513 1.103683 1.103641 1.103639 1.103639 2.500 1.359558 1.359517 1.359515 1.359514 3.000 1.670186 1.669431 1.669393 1.669391 1.669391 2. Euler's Method : Though in principle it is possible to use Taylor's method of any order for the given initial value problem to get good approximations, it has few draw backs like   The scheme assumes the existence of all higher order derivatives for the given function f(x,y) which is not a requirement for the existence of the solution for any first order initial value problem.  Even the existence of these higher derivatives is guaranteed it may not be easy to compute them for any given f(x,y).  Because of the usage of higher order derivatives in the formula it is not convenient to write computer programs , that is the method is more suited for hand calculations. To overcome these difficulties, Euler developed a scheme by approximating y' in the givenivp. The scheme is as follows: The derivative term in the first order ivp y' = f(x, y) , y(x0) = y0 is approximated by making use of Taylor series approximation of the dependent variable y(x) at the point xi+1. Find y(0.5) if y is the solution of IVP y' = -2x-y, y(0) = -1 using Euler's Example 1 method with step length 0.1. Also find the error in the approximation. Solution: f(x, y) = -2x - y, y1 = y0 + h f(x0, y0) = -1 + 0.1* (-2*0 - (-1)) = -0.8999 y2 = y1 + h f(x1, y1) = -0.8999 + 0.1* (-2*0 - (-0.8999)) = -0.8299 y3 = y2 + h f(x2, y2) = -0.8299 + 0.1* (-2*0 - (-0.8299)) = -0.7869 y4 = y3 + h f(x3, y3) = -0.7869 + 0.1* (-2*0 - (-0.7869)) = -0.7683 y5 = y4 + h f(x4, y4) = -0.7683 + 0.1* (-2*0 - (-0.7683)) = -0.7715 Truncation error in the approximation = ( h2/2 ) f''( ), where 0 <