Question 10.6: A user-defined function for solving a first-order ODE using ......
A user-defined function for solving a first-order ODE using the fourthorder Runge-Kutta method.
Write a user-defined MATLAB function that solves a first-order ODE using the classical fourthorder Runge-Kutta method. Name the function [x, y]= odeRK 4(O D E, a, b, h), where ODE is the name of a user-defined function that calculates the derivative d y / d x.
Use the function odeRK4 to solve: using h=0.5.
\frac{d y}{d x}=-1.2 y+7 e^{-0.3 x} \text { from } x=0 \text { to } x=2.5 \text { with the initial condition } y=3 \text { at } x=0 .
Compare the results with the exact (analytical) solution: y=\frac{70}{9} e^{-0.3 x}-\frac{43}{9} e^{-1.2 x}.
Script File
To solve the problem, a user-defined MATLAB function called odeRK4, which solves a first-order initial value ODE, is written. The function is then used in a script file, which also generates a plot that shows a comparison between the numerical and the exact solutions. The ODE itself is written in a separate user-defined function that is used by the odeRK4 function.
Program 10-4: User-defined function. Solving first-order ODE using Runge-Kutta fourth-order metho4.
function [x, y] = odeRK4(0DE,a,b,h,yini)
% odeRK4 solves a first order initial value ODE usinq Runqe-Kutta fourth order method.
% Input variables:
% ODE Name for the function that calculates dy/dx.
% a The first value of x.
% b The last value of x.
% h Step size.
% yIni The value of the solution y at the first point (initial value).
% Output variables:
% x A vector with the x coordinate of the solution points.
% y A vector with the y coordinate of the solution points.
x(1) = a; y(1) = yIni
n = (b - a)/h;
for i = 1:n
x(i + 1) = x(i) + h;
K1 = ODE(x(i) ,y(i));
xhalf = x(i) + h/2;
yK1 = y(i) + K1*h/2;
K2 = ODE(xhalf,yK1);
yK2 = y(i) + K2*h/2;
K3 = ODE(xhalf,yK2);
yK3 = y(i) + K3*h;
K4 = ODE(x(i + 1),yK3);
y (i + 1) =y (i) + (K1+ 2*K2 + 2*K3 + K4) *h/6;
end
Calculate K_1, Eq. (10.87).
\begin{aligned} & K_1=f\left(x_i, y_i\right) \\ & K_2=f\left(x_i+\frac{1}{2} h, y_i+\frac{1}{2} K_1 h\right) \\ & K_3=f\left(x_i+\frac{1}{2} h, y_i+\frac{1}{2} K_2 h\right) \\ & K_4=f\left(x_i+h, y_i+K_3 h\right) \end{aligned} (10.87)
Calculate K_2, Eq. (10.87).
Calculate K_3, Eq. (10.87).
Calculate K_4, Eq. (10.87).
Calculate the next value of the dependent I variable, Eq. (10.86).
y_{i+1}=y_i+\frac{1}{6}\left(K_1+2 K_2+2 K_3+K_4\right) h (10.86)
The following script file uses the function odeRK4 for solving the ODE.
a = 0; b = 2.5;
h = 0.5; yIni = 3;
[x,y]=odeRK4(@Chap10Exmp6ODE,a,b,h,yIni)
xp=a:0.1:b;
yp=70/9*exp(-0.3*xp) - 43/9*exp(-1.2*xp);
plot(x,y, '*r' ,xp,yp)
yExact=70/9*exp(-0.3*x) - 43/9*exp(-1.2*x)
error= yExact - y
The user-defined function Cha pl 0Exmp60DE used in the argument of the function odeRK4 calculates the value of dyldx:
function dydx = Chap10Exmp6ODE(x,y)
dydx = -1.2*y + 7*exp(-0.3*x);
When the script file is executed, the following data is displayed in the Command Window. In addition, a plot that shows the points of the numerical solution and the exact solution is displayed in the Figure Window.
The results (the error vector and the figure) show that the numerical solution has a very small error, even though the step size is quite large.
