Solve a system of equations with Runge Kutta 4: Matlab

Question

I want to solve a system of THREE differential equations with the Runge Kutta 4 method in Matlab (Ode45 is not permitted).

After a long time spent looking, all I have been able to find online are either unintelligible examples or general explanations that do not include examples at all. I would like a concrete example on how to implement my solution properly, or the solution to a comparable problem which I can build on.

I have come quite far; my current code spits out a matrix with 2 correct decimals on most of the components, which I am quite happy with.

However, when the step-size is decreased, the errors become enormous. I know the for-loop I have created is not entirely correct. I may have defined the functions incorrectly, but I am quite certain that the problem is solved if some minor changes are made to the for-loop because it appears to be solving the equation-system fairly well already in its current state.

clear all, close all, clc

%{
____________________TASK:______________________
Solve the system of differential equations below 
in the interval 0<t<1, with stepsize h = 0.1.
x'= y                 x(0)=1
y'= -x-2e^t+1         y(0)=0   ,   where x=x(t), y=y(t),  z=z(t)  
z'= -x - e^t + 1      z(0)=1

THE EXACT SOLUTIONS for x y and z can be found in this pdf:
archives.math.utk.edu/ICTCM/VOL16/C029/paper.pdf
_______________________________________________

%}

h = 0.1;
t    = 0:h:1
N = length(t);

%Defining the functions
x    = zeros(N,1);%I am not entierly sure if x y z are supposed to be defined in this way.
y    = zeros(N,1)
z    = zeros(N,1)

f = @(t, x, y, z) -x-2*exp(t)+1;%Question: Do i need a function for x here as well??
g = @(t, x, y, z) -x - exp(t) + 1;

%Starting conditions
x(1) = 1; 
y(1) = 0;
z(1) = 1;

for i = 1:(N-1)
    K1     = h * ( y(i));%____I think z(i) is supposed to be here, but i dont know in what way. 
    L1     = h * f( t(i)          , x(i)        , y(i) ,      z(i));
    M1     = h * g( t(i)          , x(i)        , y(i) ,      z(i));

    K2     = h *  (y(i) + 1/2*L1 + 1/2*M1);%____Again, z(i) should probably be here somewhere. 
    L2     = h * f(t(i) + 1/2*h,  x(i)+1/2*K1 , y(i)+1/2*L1 , z(i)+1/2*M1);
    M2     = h * g(t(i) + 1/2*h,  x(i)+1/2*K1 , y(i)+1/2*L1 , z(i)+1/2*M1);

    K3     = h *  (y(i) + 1/2*L2 + 1/2*M2);%____z(i). Should it just be added, like "+z(i)" ? 
    L3     = h * f(t(i) + 1/2*h,  x(i) + 1/2*K2 , y(i) + 1/2*L2 , z(i) + 1/2*M2);
    M3     = h * g(t(i) + 1/2*h,  x(i) + 1/2*K2 , y(i) + 1/2*L2 , z(i) + 1/2*M2);

    K4     = h *  (y(i) + L3 + M3);%_____z(i) ... ? 
    L4     = h * f( t(i)+h    ,  x(i)+K3     , y(i)+L3,     z(i)+M3);
    M4     = h * g( t(i)+h    ,  x(i)+K3     , y(i)+L3,     z(i)+M3);

    x(i+1) = x(i)+1/6*(K1+2*K2+2*K3+K4);                                                                             
    y(i+1) = y(i)+1/6*(L1+2*L2+2*L3+L4);
    z(i+1) = z(i)+1/6*(M1+2*M2+2*M3+M4);
end

Answer_Matrix = [t' x y z]

your functions `f` and `g` have to have `t,x,y,z` passed to them, but are actually only functions of `x` and `t`... are those functions correct? — Wolfie, Apr 14 '17 at 09:39
Yes, f and g are only functions of x and t in this particular example, but my understanding is that if i make them functions of y and z as well, this code can be applied to any equation system. Love your edit by the way :) — gelbrekt, Apr 14 '17 at 09:43
Again, I am not sure if the functions are correctly defined, but the lines "-x-2*exp(t)+1" and "-x - exp(t) + 1" correctly reflect the given functions y' and z' in the task. — gelbrekt, Apr 14 '17 at 09:45
I think you're missing information about `x`, your derivative functions are governed by it, but you are treating it like a dependent variable? You need a function defining its derivative or to just define it as a constant vector. — Wolfie, Apr 14 '17 at 09:53
You were right, all i needed was a function for x' that says x'=y. How do i mark this as post solved? Also, I think it would be very useful to post the correct code here, but i dont know how to do that... — gelbrekt, Apr 14 '17 at 10:07
The `1/2*M1`, `1/2*M2`, and `M3` are not correct. They shouldn't be included since only `y` is stepped in `x`'s equations. — TroyHaskin, Apr 14 '17 at 10:11
Well, the issue really was that i felt that concrete examples on how to use Rk4 to solve differential equations with three equations were severely lacking on the web. Again; I only managed to find really bad concrete examples and general descriptions. — gelbrekt, Apr 14 '17 at 10:11
Not a problem. In truth, I blame the text for giving such a terrible example of how to code an RK method. Treating the equations individually instead of like a vector makes the code extremely ugly, complex, and error prone. But I'm glad you figured it out. — TroyHaskin, Apr 14 '17 at 10:17
Well, I don't agree that it makes the code complex. I found this to be much more intuitive than the vector functions. Could you recommend a good page/pdf for solving a similar problem with vectors instead (that contains concrete examples)? — gelbrekt, Apr 14 '17 at 10:26
I have given an answer which implements the method in a vectorised form, see below. — Wolfie, Apr 14 '17 at 11:13
Thank you very much Wolfie; this gave me way more help than I had anticipated. Just one more question though; how does Troy's code work exactly? Where is the solution in that code, and what does the W-vector mean (w does not appear to contain the system's solution)? — gelbrekt, Apr 14 '17 at 11:31
@gelbrekt I don't know how his code works because I can't access pastebin from this PC, he should post it as an answer here if he wants it to be discussed here. If my code solves your problem and provides helpful insight, please consider marking it as accepted, thanks. — Wolfie, Apr 14 '17 at 11:48

