Help with using the Runge-Kutta 4th order method on a system of 2 first order ODE's.

Question

The original ODE I had was $$ \frac{d^2y}{dx^2}+\frac{dy}{dx}-6y=0$$ with $y(0)=3$ and $y'(0)=1$. Now I can solve this by hand and obtain that $y(1) = 14.82789927$. However I wish to use the 4th order Runge-Kutta method, so I have the system:

$$ \left\{\begin{array}{l} \frac{dy}{dx} = z \\ \frac{dz}{dx} = 6y - z \end{array}\right. $$ With $y(0)=3$ and $z(0)=1$.

Now I know that for two general 1st order ODE's $$ \frac{dy}{dx} = f(x,y,z) \\ \frac{dz}{dx}=g(x,y,z)$$ The 4th order Runge-Kutta formula's for a system of 2 ODE's are: $$ y_{i+1}=y_i + \frac{1}{6}(k_0+2k_1+2k_2+k_3) \\ z_{i+1}=z_i + \frac{1}{6}(l_0+2l_1+2l_2+l_3) $$ Where $$k_0 = hf(x_i,y_i,z_i) \\ k_1 = hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_0,z_i+\frac{1}{2}l_0) \\ k_2 = hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1,z_i+\frac{1}{2}l_1) \\ k_3 = hf(x_i+h,y_i+k_2,z_i+l_2) $$ and $$l_0 = hg(x_i,y_i,z_i) \\ l_1 = hg(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_0,z_i+\frac{1}{2}l_0) \\ l_2 = hg(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1,z_i+\frac{1}{2}l_1) \\ l_3 = hg(x_i+h,y_i+k_2,z_i+l_2)$$

My problem is I am struggling to apply this method to my system of ODE's so that I can program a method that can solve any system of 2 first order ODE's using the formulas above, I would like for someone to please run through one step of the method, so I can understand it better.

Answer 1

I will outline the process and you can fill in the calculations.

We have our system as:

$$ \left\{\begin{array}{l} \frac{dy}{dx} = z \\ \frac{dz}{dx} = 6y - z \end{array}\right. $$ With $y(0)=3$ and $z(0)=1$.

We must do the calculations in a certain order as there are dependencies between the numerical calculations. This order is:

• $k_0 = hf(x_i,y_i,z_i)$

• $l_0 = hg(x_i,y_i,z_i)$

• $k_1 = hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_0,z_i+\frac{1}{2}l_0)$

• $l_1 = hg(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_0,z_i+\frac{1}{2}l_0)$

• $k_2 = hf(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1,z_i+\frac{1}{2}l_1)$

• $l_2 = hg(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1,z_i+\frac{1}{2}l_1)$

• $k_3 = hf(x_i+h,y_i+k_2,z_i+l_2)$

• $l_3 = hg(x_i+h,y_i+k_2,z_i+l_2)$

• $y_{i+1}=y_i + \frac{1}{6}(k_0+2k_1+2k_2+k_3)$

• $z_{i+1}=z_i + \frac{1}{6}(l_0+2l_1+2l_2+l_3)$

We typically need some inputs for the algorithm:

• A range that we want to do the calculations over: $a \le t \le b$, lets use $a = 0, b = 1$.
• The number of steps $N$, say $N = 10$.
• The steps size $h = \dfrac{b-a}{N} = \dfrac{1}{10}$

The system we are solving is:

$$ \frac{dy}{dx} = f(x,y,z) = z \\ \frac{dz}{dx}=g(x,y,z) = 6y - z$$

Doing the calculations using the above order for the first time step $i= 0, t_0 = 0 = x_0$, yields:

• $k_0 = hf(x_0,y_0,z_0) = \dfrac{1}{10}(z_0) = \dfrac{1}{10}(1) = \dfrac{1}{10}$

• $l_0 = hg(x_0,y_0,z_0) = \dfrac{1}{10}(6y_0 - z_0) = \dfrac{1}{10}(6 \times 3 - 1) = \dfrac{1}{10}(17)$

• $k_1 = hf(x_0+\frac{1}{2}h,y_0+\frac{1}{2}k_0,z_0+\frac{1}{2}l_0) = \dfrac{1}{10}(1 + \dfrac{1}{2}\dfrac{1}{10}(17)) ~~$(You please continue the calcs.)

• $l_1 = hg(x_0+\frac{1}{2}h,y_i+\frac{1}{2}k_0,z_0+\frac{1}{2}l_0)$

• $k_2 = hf(x_0+\frac{1}{2}h,y_0+\frac{1}{2}k_1,z_0+\frac{1}{2}l_1)$

• $l_2 = hg(x_0+\frac{1}{2}h,y_0+\frac{1}{2}k_1,z_0+\frac{1}{2}l_1)$

• $k_3 = hf(x_0+h,y_0+k_2,z_0+l_2)$

• $l_3 = hg(x_0+h,y_0+k_2,z_0+l_2)$

• $y_{1}=y_0 + \frac{1}{6}(k_0+2k_1+2k_2+k_3)$

• $z_{1}=z_0 + \frac{1}{6}(l_0+2l_1+2l_2+l_3)$

You now have $x_1$ and $z_1$ which you need for the next time step after all of the intermediate (in order again). Now, we move on to the next time step $i = 1, t_1 = t_0 + h = \dfrac{1}{10} = x_1$, so we have:

• $k_0 = hf(x_1,y_1,z_1) = \dfrac{1}{10}(z_1)$

• $l_0 = hg(x_1,y_1,z_1) = \dfrac{1}{10}(6y_1 - z_1)$

• $k_1 = hf(x_1+\frac{1}{2}h,y_1+\frac{1}{2}k_0,z_1+\frac{1}{2}l_0)$

• $l_1 = hg(x_1+\frac{1}{2}h,y_1+\frac{1}{2}k_0,z_1+\frac{1}{2}l_0)$

• $k_2 = hf(x_1+\frac{1}{2}h,y_1+\frac{1}{2}k_1,z_1+\frac{1}{2}l_1)$

• $l_2 = hg(x_1+\frac{1}{2}h,y_1+\frac{1}{2}k_1,z_1+\frac{1}{2}l_1)$

• $k_3 = hf(x_1+h,y_1+k_2,z_1+l_2)$

• $l_3 = hg(x_1+h,y_1+k_2,z_1+l_2)$

• $y_{2}=y_1 + \frac{1}{6}(k_0+2k_1+2k_2+k_3)$

• $z_{2}=z_1 + \frac{1}{6}(l_0+2l_1+2l_2+l_3)$

Continue this for $10$ time steps. Your final result should match closely (assuming the numerical algorithm is stable for this problem) to the exact solution. You will compare $z_{10}$ to the exact result. The exact solution is:

$$y(x) = e^{-3 x}+2 e^{2 x}$$

If we find $y(1) = \dfrac{1}{e^3} + 2 e^2 = 14.8278992662291643974401973...$.

Answer 2

A Rust implementation is below:

const M: usize = 2; // Size of the system

fn main() {
    let mut y: [f64; M] = [0.0; M];
    let mut a: f64 = 0.0;
    let mut x: f64 = a;
    let mut b: f64 = 0.0;
    let n: usize = 10;
    let h: f64 = (b - a) / n as f64;

    // Initial conditions
    y[0] = 3.0; // y(0)
    y[1] = 1.0; // y'(0)

    for _i in 0..n {
        y = rk4_step(x, y, h, &derivatives);
        x += h;
    }

    println!("{:?}", y);
}

// Define the derivatives for the system of equations
fn derivatives(x: f64, y: [f64; M]) -> [f64; M] {
    let mut dy_dx: [f64; M] = [0.0; M];

    dy_dx[0] = y[1]; // dy/dx = z
    dy_dx[1] = 6.0 * y[0] - y[1]; // dz/dx = 6y - z

    dy_dx
}

// Perform one step of RK4
fn rk4_step(x: f64, y_n: [f64; M], h: f64, f: &dyn Fn(f64, [f64; M]) -> [f64; M]) -> [f64; M] {
    let mut k1: [f64; M] = [0.0; M];
    let mut k2: [f64; M] = [0.0; M];
    let mut k3: [f64; M] = [0.0; M];
    let mut k4: [f64; M] = [0.0; M];
    let mut y_n_plus_1: [f64; M] = [0.0; M];

    k1 = scale_vec(&f(x, y_n), h);
    k2 = scale_vec(&f(x + h / 2.0, add_vec(&y_n, &scale_vec(&k1, 0.5))), h);
    k3 = scale_vec(&f(x + h / 2.0, add_vec(&y_n, &scale_vec(&k2, 0.5))), h);
    k4 = scale_vec(&f(x + h, add_vec(&y_n, &k3)), h);

    y_n_plus_1 = add_vec(
        &y_n,
        &scale_vec(
            &add_vec(&k1, &add_vec(&scale_vec(&k2, 2.0), &add_vec(&scale_vec(&k3, 2.0), &k4))),
            1.0 / 6.0,
        ),
    );

    y_n_plus_1
}

// Helper function: element-wise vector addition
fn add_vec(a: &[f64; M], b: &[f64; M]) -> [f64; M] {
    let mut result: [f64; M] = [0.0; M];
    for i in 0..M {
        result[i] = a[i] + b[i];
    }
    result
}

// Helper function: element-wise scalar multiplication of a vector
fn scale_vec(a: &[f64; M], scalar: f64) -> [f64; M] {
    let mut result: [f64; M] = [0.0; M];
    for i in 0..M {
        result[i] = a[i] * scalar;
    }
    result
}

Answer 3

A Matlab implementation is given below:

% It calculates ODE using Runge-Kutta 4th order method
% Author Ido Schwartz
% Originally available form: http://www.mathworks.com/matlabcentral/fileexchange/29851-runge-kutta-4th-order-ode/content/Runge_Kutta_4.m
% Edited by Amin A. Mohammed, for 2 ODEs(April 2016)

clc;                                               % Clears the screen
clear all;

h=0.1;                                             % step size
x = 0:h:1;                                         % Calculates upto y(1)
y = zeros(1,length(x)); 
z = zeros(1,length(x)); 
y(1) = 3;                                          % initial condition
z(1) = 1;                                          % initial condition
% F_xy = @(t,r) 3.*exp(-t)-0.4*r;                  % change the function as you desire
F_xyz = @(x,y,z) z;                                  % change the function as you desire
G_xyz = @(x,y,z) 6*y-z;

for i=1:(length(x)-1)                              % calculation loop
    k_1 = F_xyz(x(i),y(i),z(i));
    L_1 = G_xyz(x(i),y(i),z(i));
    k_2 = F_xyz(x(i)+0.5*h,y(i)+0.5*h*k_1,z(i)+0.5*h*L_1);
    L_2 = G_xyz(x(i)+0.5*h,y(i)+0.5*h*k_1,z(i)+0.5*h*L_1);
    k_3 = F_xyz((x(i)+0.5*h),(y(i)+0.5*h*k_2),(z(i)+0.5*h*L_2));
    L_3 = G_xyz((x(i)+0.5*h),(y(i)+0.5*h*k_2),(z(i)+0.5*h*L_2));
    k_4 = F_xyz((x(i)+h),(y(i)+k_3*h),(z(i)+L_3*h)); % Corrected        
    L_4 = G_xyz((x(i)+h),(y(i)+k_3*h),(z(i)+L_3*h));

    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  % main equation
    z(i+1) = z(i) + (1/6)*(L_1+2*L_2+2*L_3+L_4)*h;  % main equation

end