MATLAB - converting Euler method to Heun's method

Question

I saw this Euler method that calculates errors:

function [t,le,ge] = euler_errors(h)

f=@(u) u*(2-u); % this is the function for the IVP
t0=0;
tn=1;
t=t0:h:tn;%we want to find the errors along this solutions
%here is the exact solution of the IVP
u_exact=(0.2*exp(2*t))./(2+0.1*(exp(2*t)+1)); %the with initial value u(0)=0.1

n=length(t);
u_e=zeros(1,n);
u_g=zeros(1,n);

i=1;
u_e(i)=0.1;
u_g(i)=0.1;

%u_e and u_g are both values given by Euler method
%u_e is for the local error
%u_g is for the global error
while (i<n)
    u_e(i+1)=u_e(i)+h*f(u_e(i));
    u_g(i+1)=u_g(i)+h*f(u_exact(i));
    i=i+1;
end;
%le1 is the local error
%ge1 is the global error
le=abs(u_e-u_exact);
ge=abs(u_g-u_exact);
end

I tried to convert the method into Heun's method – here is my attempt:

function [t,le,ge] = heun_errors(h)

f=@(u) u*(2-u); % this is the function for the IVP
t0=0;
tn=1;
t=t0:h:tn;%we want to find the errors along this solutions
%here is the exact solution of the IVP
u_exact=(0.2*exp(2*t))./(2+0.1*(exp(2*t)+1)); %the with initial value u(0)=0.1

n=length(t);
u_e=zeros(1,n);
u_g=zeros(1,n);

i=1;
u_e(i)=0.1;
u_g(i)=0.1;

%u_e and u_g are both values given by Euler method
%u_e is for the local error
%u_g is for the global error
while (i<n)
    u_e1(i+1)=u_e(i)+h*f(u_e(i));
u_e(i+1)=u_e(i)+(h/2)*(f(u_e(i))+f(u_e1(i+1)));
u_g1(i+1)=u_g(i)+h*f(u_exact(i));
u_g(i+1)=u_g(i)+(h/2)*(f(u_exact(i))+f(u_g1(i+1)));

    i=i+1;
end;
%le1 is the local error
%ge1 is the global error
le=abs(u_e-u_exact);
ge=abs(u_g-u_exact);
end

But now the error has actually increased. Can someone tell me what I did wrong and how I might fix it?

Answer 1

Correction of the Exact Solution

First, the exact solution has a sign error which invalidates all error computations. As Bernoulli equation the solution can be obtained via

(e^(2t)/u)'=e^(2t)*(-u'/u^2 + 2/u) = e^(2t)

e^(2t)/u(t)-1/u0 = (e^(2t)-1)/2  

u(t) = 2*u0*exp(2*t) ./ ( 2 + u0*(exp(2*t)-1) )

---

On the Computation of Local and Global Errors

The organization of the Euler code for the computation of local and global error already has something conceptually mixed up. The full, normal integration gives the global result and thus also the global error in difference to the exact solution. Thus what is called u_e should be u_g.

The local error compares the exact solution u_i(t) started in (t_i,u_i) against the method step from the same point. As there is only one exact solution in play, it should perform the step for (t_i,u_i) = (t(i), u_exact(i)) to then compare the method step whose result gets stored in u_e(i+1) against the exact value u_exact(i+1).

Both issued are addressed with the change of the loop to

while (i<n)
  u_g(i+1)=u_g(i)+h*f(u_g(i));
  u_e(i+1)=u_exact(i)+h*f(u_exact(i));
  i=i+1;
end;

Similarly, these issues also have to be addressed in the Heun code, where the loop thus changes to (using intermediary variables for the stages as in a general Runge-Kutta methos)

while (i<n)
  k1=f(u_g(i));
  k2=f(u_g(i)+h*k1);
  u_g(i+1)=u_g(i)+(h/2)*(k1+k2);
  k1=f(u_exact(i));
  k2=f(u_exact(i)+h*k1);
  u_e(i+1)=u_exact(i)+(h/2)*(k1+k2);
  i=i+1;
end;

Plots of the Error Profiles

Then plotting the error profiles le/h^(p+1) of the local and lg/h^p global error over time, with the error orders p=1 for Euler and p=2 for Heun, and overlaying the dot plots of the different step sizes, gives

which are visually invariant relative to the step size.

As a check for consistency, the local Euler error is around 0.3*h^2, and the Lipschitz constant of the ODE around 2 in the relevant region, so the global error follows by the Grönwall lemma approximately the law (exp(2*t)-1)*0.15*h for small t. Similarly, the Heun method has the local error around 0.15*h^3 giving a global error (exp(2*t)-1)*0.075*h^2 for small t, which is both compatible with the global error graph.

For larger t, the sign of the local error can change, thus compensating for a part of the earlier global error. This however is not guaranteed to happen and also not easily predictable.