Question 10.7: Solving a system of two first-order OD Es using Euler's expl......
Solving a system of two first-order OD Es using Euler’s explicit method and second-order Runge-Kutta method.
Consider the following initial value problem consisting of two first-order ODEs:
\begin{gathered}\frac{d y}{d x}=(-y+z) e^{(1-x)}+0.5 y \text { with the initial condition } y(0)=3 & (10.127)\\\frac{d z}{d x}=y-z^{2} \text { with the initial condition } z(0)=0.2 & (10.128)\end{gathered}
on the domain from x=0 to x=3.
(a) Solve the system for the first three steps by hand with Euler’s explicit method using h=0.25.
(b) Solve the system for the first two steps by hand with the second-order Runge-Kutta method (modified Euler version) using h=0.25.
(c) Write a MATLAB program in a script file that solves the system with the second-order RungeKutta method (modified Euler version) using h=0.1.
Show the results from the three parts in a plot.
(a) Solution by hand with Euler’s method: The first point of the solution is x_{1}=0, y_{1}=3, and z_{1}=0.2, which is the point where the initial conditions are given. For the first point i=1.
The rest of the solution is determined by using Eqs. (10.112)-(10.114). In the present problem these equations have the form:
\begin{gathered} x_{i+1}=x_{i}+h=x_{i}+0.25 & (10.129)\\ y_{i+1}=y_{i}+f_{1}\left(x_{i}, y_{i}, z_{i}\right) h=y_{i}+\left[\left(-y_{i}+z_{i}\right) e^{\left(1-x_{i}\right)}+0.5 y_{i}\right] h & (10.130)\\ z_{i+1}=z_{i}+f_{2}\left(x_{i}, y_{i}, z_{i}\right) h=z_{i}+\left[y_{i}-z_{i}^{2}\right] h & (10.131)\end{gathered}
Equations (10.129)-(10.131) are next applied three times with i=1,2,3.
First step: For the first subinterval i=1. Equations (10.129)-(10.131) give:
x_{2}=x_{1}+0.25=0+0.25=0.25y_{2}=y_{1}+\left[\left(-y_{1}+z_{1}\right) e^{\left(1-x_{1}\right)}+0.5 y_{1}\right] h=3+\left[(-3+0.2) e^{(1-0)}+0.5 \cdot 3\right] 0.25=1.472
z_{2}=z_{1}+\left[y_{1}-z_{1}^{2}\right] h=0.2+\left[3-0.2^{2}\right] 0.25=0.94
The second point of the solution is: x_{2}=0.25, y_{2}=1.472, and z_{2}=0.94.
Second step: For the second step i=2. Equations (10.129)-(10.131) give:
x_{3}=x_{2}+0.25=0.25+0.25=0.5y_{3}=y_{2}+\left[\left(-y_{2}+z_{2}\right) e^{\left(1-x_{2}\right)}+0.5 y_{2}\right] h=1.472+\left[(-1.472+0.94) e^{(1-0.25)}+0.5 \cdot 1.472\right] 0.25=1.374
z_{3}=z_{2}+\left[y_{2}-z_{2}^{2}\right] h=0.94+\left[1.472-0.94^{2}\right] 0.25=1.0871
The third point of the solution is: x_{3}=0.5, y_{3}=1.374, and z_{3}=1.087.
Third step: For the third step i=3. Equations (10.129)-(10.131) give:
x_{4}=x_{3}+0.25=0.5+0.25=0.75y_{4}=y_{3}+\left[\left(-y_{3}+z_{3}\right) e^{\left(1-x_{3}\right)}+0.5 y_{3}\right] h=1.374+\left[(-1.374+1.087) e^{(1-0.5)}+0.5 \cdot 1.374\right] 0.25=1.427
z_{4}=z_{3}+\left[y_{3}-z_{3}^{2}\right] h=1.087+\left[1.374-1.087^{2}\right] 0.25=1.135 The fourth point of the solution is: x_{3}=0.75, y_{3}=1.427, and z_{3}=1.135.
The solution obtained with Euler’s explicit method is shown in the figure at the end of the solution of part (c)
(b) Solution by hand with second-order Runge-Kutta method: The first point of the solution is x_{1}=0, y_{1}=3, and z_{1}=0.2, which is the point where the initial conditions are given. For the first point i=1.
The rest of the solution is determined by using Eqs. (10.124)-(10.126). In the present problem these equations have the form:
\begin{gathered} x_{i+1}=x_{i}+h=x_{i}+0.25 &(10.132)\\ K_{y, 1}=\left(-y_{i}+z_{i}\right) e^{\left(1-x_{i}\right)}+0.5 y_{i} &(10.133)\\ K_{z, 1}=y_{i}-z_{i}^{2} \\ K_{y, 2}=\left(-y_{e s t}+z_{e s t}\right) e^{\left(1-x_{i+1}\right)}+0.5 y_{e s t} & (10.134)\\ K_{z, 2}=y_{e s t}-z_{e s t}^{2} \end{gathered}
K_{n,1}=f_{n}(x_{i}, y_{1,i},…,y_{n,i} ) (10.124)
y_{n,i+1}=y_{n,i}+\frac{1}{2}(K_{n,1}+K_{n,2})h (10.126)
Note: To simplify the equations, the quantities y_{e s t} and z_{e s t} are introduced. They are calculated after K_{\mathrm{y}, 1} and K_{z, 1} are determined, and are used for evaluating K_{\mathrm{y}, 2} and K_{z, 2}.
where: y_{e s t}=y_{i}+K_{y, 1} h, and z_{e s t}=z_{i}+K_{z, 1} h.
\begin{gathered}y_{i+1}=y_{i}+\frac{1}{2}\left(K_{y, 1}+K_{y, 2}\right) h \\z_{i+1}=z_{i}+\frac{1}{2}\left(K_{z, 1}+K_{z, 2}\right) h\end{gathered} (10.135)
Equations (10.132)-(10.135) are next applied three times with i=1,2,3.
First step: For the first interval i=1. Equations (10.132)-(10.135) give:
x_{2}=x_{1}+0.25=0+0.25=0.25K_{y, 1}=\left(-y_{1}+z_{1}\right) e^{\left(1-x_{1}\right)}+0.5 y_{1}=(-3+0.2) e^{(1-0)}+0.5 \cdot 3=-6.111
K_{z, 1}=y_{1}-z_{1}^{2}=3-0.2^{2}=2.96
y_{e s t}=y_{1}+K_{y, 1} h=3+(-6.111) 0.25=1.472
z_{e s t}=z_{1}+K_{z, 1} h=0.2+2.96 \cdot 0.25=0.94
K_{y, 2}=\left(-y_{e s t}+z_{e s t}\right) e^{\left(1-x_{i+1}\right)}+0.5 y_{e s t}=(-1.472+0.94) e^{(1-0.25)}+0.5 \cdot 1.472=-0.3902
K_{z, 2}=y_{e s t}-z_{e s t}^{2}=1.472-0.94^{2}=0.5884
y_{2}=y_{1}+\frac{1}{2}\left(K_{y, 1}+K_{y, 2}\right) h=3+\frac{1}{2}[-6.111+(-0.3902)] 0.25=2.187
z_{2}=z_{1}+\frac{1}{2}\left(K_{z, 1}+K_{z, 2}\right) h=0.2+\frac{1}{2}[2.96+0.5884] 0.25=0.6436
The second point of the solution is: x_{2}=0.25, y_{2}=2.187, and z_{2}=0.6436.
Second step: For the second interval i=2. Equations (10.132)-(10.135) give:
x_{3}=x_{2}+0.25=0.25+0.25=0.5K_{y, 1}=\left(-y_{2}+z_{2}\right) e^{\left(1-x_{2}\right)}+0.5 y_{2}=(-2.187+0.6436) e^{(1-0.25)}+0.5 \cdot 2.187=-2.173
K_{z, 1}=y_{2}-z_{2}^{2}=2.187-0.6436^{2}=1.773
y_{e s t}=y_{2}+K_{y, 1} h=2.187+(-2.173) 0.25=1.644
z_{e s t}=z_{2}+K_{z, 1} h=0.6436+1.773 \cdot 0.25=1.087
K_{y, 2}=\left(-y_{e s t}+z_{e s t}\right) e^{\left(1-x_{i+1}\right)}+0.5 y_{e s t}=(-1.644+1.087) e^{(1-0.5)}+0.5 \cdot 1.644=-0.09634
K_{z, 2}=y_{e s t}-z_{e s t}^{2}=1.644-1.087^{2}=0.4624
y_{3}=y_{2}+\frac{1}{2}\left(K_{y, 1}+K_{y, 2}\right) h=2.187+\frac{1}{2}[-2.173+(-0.09634)] 0.25=1.903
z_{3}=z_{2}+\frac{1}{2}\left(K_{z, 1}+K_{z, 2}\right) h=0.6436+\frac{1}{2}[1.773+0.4624] 0.25=0.9230
The third point of the solution is: x_{3}=0.5, y_{3}=1.903, and z_{3}=0.9230.
The solution points that were obtained with the second-order Runge-Kutta method are shown in the figure at the end of the solution of part (c).
Script File
c) Solution with a computer program that uses the second-order Runge-Kutta method. First, a user-defined function, named Sys2ODEsRK2, that solves a system of two ODEs using the secondorder Runge-Kutta method is written. This function uses two other user-defined functions (listed in the function definition line as ODE1 and ODE2) that calculate the value of d y / d x and d z / d x in each step. The functions are named odeExample 7 d y d x and odeExample 7 \mathrm{dz} d x. The solution itself is done in a MATLAB program in a script file. The program also produces a plot that shows the solution from part (c) and the solutions calculated in parts (a) and (b).
Program 10-5: User-defined function. Solving a system of two first-order ODEs using a second-order Runge-Kutta method (modified Euler version).
function [x, y, z] = Sys2ODEsRK2(ODE1,ODE2,a,b,h,yINI,zINI)
% Sys2ODEsRK2 solves a system of two first-order initial value ODEs using
% second-order Runge-Kutta method.
% The independent variable is x, and the dependent variables are y and z.
% Input variables:
% ODE1 Name for the function that calculates dy/dx.
% ODE2 Name for the function that calculates dz/dx.
% a The first value of x.
% b The last value of x.
% h The size of a increment.
% yINI The initial value of y.
% zINI The initial value of z.
% Output variables:
% x A vector with the x coordinate of the solution points.
% y A vector with the y coordinate of the solution points.
% z A vector with the z coordinate of the solution points.
x(1) = a; y(1) = yINI; z(1) = zINI;
N = (b - a)/h;
for i = 1:N
x(i + 1) = x(i) + h;
Ky1 = ODE1 (x (i) , y (i) , z (i)) ;
Kz1 = ODE2 (x (i) , y (i) , z (i)) ;
Ky2=0DE1(x(i+l) ,y(i) +Ky1*h,z(i) +Kz1*h);
Kz2 = ODE2 (x (i + 1) , y (i) + Ky1*h, z (i) + Kz1*h) ;
y (i + 1) =y (i) + (Ky1 + Ky2) *h/2;
z(i+l) =z(i) + (Kz1+Kz2)*h/2;
end
Listed next are two user-defined functions that calculate the value of dyldx, Eq. (10.127), and dzldx, Eq. (10.128). The functions are named odeExample7dydx and odeExample7dzdx.
function dydx = odeExample7dydx(x,y,z)
dydx= (-y+ z) *exp (1 - x) + 0. 5*y;
function dzdx = odeExample7dzdx(x,y,z)
dzdx=y - z^2;
The MATLAB program in a script file that uses the above functions for solving the problem is listed next. The program also displays a plot that shows the solution from part (c) and the solutions that were obtained in parts (a) and (b).
% Solving Chapter 10 Example 7
clear
a = 0; b = 3; yINI = 3; zINI = 0.2; h = 0.1;
% Data from part (a)
xa = [0 0.25 0.5 0.75];
ya = [3 1.472 1.374 1.427);
za = [0.2 0.94 1.087 1.135);
% Data from part (b)
xb = [0 0.25 0.5];
yb = [3 2.187 1.903);
zb = [0.2 0.6436 0.9230);
plot (x, y, '-k' , x, z, '-r' , xa, ya, '*k' , xa, za, '*r' , xb, yb, 'ok' , xb, zb, 'or')
The figure on the right shows the numerical solution obtained for the whole domain in part (c) (solid lines), the first four points from part (a) (star markers), and the three points calculated in part (b) (circle markers).