Wolfie · Accepted Answer · 2017-04-14T11:19:26.480

So your main issue was not defining x properly. You were propagating its value using the Runge Kutta 4 (RK4) method, but never actually defined what its derivative was!

At the bottom of this answer is a function which can take any given number of equations and their initial conditions. This has been included to address your need for a clear example for three (or more) equations.

For reference, the equations can be directly lifted from the standard RK4 method described here.

Working Script

This is comparable to yours, but uses slightly clearer naming conventions and structure.

% Initialise step-size variables
h = 0.1;
t = (0:h:1)';
N = length(t);

% Initialise vectors
x = zeros(N,1);    y = zeros(N,1);    z = zeros(N,1);
% Starting conditions
x(1) = 1;     y(1) = 0;    z(1) = 1;

% Initialise derivative functions
dx = @(t, x, y, z) y;                  % dx = x' = dx/dt
dy = @(t, x, y, z) - x -2*exp(t) + 1;  % dy = y' = dy/dt
dz = @(t, x, y, z) - x -  exp(t) + 1;  % dz = z' = dz/dt

% Initialise K vectors
kx = zeros(1,4); % to store K values for x
ky = zeros(1,4); % to store K values for y
kz = zeros(1,4); % to store K values for z
b = [1 2 2 1];   % RK4 coefficients

% Iterate, computing each K value in turn, then the i+1 step values
for i = 1:(N-1)        
    kx(1) = dx(t(i), x(i), y(i), z(i));
    ky(1) = dy(t(i), x(i), y(i), z(i));
    kz(1) = dz(t(i), x(i), y(i), z(i));

    kx(2) = dx(t(i) + (h/2), x(i) + (h/2)*kx(1), y(i) + (h/2)*ky(1), z(i) + (h/2)*kz(1));
    ky(2) = dy(t(i) + (h/2), x(i) + (h/2)*kx(1), y(i) + (h/2)*ky(1), z(i) + (h/2)*kz(1));
    kz(2) = dz(t(i) + (h/2), x(i) + (h/2)*kx(1), y(i) + (h/2)*ky(1), z(i) + (h/2)*kz(1));

    kx(3) = dx(t(i) + (h/2), x(i) + (h/2)*kx(2), y(i) + (h/2)*ky(2), z(i) + (h/2)*kz(2));
    ky(3) = dy(t(i) + (h/2), x(i) + (h/2)*kx(2), y(i) + (h/2)*ky(2), z(i) + (h/2)*kz(2));
    kz(3) = dz(t(i) + (h/2), x(i) + (h/2)*kx(2), y(i) + (h/2)*ky(2), z(i) + (h/2)*kz(2));

    kx(4) = dx(t(i) + h, x(i) + h*kx(3), y(i) + h*ky(3), z(i) + h*kz(3));
    ky(4) = dy(t(i) + h, x(i) + h*kx(3), y(i) + h*ky(3), z(i) + h*kz(3));
    kz(4) = dz(t(i) + h, x(i) + h*kx(3), y(i) + h*ky(3), z(i) + h*kz(3));

    x(i+1) = x(i) + (h/6)*sum(b.*kx);       
    y(i+1) = y(i) + (h/6)*sum(b.*ky);       
    z(i+1) = z(i) + (h/6)*sum(b.*kz);        
end    

% Group together in one solution matrix
txyz = [t,x,y,z];

Implemented as function

You wanted code which can "be applied to any equation system". To make your script more usable, let's take advantage of vector inputs, where each variable is on its own row, and then make it into a function. The result is something comparable (in how it is called) to Matlab's own ode45.

% setup
odefun = @(t, y) [y(2); -y(1) - 2*exp(t) + 1; -y(1) - exp(t) + 1];
y0 = [1;0;1];
% ODE45 solution
[T, Y] = ode45(odefun, [0,1], y0);
% Custom RK4 solution
t = 0:0.1:1;
y = RK4(odefun, t, y0);
% Compare results
figure; hold on;
plot(T, Y); plot(t, y, '--', 'linewidth', 2)

