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 que…

Accepted Answer

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.

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

Solution

Problem