Topic: Numerical Solution to Ordinary Differential Equation
Subject: Numerical Methods
Tool: Scilab
The 4th -order Runge-Kutta methods (RK4) is a popular and most widely used numerical technique of approximating solution to the initial value problem (IVP) of 1st -order ordinary differential equation (ODE). This article presents the equations and algorithm to apply the RK4 to a problem using Scilab scripts.
The RK4 method estimates the solution \(y=f(x)\) of an ODE with initial value condition (IVC) of the form:
$$\frac{dy}{dx}\,=\,f(x,y), \quad y(x_0)\,=\,y_0 \tag{1}$$
The RK4 algorithm requires a set of approximate values \((k_1, \; k_2, \; k_3, \; \text{and} \; k_4)\) of the slope of the solution curve at different points within the current step over an interval, \(h\).
$$\begin{aligned}k_1 &= h \cdot f \left( x_n, y_n \right) \\ k_2 &= h \cdot f \left( x_n + \frac{h}{2}, y_n+\frac{k_1}{2} \right) \\ k_3 &= h \cdot f \left( x_n + \frac{h}{2}, y_n + \frac{k_2}{2} \right) \\ k_4 &= h \cdot f \left( x_n + h , y_n + k_3 \right) \end{aligned} \tag{2}$$
The \(k\) values from Equation (2) are then weighted and combined using Equation (3) to determine the estimates for \(y_{n + 1}\) of the numerical solution and \(x_{n + 1}\) is updated as per Equation (4).
$$y_{n+1} = y_n + \frac{1}{6} \left( k_1 + 2 k_2 + 2 k_3 + k_4 \right) \tag{3}$$
$$x_{n+1} = x_n + h \tag{4}$$
Let us apply using Scilab scripting the implementation of approximating solution to the ODE with initial condition in Equation (5) using the RK4 solution algorithm.
$$\frac{dy}{dx} = \frac{2x}{y} - xy , \quad y(0) = 1 \tag{5}$$
- Scilab Script
The first block of the script is the
function definition in Scilab of Equation (5).
clear;clc; // (1) 1st-order ODE (dy/dx = f(x, y)) in Scilab form function dydx=f(x, y) dydx = (2*x./y)-(x.*y); // Example: dy/dx = 2x/y - xy endfunction
The RK4 solution algorithm is defined in Scilab
and shown in the second block of script.
// (2) Scilab function of Runge-Kutta 4th order method (RK4) function [x, y]=rk4(f, x0, y0, h, n) x = zeros(1, n+1); y = zeros(1, n+1); x(1) = x0; y(1) = y0; for i = 1:n k1 = h * f(x(i), y(i)); k2 = h * f(x(i) + h/2, y(i) + k1/2); k3 = h * f(x(i) + h/2, y(i) + k2/2); k4 = h * f(x(i) + h, y(i) + k3); x(i+1) = x(i) + h; y(i+1) = y(i) + (k1 + 2*k2 + 2*k3 + k4) / 6; end endfunction
Initial conditions and RK4 parameters
such as step size and number of intervals are shown below.
// (3) initial conditions and parameters x0 = 0; // initial value of x y0 = 1; // initial value of y h = 0.1; // step size n = 50; // no. of interval
Scilab function for the RK4 algorithm is executed
to find the solution vectors for both \(x\) and \(y\) variables. Displaying on the Scilab console
the last element of the variables provides an insight for the domain and range
of the graph.
// (4) numerical solution to the 1st-order ODE IVP using RK4 [x, y] = rk4(f, x0, y0, h, n); disp(x($),y($)); // console display of the last value of x & y
Finally, visualization of numerical solutions requires plotting values of variables using appropriate graph. Shown in Figure 1 is the solution to the given ODE with IVC at \((0,1)\).
// (5) plot of the solution clf; f = gcf(); f.figure_size = [600,600]; plot(x, y,"b-","linewidth",3); xlabel("$\large\textbf{x-value}$"); ylabel("$\large\textbf{y-value}$"); title("$\large\text{Solution to}\;\frac{dy}{dx}=\frac{2x}{y}-{xy}$"); legend("$\LARGE{y\,=\,f(x)}\;\text{with RK4}$",with_box=%F); xstring(2.75,1.05,"https://gani-mech-toolbox.blogspot.com"); xgrid(3,1); ax = gca(); ax.data_bounds = [0 1; 5 1.5];
- Scilab Output (Figure 1)
Figure 1. Solution Plot of 1st-order ODE with IVC using RK4.
RK4 solution to 2nd-order ODE
with IVC shall be presented in our next article.
Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the soft copy of the file.