You can see that the RK4 function (below) gives the same result of the ode45 function.

The function RK4 is simply a "condensed" version of the above script, it will work for however many equations you want to use. For broad use, you would want to include input-checking in the function. I have left this out for clarity.

function y = RK4(odefun, tspan, y0)
% ODEFUN contains the ode functions of the system
% TSPAN  is a 1D vector of equally spaced t values
% Y0     contains the intial conditions for the system variables

    % Initialise step-size variables
    t = tspan(:); % ensure column vector = (0:h:1)';
    h = t(2)-t(1);% define h from t
    N = length(t);

    % Initialise y vector, with a column for each equation in odefun
    y = zeros(N, numel(y0));
    % Starting conditions
    y(1, :) = y0(:)';  % Set intial conditions using row vector of y0

    k = zeros(4, numel(y0));              % Initialise K vectors
    b = repmat([1 2 2 1]', 1, numel(y0)); % RK4 coefficients

    % Iterate, computing each K value in turn, then the i+1 step values
    for i = 1:(N-1)        
        k(1, :) = odefun(t(i), y(i,:));        
        k(2, :) = odefun(t(i) + (h/2), y(i,:) + (h/2)*k(1,:));        
        k(3, :) = odefun(t(i) + (h/2), y(i,:) + (h/2)*k(2,:));        
        k(4, :) = odefun(t(i) + h, y(i,:) + h*k(3,:));

        y(i+1, :) = y(i, :) + (h/6)*sum(b.*k);    
    end    
end

score 1 · Answer 2 · answered Apr 14 '17 at 10:17

Ok, turns out it was just a minor mistake where the x-variable was not defined as a function of y (as x'(t)=y according to the problem.

So: Below is a concrete example on how to solve a differential equation system using Runge Kutta 4 in matlab:

clear all, close all, clc

%{
____________________TASK:______________________
Solve the system of differential equations below 
in the interval 0<t<1, with stepsize h = 0.1.
x'= y                 x(0)=1
y'= -x-2e^t+1         y(0)=0   ,   where x=x(t), y=y(t),  z=z(t)  
z'= -x - e^t + 1      z(0)=1

THE EXACT SOLUTIONS for x y and z can be found in this pdf:
archives.math.utk.edu/ICTCM/VOL16/C029/paper.pdf
_______________________________________________

%}

%Step-size
h = 0.1;
t    = 0:h:1
N = length(t);

%Defining the vectors where the answer is stored. 
x    = zeros(N,1);
y    = zeros(N,1)
z    = zeros(N,1)

%Defining the functions
e = @(t, x, y, z) y;
f = @(t, x, y, z) -x-2*exp(t)+1;
g = @(t, x, y, z) -x - exp(t) + 1;

%Starting/initial conditions
x(1) = 1; 
y(1) = 0;
z(1) = 1;

for i = 1:(N-1)
    K1     = h * e( t(i)          , x(i)        , y(i) ,      z(i));
    L1     = h * f( t(i)          , x(i)        , y(i) ,      z(i));
    M1     = h * g( t(i)          , x(i)        , y(i) ,      z(i));

    K2     = h * e(t(i) + 1/2*h,  x(i)+1/2*K1 , y(i)+1/2*L1 , z(i)+1/2*M1);
    L2     = h * f(t(i) + 1/2*h,  x(i)+1/2*K1 , y(i)+1/2*L1 , z(i)+1/2*M1);
    M2     = h * g(t(i) + 1/2*h,  x(i)+1/2*K1 , y(i)+1/2*L1 , z(i)+1/2*M1);

    K3     = h * e(t(i) + 1/2*h,  x(i) + 1/2*K2 , y(i) + 1/2*L2 , z(i) + 1/2*M2);
    L3     = h * f(t(i) + 1/2*h,  x(i) + 1/2*K2 , y(i) + 1/2*L2 , z(i) + 1/2*M2);
    M3     = h * g(t(i) + 1/2*h,  x(i) + 1/2*K2 , y(i) + 1/2*L2 , z(i) + 1/2*M2);

    K4     = h * e( t(i)+h    ,  x(i)+K3     , y(i)+L3,     z(i)+M3);
    L4     = h * f( t(i)+h    ,  x(i)+K3     , y(i)+L3,     z(i)+M3);
    M4     = h * g( t(i)+h    ,  x(i)+K3     , y(i)+L3,     z(i)+M3);

    x(i+1) = x(i)+1/6*(K1+2*K2+2*K3+K4);                                                                             
    y(i+1) = y(i)+1/6*(L1+2*L2+2*L3+L4);
    z(i+1) = z(i)+1/6*(M1+2*M2+2*M3+M4);
end

Answer_Matrix = [t' x y z]

I'd like to point out that `e` doesn't need to be defined for this to work; using raw values of the solution vector works too. However, having `e` defined does significantly lower the possibility of an error since the actual functional form is hidden beneath a function call. — TroyHaskin, Apr 14 '17 at 10:26

Solve a system of equations with Runge Kutta 4: Matlab

2 Answers2

Working Script

Implemented as function

Linked

Related