Numerical Analysis Past Paper PDF

# BAS 212 Part 1_ Numerical Analysis Second year (Comm.& power) ## Chapter (5)_Numerical Solution of Ordinary Differential Equations A differential equation is an equation involving an unknown function and one or more of its derivatives. The equation is an ordinary differential equation (ODE) if t...

# BAS 212 Part 1_ Numerical Analysis Second year (Comm.& power) ## Chapter (5)_Numerical Solution of Ordinary Differential Equations A differential equation is an equation involving an unknown function and one or more of its derivatives. The equation is an ordinary differential equation (ODE) if the unknown function depends only on one independent variable. In the examples the independent variable is x, and y(x) is the dependent variable. We will consider methods for solving first order differential equations (if the equation contains derivatives of the nth order, then it is said to be an nth order differential equation). Our typical 1st order differential equation is of the form $y' = f(x, y)$ (1) $y(x0) = yo, a ≤ x ≤ b$ (2) $y(x) = yo$ is known from initial condition. In the following numerical methods we shall divides the interval $[a, b]$ into subintervals by the points $a = X0, X1, X2,..., Xkxk+1,..., xn = b$ where $Xk+1-xk =h$ and $X₁ = x0 + ih$ ## 5.1. Taylor's method: In This method we develop the relation between y and x by expanding y(x) about the point $x = x0$, using the Taylor series: $y(x) = y(x0) + (x – xo) y'(x0) + \frac{(x-x0)^2}{2!}y"(x0)+ \frac{(x-x0)^3}{3!} y'''(x0)$ $+……+ \frac{(x-xo)^n}{n!}y^{(n)}(x0) + Rn$ (3) Letting $(x-x0) = h$, we can write (3) as: $y(xo + h) = y(x0) + hy'(x0) + \frac{h^2}{2!} y"(x0) + \frac{h^3}{3!}y'''(x0)$ $+……+ \frac{h^n}{n!}y^{(n)}(x0)+Rn$ (4) where $Rn$ is the local truncation error and given by $Rn= \frac{(x-xo)^{n+1}}{(n+1)!}y^{(n+1)}(\theta), x_0 < \theta < x$ (5) The solution of (1) over $[a,b]$ is found by adapting formula (3) or (4) to each subinterval $[x,y]$. The derivatives $y^{(i)}(x)$ are obtained by differentiating f(x, y). Note that equation (3) gives y as a function of x but equation (4) gives the value of y numerically. ## Example (1): Solve by Taylor's method the differential equation $y'=1+xy, y(0)=1, 0≤ x ≤ 1$ (6) **Solution:** We may truncate after the third derivative term, then $y(x) = y(a) + (x -a)y'(a)+ \frac{(x-a)^2}{2!}y"(a) + \frac{(x-a)^3}{3!} y'''(a)$ (7) and using $y' = 1+ x y$ we have: $\therefore y" = xy' + y$, $y''' = xy" + y' + y'$ (9) (10) Consider a partition to the interval $[a,b]$ by the points 0,0.25,0.5,1, a=0 **The solution on [0,0.25]** Substitute a = 0 in (9),(10) then we have $y'(0) = 1, y"(0) = 1, y'''(0) = 2$ from equation (7) we have $y(x)=1+x+ \frac{x^2}{2} + \frac{x^3}{3}$, $0<x<0.25$ (11) we find the initial condition for the second interval [0.25, 0.5] by calculate y(0.25) using (11) as following $y(0.25) = 1 +0.25 + \frac{(0.25)^2}{2} + \frac{(0.25)^3}{3} = 1.2865$ Then, $y' = (1 + xy) = 1+ (0.25)(1.2865) = 1.3216$ $y" = xy' + y = (0.25)(1.3216) +1.2865 = 1.6169$ $y''' = xy" + 2y' = (0.25)(.6169) + 2(1.6169) = 3.0474$ And using (7) to get the solution on [0.25, 0.5] $y(x) = 1.2865+ (x -0.25)(1.3216) + \frac{(x-0.25)^2}{2}(1.6169)$ $+ \frac{(x-0.25)^3}{6}(3.0474), 0.25<x<0.5 $ (12) Repeat the previous calculation to find the solution on [0.5, 0.75] $y(0.5) = 1.2865+ (0.5-0.25)(1.3216) + \frac{(0.5-0.25)^2}{2}(1.6169)$ $+ \frac{(0.5-0.25)^3}{6}(3.0474) = 5.0465x10^{-2}$ $y(0.5)= 5.0465x10^{-2}$ $y' = (1 + xy) = 1 + (0.0.5)(5.046510^{-2}) = 1.0252$ $y" = xy' + y = (0.25)(1.0252) + 5.0465x10^{-2} = 0.56307 $ $y''' = xy" + 2y' = (0.25)(0.56307) + 2(1.252) = 2.3319 $ $\therefore y(x) = 5.0465x10^{-2} +(x−0.5)(1.0252)$ $+ (x-0.5)^2 (0.56307) + (x -0.5)^3 (2.3319),$ $0.5<x<0.75$ (13) $y(0.75) = 5.0465x10^{-2} + (0.75-0.5)(1.0252)$ $+ (0.75-0.5)^2 (0.56307)+ (0.75-0.5)^3 (2.3319) = 6.6177×10^{-2}$ $y(0.75) = 6.6177x10^{-2}$ $y' = (1 + xy) = 1 + (0.75)(6.6177x10^{-2}) = 1.0252$ $y" = xy' + y = (0.25)(1.0252) + 5.0465x10^{-2} = 0.56307 $ $y''' = xy" + 2y' = (0.25)(0.56307) + 2(1.252) = 2.3319 $ $y(x) = 6.6177x10^{-2} + (x -0.75)(1.0252)$ $+ (x-0.75)^2 (0.56307) + (x-0.75)^3 (2.3319),$ $0.75<x<1.0$ (14) **Note:** If we can find $y^{(n)} (a)$ we can express the solution as a polynomial of degree n $y(x) = \sum_{k=0}^{n} \frac{(x-a)^k}{k!}y^{(k)}(a)$ ## 5.2. Picard's method: To solve the initial deferential equation: $y' = f(x, y)$ $y(x0) = yo, a ≤x≤b$ We integrate the deferential equation $y' = f(x, y)$ $\int_a^X dy = \int_a^X f(x, y)dx$ $\int_a^X dy = \int_a^X f(x, y)dx$ $y(x) - y(a) = \int_a^X f(x, y)dx$ $y(x) = y(a) + \int_a^X f(x, y)dx$ We cant integrand the function $f(x,y)$ before replacing y as a function of x. We use the initial condition $y(xo) = yo$ to find the first approximation as follows $y₁(x) = y(a) + \int_a^X f(x, y(x_0))dx$ And we use $y₁(x)$ to find the second approximation $y2(x) = y(a) + \int_a^X f(x, y₁(x))dx$ Then $n^{th}$ approximation given by: $y_n(x) = y(a) + \int_a^X f(x,y_{n-1}(x))dx$ (15) By applying this procedure we the approximate solution as a polynomial of degree depend on the number of iteration. To accurate the solution we divide the interval of solution into subinterval $[X0, X1],[X1,X2],..., [Xn-1,xn]$ and apply the method in each subinterval to find the approximate solution as a piecewise polynomial. ## Example(2) By using Picard's method approximate $y'= 1 + xy, y(0) = 1, 0 < x < 1$ at the points $x = 1, x = 0.5$ **Solution** We divides the interval [0, 1] into two subintervals [0, 0.5],[0.5, 1] On the interval [0, 0.5] $y₁(x) = y(0) + \int_a^X f(x,y(0))dx$ $y(0) =1$ Since f(x,y)=1+xy then f(x, y(0)) = 1+ x and the first iteration in the form $y₁(x) =1+\int_{0}^{x} (1+x)dx = 1+x+ \frac{x^2}{2} $ Then $f(x,y₁(x)) = 1 + x(y)= 1 + x(1+x+\frac{x^2}{2}) = 1 + x + x^2 + \frac{x^3}{2} $ And the second iteration in the same interval $y_2(x) = y(0) + \int_{0}^{x} (1 + x + x^2 + \frac{x^3}{2})dx $ $y_2(x)= 1 + x+ \frac{x^2}{2} + \frac{x^3}{3} + \frac{x^4}{8} $ $\therefore y(x) = 1+x+\frac{x^2}{2} + \frac{x^3}{3} + \frac{x^4}{8} $, $0≤x≤0.5$ **On the interval [0.5,1]** $y(x_0)= y(0.5) = 1 +0.5 + \frac{(0.5)^2}{2} + \frac{(0.5)^3}{3} + \frac{(0.5)^4}{8} = 1.6745$ The initial condition for this interval is $y(0.5) ≈ 1.6745$ Using (15) and the data $x = 0.5 y(x) = y(0.5) = 1.6745$ then $y₁(x) = y(a) + \int_a^X f(x, y(x_0))dx$ $= y(0.5)+ \int_{0.5}^{x} 1+ x(y(0.5))dx = 1.6745+ \int_{0.5}^{x} (1+1.6745x)dx$ $= 0.96519+ x +0.83725x^2 $ $0.5<x<1$ $y(x) = \left\{ \begin{array}{ll} 1+x+\frac{x^2}{2} + \frac{x^3}{3} + \frac{x^4}{8}, & 0≤x≤0.5\\ 0.96519 + x + 0.83725x^2, & 0.5<x<1 \end{array} \right. $ Not $\lim_{x \to 0.5^{-}} y(x) = \lim_{x \to 0.5^{+}} y(x)$ which show that $y(x)$ is continuous function $f(x, y₁) = 1 + xy₁ = 1 + x(0.96519+ x +0.83725x^2)$ $\therefore y_2(x) = 1.6745+ \int_{0.5}^{x} (1+ xy_1)dx$ $= 1.6745+ \int_{0.5}^{x} (1+0.96519x + x^2 +0.83725x^3)dx$ $y_2(x) = 0.20931x^4 +0.333x^3 +0.4826x^2 + x +0.9991, 0.5<x<1$ $\therefore y(1) = 3.0243$ ## 5.3. Euler's method: To obtain an approximation to the initial value problem $y' = f(x, y), y(a) = α$ $a≤x≤b$ (16) By choosing a positive integer N and selecting a mesh points $\{x0,x1,...,xy\}$ where $x₁ = a+ih$ for each i = 0,1,2,..., N, the common distance between the points $h = \frac{b-a}{N}$ is called the step size. Suppose that $y(x)$ has two continuous derivative on the interval $[a,b]$ so for each i = 0,1,2,..., N-1 $y(xi+1) = y(xi) + (xi+1 - x₁)y'(x) + \frac{(xi+1-x)^2}{2}y"(ξ_1)$ (17) for some number ξ; where $x₁ <i <xi+1$ If $h=xi+1-xi$ then $y(xi+1) = y(xi) + hy'(x)+ \frac{h^2}{2}y"(ξ_1)$ (18) and since $y(x)$ satisfy the differential equation (16) then $y(xi+1) = y(xi) + hf (xi, y(xi)) + \frac{h^2}{2}y"(ξ_1)$ (19) Euler's method construct that $ y_i = y(x_i)$ for each i = 1,2,3,...,N By deleting the error term in (19) thus $y_0 = a$ $y_{i+1} = y(xi)+h f (xi, y_i), i = 1, 2, ...., n.$ (20) Equation (20) called the difference equation associated with Euler's method. ## Example (3) Approximate the initial value problem $y'=x-y+1,y(0)=1 $(21) in the interval [0, 1] **Solution:** Suppose N = 5 so that h=0.2 and $x₁ = (0.2)i$ Using the fact that $f(x,y)=x-y+1$ gives $y_0 = 1$ $y(xi+1) = y(xi) + hf (xi, Vi)$ $y(xi+1) = yi +0.2(xi - yi + 1)$ $Vi+1 = Vi + 0.2xi – 0.2 y; + 0.2$ $\therefore x₁ = a + ih = 0 + 0.2i$ $Vi+1 = i + (0.2)2i – 0.2y; +0.2$ $= (0.8) y; +0.04i +0.2, i = 0,1,2,3,4.$ Thus the difference equation in the form $\therefore Vi+1 = (0.8)y; +0.04i+0.2, i=0,1,2,3,4.$ $y_1 = 0.2 + (0)(0.04) + (1)(0.8) = ₁ע$ $V_2 = (0.8)y_1 +0.04(1) + 0.2 = 1.04$ $V_3 = (0.8)y_2 + 0.04(2) + 0.2 = 1.11$ $y_4 = (0.8) y_3 + 0.04(3) + 0.2 = 1.20$ $y_5 = (0.8) y_4 +0.04(4) + 0.2 = 1.32$ The exact solution to (21) is $y(x)=x+e^{-x}$ Table (1) shows the comparison between the approximate value at $x_i$ and the actual values | n | x | Yapproximate | Vexact | error | |:---:|:---:|:---:|:---:|:---:| | 0 | 0 | 1 | 1 | 0 | | 1 | 0.2 | 1 | 1.018730753 | 0.018730753 | | 2 | 0.4 | 1.04 | 1.070320046 | 0.030320046 | | 3 | 0.6 | 1.112 | 1.148811636 | 0.036811636 | | 4 | 0.8 | 1.2096 | 1.249328964 | 0.039728964 | | 5 | 1 | 1.32768 | 1.367879441 | 0.040199441 | If N=10 then h = 0.1 and $x₁ = (0.1)i$, table (2) shows the comparison between the approximate value at $x_i$ and the actual values | nx | Yapproximate | Vexact | error | |:---:|:---:|:---:|:---:| | 0 | 0 | 1 | 0 | | 1 | 0.1 | 1 | 1.004837418 | 0.00483741 | | 2 | 0.2 | 1.01 | 1.018730753 | 0.00873075 | | 3 | 0.3 | 1.029 | 1.040818221 | 0.01181822 | | 4 | 0.4 | 1.0561 | 1.070320046 | 0.01422004 | | 5 | 0.5 | 1.09049 | 1.10653066 | 0.01604066 | | 6 | 0.6 | 1.131441 | 1.148811636 | 0.01737063 | | 7 | 0.7 | 1.1782969 | 1.196585304 | 0.01828840 | | 8 | 0.8 | 1.23046721 | 1.249328964 | 0.01886175 | | 9 | 0.9 | 1.28742048 | 1.30656966 | 0.01914917 | | 10 | 1 | 1.34867844 | 1.36787944 | 0.01920100 | ## Example (4): Approximate the initial value problem $y'=1+xy, y(0) = 1,$ $0<x<1$ $h=0.2$ **Solution:** $Vi+1 = Vi + hf (xi, Vi)$ $Vi+1 = yi + hf (xi, vi)$ $= y; + 0.2[1 + xiyi]$ $= y; + 0.2[1 + (0+ ih)y;]$ $= y; + 0.2[1 + (0.2i)yi]$ $Vi+1 = yi +0.2+0.04i yi$, $i = 0,1,2,3,4.$ $y_1 = yo +0.2+0.04(0)yo = 1 +0.2 = 1.2 $ $V_2 = 1 +0.02 +0.04(1)y_1$ $2 = 1.2 +0.2 +0.04(1.2) = 1.448$ $V_3 = 2 +0.2+0.04(2) y_2$ $Уз = 1.448 +0.2 +0.04(2)(1.448) = 1.076384$ $4 = 3 + 0.2 +0.04(3)y_3$ $4 = 1.76384 +0.2 +0.04(3)(1.76384) = 2.1755008$ $V_5 = 5 +0.2+0.04(4)y_4$ $y_5 = 0 = 2.1755008 +0.2 + 0.04(4)(2.1755008) = 2.723580928$ ## 5.4. Runge-Kutta methods: The most common Runge-Kutta method in use is of order four and in difference equation form $y_0 = α$ $Vi+1 = Vi+ \frac{h}{6}[k₁ + 2k2 + 2k3 +k4],$ $k₁ = f(xi, Vi)$ $k₂ = f(xi + \frac{h}{2}, Vi+ \frac{h}{2}k_1)$ $k3 = f(x + \frac{h}{2}, Vi+ \frac{h}{2}k_2)$ $k₁ = f(xi + h, yi + k3)$ $i = 1, 2, ..., N-1$ ## Example (5): Solve $y'=1+y^2$, $0<x<0.3$, $y(0) = 0$ By using Runge-Kutta method considering $h = 0.1$ **Solution** $f(xi, yi)=1+y^2$, $h=0.1$ $k₁ = f(x_0,y_0) = 1 + 0 = 0 = 1$ $k₂ = f(x_0+ \frac{0.1}{2}, y_0 + \frac{0.1}{2}k_1) = f (0.05 +0.05) = 1 + (0.05)^2 = 1.0025$ $k_3 = f(x_0+ \frac{0.1}{2}, y_0 + \frac{0.1}{2}k_2)$ $= f(x_0 + \frac{0.1}{2}y_0 +1.0025) = 1+ (0.050125)^2 = 1.002512516$ $k_4 = f(x_0 +0.1, y_0 +0.1(1.002512516))$ $= 1+ (0.1002512516)^2 = 1.010050313$ $\therefore y_1=y_0+[k₁+2k2+2k3+k₄ ] = 1 [k₁ + 2k2 + 2k3 + k4]$ $\therefore 0.100334589 = ₁ע$ $x_1 = 0.1$ $0.100334589 = ₁ע, h=0.1$ and $f(x,y)=1+ y^2$ at $x = 0.2$ $k₁ = f(x_1,y_1)=1+ y₁^2 = 1+ (0.100334589)^2 = 1.01006703$ $k₂ = f(x_1+ \frac{h}{2}, y_1 + \frac{h}{2}k_1) = 1.022752084$ $k3 = f(x_1+ \frac{h}{2}, y_1 + \frac{h}{2}k_2) = 1.022943825$ $k₁ = f(x_1+h, y_1 + hk_3) = 1.0410585$ $\therefore y_2=1+[k₁+2k2+2k3 + k4]$ $= 0.100334589 + 0.1 [1.01006703 + 2(1.022752084) $ $+2(1.022943825) +1.0410585] = 0.202709878$ $\therefore y_2 = 0.202709878, x_2 = 0.2$ | i | $y_i$ | $y_{i+1}$| exact solution x | error | |---|---|---|---|---| | 0 | 0 | 0.100334589 | 0.099668652 | 0.00065937 | | 1 | 0.100334589 | 0.202709878 | 0.19739556 | 0.005314318 | | 2 | 0.202709878 | 0.309336039 | 0.291456794 | 0.017879245 | ## Example (6): Solve $y' = 1 + xy, y(0)=1,$ $0≤x≤1$ By using Runge-Kutta method. **Solution:** Let $h = 0.2$ | i | $x_i$ | $k_1$ | $k_2$ | $k_3$ | $k_4$ | $y_{i+1}$ | |---|---|---|---|---|---|---| | 0 | 0 | 1 | 1.01 | 1.0101 | 1.040404 | 0.2026868 | | 1 | 0.2 | 1.04053736 | 1.092022161 | 1.093566705 | 1.168560056 | 0.422029305 | | 2 | 0.4 | 1.168811722 | 1.269455239 | 1.274487414 | 1.406156073 | 0.677457742 | | 3 | 0.6 | 1.406474645 | 1.572673644 | 1.584307574 | 1.795455405 | 0.994654158 | | 4 | 0.8 | 1.795723326 | 2.056803841 | 2.080301088 | 2.410714375 | 1.410675743 | If h=0.1 | i | $x_i$ | $k_1$ | $k_2$ | $k_3$ | $k_4$ | $y_{i+1}$ | |---|---|---|---|---|---|---| | 0 | 0 | 1 | 1.0025 | 1.00250625 | 1.010025063 | 0.100333959 | | 1 | 0.1 | 1.010033396 | 1.022625344 | 1.022719784 | 1.040521188 | 0.20268804 | | 2 | 0.2 | 1.040537608 | 1.06367873 | 1.063967994 | 1.092725452 | 0.309163982 | | 3 | 0.3 | 1.092749195 | 1.127330505 | 1.127935677 | 1.16878302 | 0.422031725 | | 4 | 0.4 | 1.16881269 | 1.216212562 | 1.217279059 |1.271879815 | 0.543826321 | | 5 | 0.5 | 1.27191316 | 1.334082088 | 1.335791734 | 1.406443296 | 0.677461389 | | 6 | 0.6 | 1.406476833 | 1.4860604 | 1.488646866 | 1.578428253 | 0.826366716 | | 7 | 0.7 | 1.578456701 | 1.678967163 | 1.682736306 | 1.795712277 | 0.994659648 | | 8 | 0.8 | 1.795727718 | 1.921779129 | 1.927136314 | 2.068635951 | 1.18736289 | | 9 | 0.9 | 2.068626601 | 2.226254509 | 2.233741835 | 2.410737074 | 1.410685496 | ## Example (7): Solve $y' = x-y+1, y(0)=1, 0≤x≤1$ By using Runge-Kutta method. **Solution:** h=0.2 | i | $x_i$ | $k_1$ | $k_2$ | $k_3$ | $k_4$ | $y_{i+1}$ | |---|---|---|---|---|---|---| | 0 | 0 | 0 | 0.1 | 0.09 | 0.182 | 1.018733333 | | 1 | 0.2 | 0.181266667 | 0.26314 | 0.254952667 | 0.33027613