Thread: Comparison of ODE solvers

Hi,

I was reading about Runge-Kutta methods and wanted to do a quick comparison of Euler's method, fourth order Runge-Kutta and MATLAB's ODE45.

I used the following differential equation:
yprime = @(t,y) vanderpoldemo(t,y,1);

Click here for more info on vanderpoldemo.

I was expecting the Runge-Kutta to be closer to ODE45 than Euler, but for some reason it is not. As can be seen from the plots, the Runge-Kutta solution is further away form ODE45 (which I'm assuming is the most correct one) than Euler's.

By changing the constant term in Runge-Kutta from (1/6) to about (1/3.8) I get quite close to ODE45. Here's a copy-paste of Runge-Kutta from MATLAB so you know which constant I'm talking about;

y(:,i+1) = y(:,i)+(1/6)*h*(k1+2*k2+2*k3+k4);

Should I just assume that my implementation of Runge-Kutta is wrong? I can't find any errors.. Here's the implementation:

Code:

function [t,y] = runge_kutta_4order(yprime, tspan, y0, h)

t0 = tspan(1); tfinal = tspan(end); 

% set up the t values at which we will approximate the solution
t = [t0:h:tfinal]';

% include tfinal even if h does not evenly divide tfinal-t0
if t(end)~=tfinal, t=[t tfinal]; end

m = length(t);
y = [y0 zeros(length(y0), length(t)-1)];

    for i=1:(m-1)
        k1 = feval(yprime,t(i),y(:,i));
        k2 = feval(yprime,t(i)+0.5*h, y(:,i)+0.5*h.*k1);
        k3 = feval(yprime,t(i)+0.5*h, y(:,i)*0.5*h.*k2);
        k4 = feval(yprime,t(i)+h, y(:,i)+h.*k3);
    
        y(:,i+1) = y(:,i)+(1/6)*h*(k1+2*k2+2*k3+k4);
    end

All comments are welcome, thanks!

Hi there,

ODE45 is, essentially, an RK4 method, and so you should see almost no error between the results that ode45 give, and the results that RK4 give.

I just ran my own code, with my own rk4 method, and the from rk4 and ode45 are identical, almost no difference whatsoever.

With that said, I can't see any particular problem with your implementation.

So, what I'll do is give you my code, and allow you to spot any errors yourself.

It could be that your stepsize is too large? What stepsize are you using?

Here is my code. I haven't commented it, but hopefully it is self-explanatory.

Code:

clear
clc
StartTime = 0;
EndTime = 20;
Stepsize = 0.01;
x = [2,2];
xdot = [0,0];
tout=zeros(1,(EndTime/Stepsize)+1);
xout=zeros((EndTime/Stepsize)+1, length(x));
xdout = xout;
i=0;

for time = StartTime:Stepsize:EndTime,
    
        i          = i+1  ;			
        tout(i)    = time ;	 	    
        xout(i,:)  = x    ;		    
        xdout(i,:) = xdot ;	        
    
    xdot = vanderpoldemo(time,x,1);         
    x = rk4('vanderpoldemo',Stepsize,time,x); 
end

[TOUT,YOUT] = ode45(@vanderpoldemo,[StartTime EndTime],[2;2]);

figure(1)
hold on
plot(tout,xout(:,1),'r--')
plot(TOUT,YOUT(:,1),'b--')
title('Comparison of ODE Solvers')
xlabel('time (s)')
ylabel('Function value')
legend('Runge-Kutta 4','ODE45')

figure(2)
hold on
plot(tout,xout(:,2),'r--')
plot(TOUT,YOUT(:,2),'b--')
title('Comparison of ODE Solvers')
xlabel('time (s)')
ylabel('Function value')
legend('Runge-Kutta 4','ODE45')

Code:

function x = rk4(model,h,t,x)
k1 = feval(model,t,x,1)        ;  
k2 = feval(model,t+0.5*h,x+h*k1'/2,1) ;   
k3 = feval(model,t+0.5*h,x+h*k2'/2,1) ;   
k4 = feval(model,t+h,x+h*k3',1)   ;  
x = x + h*(k1' + 2*k2' + 2*k3' + k4')/6    ;   
end

Code:

function dydt = vanderpoldemo(t,y,Mu)
Mu = 1;
dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];

Let me know how you get on.

Thank you so much. After looking at your code I noticed that I had a * instead of a + at one place in my implementation of rk4. I must be going blind.
Again, thanks!