# Comparison of ODE solvers

• Mar 25th 2011, 07:05 AM
Mollier
[SOLVED] 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);

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

http://i163.photobucket.com/albums/t...ison_sys_1.jpg

http://i163.photobucket.com/albums/t...ison_sys_2.jpg

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!
• Mar 26th 2011, 10:32 AM
Phugoid
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.
• Mar 26th 2011, 11:10 PM
Mollier
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!