Implementing Runge-Kutta for Systems of DE

Question

I'm trying to implement the Runge-Kutta Method for Systems of DEs in MATLAB. I'm not getting the correct answers, I'm not sure if there is something wrong in the code or the commands I use to run it.

Here's my code:

function RKSystems(a, b, m, N, alpha, f)
    h = (b - a)/N;
    t = a;
    w = zeros(1, m);

    for j = 1:m
        w(j) = alpha(j);
    end

    fprintf('t = %.2f;', t);
    for i = 1:m
        fprintf(' w(%d) = %.10f;', i, w(i));
    end
    fprintf('\n');

    k = zeros(4, m);
    for i = 1:N
        for j = 1:m
           k(1, j) = h*f{j}(t, w);
        end

        for j = 1:m
            k(2, j) = h*f{j}(t + h/2, w + (1/2)*k(1));
        end

        for j = 1:m
            k(3, j) = h*f{j}(t + h/2, w + (1/2)*k(2));
        end

        for j = 1:m
            k(4, j) = h*f{j}(t + h, w + k(3));
        end

        for j = 1:m
            w(j) = w(j) + (k(1, j) + 2*k(2, j) + 2*k(3, j) + k(4, j))/6;
        end

        t = a + i*h;

        fprintf('t = %.2f;', t);
        for k = 1:m
            fprintf(' w(%d) = %.10f;', k, w(k));
        end
        fprintf('\n');

    end 
end

I'm trying to test it out on this problem. Here are my commands and output:

>> U1 = @(t, u) 3*u(1) + 2*u(2) - (2*t^2 + 1)*exp(2*t);

>> U2 = @(t, u) 4*u(1) + u(2) + (t^2 + 2*t - 4)*exp(2*t);

>> a = 0; b = 1; alpha = [1 1]; m = 2; h = 0.2; N = (b - a)/h;

>> RKSystems(a, b, m, N, alpha, {U1 U2});

t = 0.00; w(1) = 1.0000000000; w(2) = 1.0000000000;

t = 0.20; w(1) = 2.2930309680; w(2) = 1.6186020410;

t = 0.40; w(1) = 5.0379629523; w(2) = 3.7300162165;

t = 0.60; w(1) = 11.4076339762; w(2) = 9.7009491301;

t = 0.80; w(1) = 27.0898576892; w(2) = 25.6603894354;

t = 1.00; w(1) = 67.1832886708; w(2) = 67.6103125539;

Answer 1

My problems were in the following lines of code:

k(2, j) = h*f{j}(t + h/2, w + (1/2)*k(1));

k(3, j) = h*f{j}(t + h/2, w + (1/2)*k(2));

k(4, j) = h*f{j}(t + h, w + k(3));

I was expecting k(1) to add the entire first row of k (a 4 by m matrix) to the matrix w (a 1 by m matrix), but it was only adding the first element of the first row. To fix it, I modified the lines as follows:

k(2, j) = h*f{j}(t + h/2, w + (1/2)*k(1, :));

k(3, j) = h*f{j}(t + h/2, w + (1/2)*k(2, :));

k(4, j) = h*f{j}(t + h, w + k(3, :));

Answer 2

So... here is how I would do it, it is very difficult for me to read what is going on in your code snippet, but I hope this helps you out. Additionally, matlab has built in rk solvers I would suggest becoming familiar with those.

%rk4 solver
dt = .2;
t = 0:dt:1;
u = zeros(2,numel(t));
u(:,1) = 1;

for jj = 2:numel(t)
    u_ = u(:,jj-1);
    t_ = t(jj-1);
    fa = rhs(u_,t_);
    fb = rhs(u_+dt/2.*fa,t_+dt/2);
    fc = rhs(u_+dt/2.*fb,t_+dt/2);
    fd = rhs(u_+dt.*fc,t_+dt);
    u(:,jj) = u(:,jj-1) + dt/6*(fa+2*fb+2*fc+fd);
end
disp([t;u]')

and rhs.m is the following:

function dudt = rhs(u,t)
dudt = [(3*u(1) + 2*u(2) - (2*t^2 + 1)*exp(2*t));
        (4*u(1) + u(2) + (t^2 + 2*t - 4)*exp(2*t))];

This returns the following:

>> rkorder4
     0    1.0000    1.0000
0.2000    2.1204    1.5070
0.4000    4.4412    3.2422
0.6000    9.7391    8.1634
0.8000   22.6766   21.3435
1.0000   55.6612   56.0305

Alternatively, you could invoke ode45 with something like:

dt = .2;
t = 0:dt:1;
rhs=@(t,u)[(3*u(1) + 2*u(2) - (2*t^2 + 1)*exp(2*t));
        (4*u(1) + u(2) + (t^2 + 2*t - 4)*exp(2*t))];

[ts,us]=ode45(@(t,u) rhs(t,u),t,[1 1],[]);
disp([ts us]);

Which gives you:

                   0   1.000000000000000   1.000000000000000
   0.200000000000000   2.125018862140859   1.511597928697035
   0.400000000000000   4.465156492097592   3.266022178547346
   0.600000000000000   9.832481410895053   8.256418221678395
   0.800000000000000  23.003033457636558  21.669270713784108
   1.000000000000000  56.738351357036301  57.106230777717208

Which is pretty close to what you get from my code. The agreement between the two can be increased by decreasing the time step, dt. They will always be in greatest agreement at the low values of t with the difference between the two increasing as t increases. Both implementations are also in pretty close agreement with the analytic solutions.