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;
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.