Absolute error of ODE45 and Runge-Kutta methods compared with analytical solution

Question

I would appreciate if someone can help with the following issue. I have the following ODE:

dr/dt = 4*exp(0.8*t) - 0.5*r   ,r(0)=2, t[0,1]       (1)

I have solved (1) in two different ways. By means of the Runge-Kutta method (4th order) and by means of ode45 in Matlab. I have compared the both results with the analytic solution, which is given by:

r(t) = 4/1.3 (exp(0.8*t) - exp(-0.5*t)) + 2*exp(-0.5*t)

When I plot the absolute error of each method with respect to the exact solution, I get the following:

For RK-method, my code is:

h=1/50;                                            
x = 0:h:1;                                        
y = zeros(1,length(x)); 
y(1) = 2;    
F_xy = @(t,r) 4.*exp(0.8*t) - 0.5*r;                   
for i=1:(length(x)-1)                              
    k_1 = F_xy(x(i),y(i));
    k_2 = F_xy(x(i)+0.5*h,y(i)+0.5*h*k_1);
    k_3 = F_xy((x(i)+0.5*h),(y(i)+0.5*h*k_2));
    k_4 = F_xy((x(i)+h),(y(i)+k_3*h));
    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  % main equation
end

And for ode45:

tspan = 0:1/50:1;
x0 = 2;
f = @(t,r) 4.*exp(0.8*t) - 0.5*r;
[tid, y_ode45] = ode45(f,tspan,x0);

My question is, why do I have oscillations when I use ode45? (I am referring to the absolute error). Both solutions are accurate (1e-9), but what happens with ode45 in this case?

When I compute the absolute error for the RK-method, why does it looks nicer?

Depside

Answer 1

Your RK4 function is taking fixed steps that are much smaller than those that ode45 is taking. What you're really seeing is the error due to polynomial interpolation that is used to produce the points in between the true steps that ode45 takes. This is often referred to as "dense output" (see Hairer & Ostermann 1990).

When you specify a TSPAN vector with more than two elements, Matlab's ODE suite solvers produce fixed step size output. This does not mean that they that they actually use a fixed step size or that they use the step sizes specified in your TSPAN however. You can see the actual step sizes used and still get your desired fixed step size output by having ode45 output a structure and using deval:

sol = ode45(f,tspan,x0);
diff(sol.x) % Actual step sizes used
y_ode45 = deval(sol,tspan);

You'll see that after an initial step of 0.02, because your ODE is simple it converges to 0.1 for the subsequent steps. The default tolerances combined with the default maximum step size limit (one tenth the integration interval) determine this. Let's plot the error at the true steps:

exactsol = @(t)(4/1.3)*(exp(0.8*t)-exp(-0.5*t))+2*exp(-0.5*t);
abs_err_ode45 = abs(exactsol(tspan)-y_ode45);
abs_err_ode45_true = abs(exactsol(sol.x)-sol.y);
abs_err_rk4 = abs(exactsol(tspan)-y);
figure;
plot(tspan,abs_err_ode45,'b',sol.x,abs_err_ode45_true,'k.',tspan,abs_err_rk4,'r--')
legend('ODE45','ODE45 (True Steps)','RK4',2)

As you can see, the error at the true steps grows more slowly than the error for RK4 (ode45 is effectively a higher order method than RK4 so you'd expect this). The error grows in between the integration points due to the interpolation. If you want to limit this, then you should adjust the tolerances or other options via odeset.

If you wanted to force ode45 to use a step of 1/50 you can do this (works because your ODE is simple):

opts = odeset('MaxStep',1/50,'InitialStep',1/50);
sol = ode45(f,tspan,x0,opts);
diff(sol.x)
y_ode45 = deval(sol,tspan);

For another experiment, try enlarging the integration interval to integrate out to t = 10 maybe. You'll see lots of interesting behavior in the error (plotting relative error is useful here). Can you explain this? Can you use ode45 and odeset to produce results that behave well? Integrating exponential functions over large intervals with adaptive step methods is challenging and ode45 is not necessarily the best tool for the job. There are alternatives however, but they may require some programming.

Edwinaedwine

Hi @horchler. It is a pity I can onlt give one score to your answer. At this moment I am not sure what I am enviyng more: your programming skills or your mathematical understanding. – Depside

I did not know that the actually steps from ode45 where the points from below. If this is the case, I am 100% agree with you, the ode is more accurate. Could you take a 'tic toc' to both methods? I have experienced that the RK-method uses 0.018 seconds while ode45 uses 0.5 seconds. Can I conclude that ode45 is more precise but slower? And, I believed that the adaptive step-size controller from ode45 was very good(?) – Depside

I have not tried to enlarge the integration interval yet, but I am really curious about the output. At the same time, it intrigues me your comments about how good results can I expect trom ode45 and odeset. I will try this as son as I can! Thanks for sharing! – Depside

@SergioHaram: Regrading timing. ode45 requires some initialization. For example, it needs to figure out what step sizes to use. This takes time. You can give it a hint via opts = odeset('InitialStep',0.1); and then passing in opts as the last argument to ode45. Also, because the ODE is so simple (at least it looks simple to the integrator, but as I explain at the end exponential growth can be challenging) you could try using ode23`. Simpler fixed step size methods can be faster in many cases, but usually not when the ODE is more complex, e.g., many oscillators. – Edwinaedwine

@ZheyuanLi: I'd suggest heading to the Matlab Chat Room for any further discussion/questions on this. There are many there who might help you out. – Edwinaedwine

Answer 2

ode45 is coupled rk4-rk5. Personally I think the ODE45 error is nicer. Notice that it stays bounded. The ode4 gets corrected when error magnitude gets too big, and the minimum error per cycle is about 1e-10. The rk4 is "running away" and nothing is stopping it.

Illailladvised

That is not what the error plot shows actually. ode45 converges to a constant step size by the second step in this case. And the plot doesn't show the error being "bounded" either, it's just growing more slowly as you'd expect. – Edwinaedwine

you are right and I gave a positive score on your answer for it. A more accurate way, and the one that I should have used in describing is "running away faster" and "error rate stays more bounded". – Illailladvised