Results 1 to 3 of 3

Math Help - Comparison of ODE solvers

  1. #1
    Member Mollier's Avatar
    Joined
    Nov 2009
    From
    Norway
    Posts
    234
    Awards
    1

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

    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!
    Last edited by Mollier; March 30th 2011 at 07:44 PM.
    Follow Math Help Forum on Facebook and Google+

  2. #2
    Junior Member
    Joined
    Jul 2010
    Posts
    41
    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.
    Follow Math Help Forum on Facebook and Google+

  3. #3
    Member Mollier's Avatar
    Joined
    Nov 2009
    From
    Norway
    Posts
    234
    Awards
    1
    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!
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. [SOLVED] Comparison test and Limit Comparison test for series
    Posted in the Calculus Forum
    Replies: 5
    Last Post: November 25th 2010, 12:54 AM
  2. Comparison or Limit Comparison Test Problem
    Posted in the Calculus Forum
    Replies: 2
    Last Post: March 12th 2010, 07:46 AM
  3. Limit comparison/comparison test series
    Posted in the Calculus Forum
    Replies: 2
    Last Post: March 25th 2009, 08:27 PM
  4. Comparison & Limit Comparison test for series
    Posted in the Calculus Forum
    Replies: 4
    Last Post: March 25th 2009, 04:00 PM
  5. Story promlem solvers please help
    Posted in the Math Topics Forum
    Replies: 1
    Last Post: November 29th 2006, 12:13 PM

Search Tags


/mathhelpforum @mathhelpforum