Numerical Solutions with Runge-Kutta Method using Scilab

Numerical Solutions with Runge-Kutta Method using Scilab

The Runge-Kutta method is a numerical technique used for solving ordinary differential equations (ODEs). It was developed by the German mathematicians Carl Runge and Martin Kutta in the early 20th century. The method is based on the idea of approximating the solution of an ODE at discrete time steps, using a weighted average of the slope of the solution curve at those time steps.

The basic idea of the Runge-Kutta method is to approximate the solution of an ODE by using a weighted average of the slope of the solution curve at various intermediate points. The method is an iterative process that involves computing the slope of the solution curve at each of these intermediate points, and then combining these slopes to obtain an estimate of the solution at the next time step.

The most commonly used form of the Runge-Kutta method is the fourth-order method, which is also known as RK4. This method involves computing four intermediate points to estimate the solution at the next time step. The formula for the RK4 method is as follows:

k1 = f(t, y)

k2 = f(t + h/2, y + h*k1/2)

k3 = f(t + h/2, y + h*k2/2)

k4 = f(t + h, y + h*k3) 

y(t + h) = y(t) + h*(k1 + 2*k2 + 2*k3 + k4)/6

Getting started with Scilab

Scilab is an open-source, cross-platform numerical computing software that is widely used for scientific and engineering applications. It offers a powerful set of tools for numerical analysis, data visualization, and scientific programming. Before we begin, make sure that Scilab is installed on your computer. You can download it from the official website and follow the installation instructions.

In this formula, t and y represent the current time and the current value of the solution, respectively. f(t, y) represents the derivative of the solution at time t, and h represents the step size. The k variables represent the slopes of the solution curve at various intermediate points, and they are computed using the derivative f(t, y) and the current value of the solution y.

Examlple 1:

y''+4*y'+5*y=0 ,Define the initial conditions 

h = 0.1; x0 = 0; y0 = 1; z0 = -4; xn = 1;

Code:

// Define the differential equations
function f=eq1(x, y, z)
    f = z;
endfunction

function g=eq2(x, y, z)
    g = -4*z - 5*y;
endfunction

// Define the Runge-Kutta 4 method
function [x, y, z]=rk4(f, g, x0, y0, z0, h, xn)
    n = round((xn-x0)/h);  // number of steps
    x = zeros(1,n+1);
    y = zeros(1,n+1);
    z = zeros(1,n+1);

    x(1) = x0;
    y(1) = y0;
    z(1) = z0;

    for i=1:n
        k1 = h * f(x(i), y(i), z(i));
        l1 = h * g(x(i), y(i), z(i));

        k2 = h * f(x(i) + h/2, y(i) + k1/2, z(i) + l1/2);
        l2 = h * g(x(i) + h/2, y(i) + k1/2, z(i) + l1/2);

        k3 = h * f(x(i) + h/2, y(i) + k2/2, z(i) + l2/2);
        l3 = h * g(x(i) + h/2, y(i) + k2/2, z(i) + l2/2);

        k4 = h * f(x(i) + h, y(i) + k3, z(i) + l3);
        l4 = h * g(x(i) + h, y(i) + k3, z(i) + l3);

        y(i+1) = y(i) + (1/6) * (k1 + 2*k2 + 2*k3 + k4);
        z(i+1) = z(i) + (1/6) * (l1 + 2*l2 + 2*l3 + l4);

        x(i+1) = x(i) + h;
    end

endfunction

// Define the initial conditions
x0 = 0;
y0 = 1;
z0 = -4;

// Define the step size and the final value of x
h = 0.1;
xn = 1;

// Call the Runge-Kutta 4 method to solve the differential equation
[x, y, z] = rk4(eq1, eq2, x0, y0, z0, h, xn);

// Print the result
disp("x    y    z")
for i=1:length(x)
    printf("%-4.1f %-4.3f %-4.3f\n", x(i), y(i), z(i))
end
Example 2:
Solve x^2*y''-4x(1+x)y''+2(1+x)y=x^3 with the boundary condition at x=1, y=1/2*e^2, 
y'=-3/2*e^2-0.5; in the range 1<=x<=3. Plot y and y' against x in the given range of the same graph.
Code:
function f=eq1(x, y, z)
    f = z;
endfunction

function g=eq2(x, y, z)
    g = (4*x*(1+x)*z - 2*(1+x)*y + x^3) / x^2;
endfunction
x0 = 1;
y0 = 1/2*exp(2);
z0 = -3/2*exp(2) - 0.5;
h = 0.01;
x = x0:h:3;
y = zeros(1, length(x));
z = zeros(1, length(x));
y(1) = y0;
z(1) = z0;

for i = 1:length(x)-1
    k1 = h*eq1(x(i), y(i), z(i));
    l1 = h*eq2(x(i), y(i), z(i));
    k2 = h*eq1(x(i)+h/2, y(i)+k1/2, z(i)+l1/2);
    l2 = h*eq2(x(i)+h/2, y(i)+k1/2, z(i)+l1/2);
    k3 = h*eq1(x(i)+h/2, y(i)+k2/2, z(i)+l2/2);
    l3 = h*eq2(x(i)+h/2, y(i)+k2/2, z(i)+l2/2);
    k4 = h*eq1(x(i)+h, y(i)+k3, z(i)+l3);
    l4 = h*eq2(x(i)+h, y(i)+k3, z(i)+l3);
    y(i+1) = y(i) + 1/6*(k1 + 2*k2 + 2*k3 + k4);
    z(i+1) = z(i) + 1/6*(l1 + 2*l2 + 2*l3 + l4);
end
plot(x, y, "-b", x, z, "-r");
xlabel("x");
ylabel("y, dy/dx");
title("Solution of x^2*d2y/dx2 - 4*x*(1+x)*dy/dx + 2*(1+x)*y = x^3");
legend(["y", "dy/dx"], "location", "northwest");
// Print the result
disp("x    y    z")
for i=1:length(x)
    printf("%-4.1f %-4.3f %-4.3f\n", x(i), y(i), z(i))
end