First-Order IVP ODE Solution using Runge-Kutta RK4 Method

Topic: Numerical Solution to Ordinary Differential Equation

Subject: Numerical Methods

Tool: Scilab

By: Gani Comia, Jan. 2025

The 4^th -order Runge-Kutta methods (RK4) is a popular and most widely used numerical technique of approximating solution to the initial value problem (IVP) of 1^st -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 2^nd-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.