1 / 57

Numerical Solutions of Ordinary Differential Equations

Numerical Solutions of Ordinary Differential Equations. Ordinary Differential Equations. Equations which are composed of an unknown function and its derivatives are called differential equations .

rosewilliams
Download Presentation

Numerical Solutions of Ordinary Differential Equations

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Numerical Solutions of Ordinary Differential Equations Dr.G.SureshKumar@KL University

  2. Ordinary Differential Equations Equations which are composed of an unknown function and its derivatives are called differential equations. Differential equations play a fundamental role in engineering because many physical phenomena are best formulated mathematically in terms of their rate of change. v - dependent variable t - independent variable Dr.G.Suresh Kumar@KL University

  3. Types of Differential Equations • When a function involves one independent variable, the equation is called an ordinary differential equation (ODE). • A partial differential equation(PDE)involves two or more independent variables. Dr.G.Suresh Kumar@KL University

  4. Model Ordinary Differential Equations Dr.G.Suresh Kumar@KL University

  5. Why Numerical solutions of ODE ? Many differential equations cannot be solved analytically, in which case we have to satisfy ourselves with an approximation to the solution. Dr.G.Suresh Kumar@KL University

  6. Initial value problem Initial value problem associated with first order differential equations is Dr.G.Suresh Kumar@KL University

  7. Methods of solving IVP • Taylor series • Picard successive approximation • Euler • Modified Euler • Runge-Kutta (RK) Dr.G.Suresh Kumar@KL University

  8. Taylor Series Method • The Taylor series method is a straight forward adaptation of classic calculus to develop the solution as an infinite series. • The method is not strictly a numerical method but it is used in conjunction with numerical schemes. Dr.G.Suresh Kumar@KL University

  9. Taylor Series Method The Taylor series expansion of y(x) about x = x0 is Dr.G.Suresh Kumar@KL University

  10. Problem#1 The differential equation y'(x) = x + y satisfying y(0)=1. Determine the Taylor series solution. Also compare with the exact solution. {Analytical solution: y(x)=2ex - x -1(Cauchy linear differential equation).} Dr.G.Suresh Kumar@KL University

  11. Solution First, we find derivatives of y' = x + y, Dr.G.Suresh Kumar@KL University

  12. Solution The Taylor series expansion of y(x) about x = x0 is Dr.G.Suresh Kumar@KL University

  13. Comparison of solutions The exact solution is Therefore, Taylor series solution and exact solution are equal. Dr.G.Suresh Kumar@KL University

  14. Picard’s Method The initial value Problem is equivalent to the following integral equation Dr.G.Suresh Kumar@KL University

  15. Picard’s Method The Picard’s nth approximate solution is Dr.G.Suresh Kumar@KL University

  16. Problem#1 Using Picard’s process of successive approximation, obtain a solution up to the fifth approximation of the equation such that y = 1 when x = 0. Check your answer by finding the exact solution. Ans. Exact solution y(x) = 2ex – (x+1). Dr.G.Suresh Kumar@KL University

  17. Application Problem For the free-falling bungee jumper with linear drag, compute the velocity at 10 seconds of a free-falling parachutist of mass 80kg and drag coefficient 10 kg/s at 10 seconds, with initial velocity 20 m/s, using (i) Taylor series method (ii) Picard’s method. Dr.G.Suresh Kumar@KL University

  18. Solution Let v(t) be the velocity and ‘a’ be the acceleration of the parachutist at any time t. Then the total force acting on the parachutist is equal to the sum of the • downward force FD due to gravity and • upward force FU due to air resistance. i.e. F = FD + FU FD = mg FU αv (linear drag) FU = - Cdv(t). Dr.G.Suresh Kumar@KL University

  19. Solution By Newton second law of motion F = ma. ma = mg – Cd v Therefore, Given that m = 80kg, Cd = 10 kg/s, and initial velocity v(0)=20 m/s. Now solve the above IVP using Taylor and Picard’s method. Dr.G.Suresh Kumar@KL University

  20. Euler’s Method Dr.G.Suresh Kumar@KL University

  21. Interpretation of Euler Method y2 y1 y0 x0 x1 x2 x Dr.G.Suresh Kumar@KL University

  22. Interpretation of Euler Method Slope=f(x0,y0) y1 y1=y0+hf(x0,y0) hf(x0,y0) y0 x0 x1 x2 x h Dr.G.Suresh Kumar@KL University

  23. Interpretation of Euler Method y2 y2=y1+hf(x1,y1) Slope=f(x1,y1) hf(x1,y1) Slope=f(x0,y0) y1=y0+hf(x0,y0) y1 hf(x0,y0) y0 x0 x1 x2 x h h Dr.G.Suresh Kumar@KL University

  24. Interpretation of Euler Method Dr.G.Suresh Kumar@KL University

  25. Problem#1 Using Euler’s method, find an approximate value of y corresponding to x=1, given that yʹ = x + y such that y=1 when x=0. Ans: Taking h = 0.1, then y(0.1) = 1.1, y(0.2) = 1.22, y(0.3) = 1.36 y(0.4) = 1.53, y(0.5) = 1.72, y(0.6) = 1.94 y(0.7) = 2.19, y(0.8) = 2.48, y(0.9) = 2.81 y(1.0) = 3.18. Dr.G.Suresh Kumar@KL University

  26. Problem#2 A cup of a coffee originally has temperature of 68oC. The temperature of the ambient is 21oC and the thermal constant is 0.017oC/min. Determine the temperature of the coffee from t = 0 to 10 minutes insteps of 2 minutes. Sol. The differential equation is Tʹ(t) = -k(T - Ta) = -0.017(T - 21) T(0) = 68 Dr.G.Suresh Kumar@KL University

  27. Improvements of Euler’s method • A fundamental source of error in Euler’s method is that the derivative at the beginning of the interval is assumed to apply across the entire interval. Two simple modifications are available to circumvent this shortcoming: • Heun’s (Modified Euler) Method • Midpoint (Improved Polygon) Method Dr.G.Suresh Kumar@KL University

  28. Modified Euler’s Method • To improve the estimate of the slope, determine two derivatives for the interval: • (1) At the initial point • (2) At the end point • The two derivatives are then averaged to obtain an improved estimate of the slope for the entire interval. Dr.G.Suresh Kumar@KL University

  29. Euler Modified method • To improve the estimate of the slope, determine two derivatives for the interval: • At the initial point • At the end point • The two derivatives are then averaged to obtain an improved estimate of the slope for the entire interval.

  30. Heun’s method (improved) Original Huen’s: Note that the corrector can be iterated to improve the accuracy of yi+1 .

  31. Problem#1 Using Modified Euler’s method, find an approximate value of y when x= 0.3, given that yʹ = x + y and y=1 when x=0. Sol. Compare the IVP with the general IVP yʹ = f(x, y), y(x0) = x0, we have f(x, y) = x + y, x0 = 0 and y0 = 1. Taking the step size h = 0.1, then x0 = 0, x1 = 0.1, x2 = 0.2, x3 = 0.3. Now we find the corresponding value of y, using modified Euler’s method.

  32. Solution Step-1. To find y1 = y(x1) = y(0.1) Predictor formula (Euler’ formula) y1(0) = y0 + h f(x0, y0) = 1 + 0.1 f(0, 1) = 1+0.1(0+1) =1.1 Corrector formula y1(k) = y0 + (h/2) [f(x0, y0) + f(x1, y1(k-1))] y1(1) = y0 + (h/2) [f(x0, y0) + f(x1, y1(0))] = 1 + (0.1/2) [f(0,1) + f(0.1, 1.1)] = 1 + 0.05 [0+1 + 0.1+1.1] =1.11

  33. Solution y1(2) = y0 + (h/2) [f(x0, y0) + f(x1, y1(1))] = 1 + (0.1/2) [f(0,1) + f(0.1, 1.11)] = 1 + 0.05 [0+1 + 0.1+1.11] =1.1105. y1(3) = y0 + (h/2) [f(x0, y0) + f(x1, y1(2))] = 1 + (0.1/2) [f(0,1) + f(0.1, 1.1105)] = 1 + 0.05 [0+1 + 0.1+1.1105] =1.1105. Therefore, y1= 1.1105.

  34. Solution Step-2. To find y2 = y(x2) = y(0.2) Predictor formula (Euler’ formula) y2(0) = y1 + h f(x1, y1) = 1.1105 + 0.1 f(0.1, 1.1105) = 1.1105+0.1(0.1 + 1.1105) =1.2316. Corrector formula y2(k) = y1 + (h/2) [f(x1, y1) + f(x2, y2(k-1))] y2(1) = y1 + (h/2) [f(x1, y1) + f(x2, y2(0))] = 1.2316 + (0.1/2) [f(0.1,1.1105) + f(0.2, 1.2316)] = 1.2316 + 0.05 [0.1+1.1105 + 0.2+1.2316] =1.2426.

  35. Solution y2(2) = y1 + (h/2) [f(x1, y1) + f(x2, y2(1))] = 1.1105 + (0.1/2) [f(0.1,1.1105) + f(0.2, 1.2426)] = 1.1105 + 0.05 [0.1+1.1105 + 0.2+1.2426] =1.2432. y2(3) = y1 + (h/2) [f(x1, y1) + f(x2, y2(2))] = 1.1105 + (0.1/2) [f(0.1,1.1105) + f(0.2, 1.2432)] = 1.1105 + 0.05 [0.1+1.1105 + 0.2+1.2432] =1.2432. Therefore, y2 = y(0.2)= 1.2432. Similarly, y3 = y(0.3) = 1.4004.

  36. Application Problem • A storage tank (Fig.) contains a liquid at depth y where y = 0 when the tank is half full. Liquid is withdrawn at a constant flow rate Q=500m3/hto meet demands. The contents are resupplied at a sinusoidal rate 3Q sin2(t). The surface area of the tank is 1200 m2. Determine the depth of the tank from t= 0 to t=6 in steps of 2 hours, using modified Euler’s method. Dr.G.Suresh Kumar@KL University

  37. Problem#1 Q =500 m3/h 3Q sin2(t) Dr.G.Suresh Kumar@KL University

  38. Solution Let y(t) be the depth of the tank at any time t and A be the surface area of the tank. The depth of the tank is changed when volume of the liquid in the tank changed. Therefore, Dr.G.Suresh Kumar@KL University

  39. Solution Since A is constant, it implies that Given that Q=500m3/h, A =1200 m2 and y(0) = 0. Substituting these value, we have Now solving the above initial value problem, using modified Euler’s method. Dr.G.Suresh Kumar@KL University

  40. Runge - Kutta Method Dr.G.Suresh Kumar@KL University

  41. Carl Runge(1856-1927) Martin Wilhelm Kutta(1867-1944) Dr.G.Suresh Kumar@KL University

  42. Introduction • Runge-Kutta methods are very popular because of their good efficiency; and are used in most computer programs for differential equations. • These methods agree with Taylor’s series solution up to the term in hr, where r differs method to method and is called the order of the method. • Euler’s and modified Euler’s are called the first and second order Runge-Kutta methods. • The fourth order R-K method is most commonly used and is often referred to as R-K method only. Dr.G.Suresh Kumar@KL University

  43. Working Rule To find increment ‘k’ of y corresponding to an increment h of x by R-K method from the initial value problem y´ = f(x, y) with y(x0) = y0 is as follows Calculate k1 = h f(x0, y0) k2 = h f(x0+0.5h, y0+0.5k1) k3 = h f(x0+0.5h, y0+0.5k2) k4 = h f(x0+h, y0+k3) Dr.G.Suresh Kumar@KL University

  44. Working Rule • Finally compute k = 1/6 (k1 + 2k2 + 2k3 + k4) is the weighted mean of k1, k2, k3 and k4. • Therefore the required approximate value is y1 = y(x1) = y(x0 + h) = y0 + k. Dr.G.Suresh Kumar@KL University

  45. Graphical Representation k2 k4 k3 k1 xi xi + h/2 xi + h Dr.G.Suresh Kumar@KL University

  46. Problem#1 • Using Runge-Kutta method of fourth order, solve yʹ = (y2 – x2)/(y2 + x2) with y(0)=1 at x = 0.2, 0.4. Sol : - Compare the given problem with the general first order initial value problem yʹ=f(x, y) satisfying y(x0)=y0, we have f(x, y) = (y2 – x2)/(y2 + x2) x0 = 0, y0 =1. Here we have to find value of y at x=0.2, 0.4. Dr.G.Suresh Kumar@KL University

  47. Solution Taking step size h = 0.2. Step(1) To find y(0.2) k1 = h f(x0, y0) = 0.2 f(0, 1) = 0.2 (12-02)/(12+02) = 0.2 k2 = h f(x0+0.5h, y0+0.5k1) = 0.2 f(0+0.5(0.2), 1+ 0.5(0.2)) = 0.2 f(0.1, 1.1) = 0.19672 Dr.G.Suresh Kumar@KL University

  48. Solution k3 = h f(x0+0.5h, y0+0.5k2) = 0.2 f(0+0.5(0.2), 1+0.5(0.19672)) = 0.2 f(0.1, 1.09836) = 0.1967 k4 = h f(x0+h, y0+k3) = 0.2 f(0+0.2, 1+0.1967) = 0.2 f(0.2, 1.1967) = 0.1891 Dr.G.Suresh Kumar@KL University

  49. Solution k = 1/6 (k1 + 2k2 + 2k3 + k4) = 1/6 (0.2 +2(0.19672)+2(0.1967)+0.1891) = 0.19599. Hence y1 = y0 + k y(0.2) = 1 + 0.19599 = 1.19599. Step(2) To find y(0.4) k1 = h f(x1, y1) = 0.2 f(0.2, 1.19599) = 0.1891 Dr.G.Suresh Kumar@KL University

  50. solution k2 = h f(x1+0.5h, y1+0.5k1) = 0.2 f(0.2+0.5(0.2), 1.19599+ 0.5(0.1891)) = 0.2 f(0.3, 1.29054) = 0.1795 k3 = h f(x1+0.5h, y1+0.5k2) = 0.2 f(0.2+0.5(0.2), 1.19599+0.5(0.1795)) = 0.2 f(0.3, 1.28574) = 0.1793 Dr.G.Suresh Kumar@KL University

More Related