Using matlab for an numerical error analysis problem in ODE

Question

(a) Consider the following differential Equation

$$Y'(t)=\frac{1}{1+t^{2}}-2[Y(t)]^2$$ $$Y(0)=0$$

The exact solution is $$Y(t)=\frac{t}{1+t^2}$$

Using the Euler method to solve the following differential equation with the help of matlab (code given below), run the program for $t=10$ and $h=0.1$ compute the error.Then run it for $h=0.05$ compute the error and note how it changes.Then use Richardson Extrapolation to improve the accuracy. Compute the error.

function y=Euler(t,h,y0,t0)

  n=floor((t-t0)/h);
  y=y0;
  for i=1:n
      y=y+h*f(t0+h*(i-1),y);
  end;


function f=f(t,y)

f=1/(1+t^2)-2*y^2;
return;

(b) Now using the Runge-Kutta method to solve the differential equation in part (a) Run the program for $t=10$ and $h=0.1$ The error should be reduced as compared to part(a). Explain the difference. Run it again for $h=0.05$ and notice how the error is changed. Use Richardson extrapolation to improve the accuracy and compute the error.

function y=Runge_Kutta(t,h,y0,t0)

  n=floor((t-t0)/h);
  y=y0;
  for i=1:n
          t=t0+h*(i-1);
          y=y+(h/2)*(f(t,y)+f(t+h,y+h*f(t,y)));
  end;


function f=f(t,y)

f=1/(1+t^2)-2*y^2;
return;

What i tried

(a) When i ran the program with a step size of $0.1$ i got an error of $0.0986$ but when i ran it with a step size of $0.05$, i got an error of $0.0988$ which dosent make sense to me beacuse a smaller step size should give a smaller error but it does not appear so in this case.(There is nothing wrong with the code because it is given by my prof). Could anyone explain this to me. Thanks

(b) Again when i run the program for this part i got an error of $0.0990$ for both step size $0.1$ and $0.05$ which is larger than that of part(a) which dosent make sense to me beacuse the Runge-Kutta method is supposed to be more accurate than the Euler method and hence have smaller error, but it dosent appear so in this case. Could anyone explain to me what is going wrong here. Thanks

Answer 1

You forgot, in your error calculation, to actually subtract the exact value 0.09900990099009901, you should get something like this for the values and errors of the Euler method:

0.1  -> val = 0.09861080536664182,  err = -0.00039909562345719074
0.05 -> val = 0.09881012084504534,  err = -0.00019978014505367403
Rich -> val = 0.09900943632344886,  err = -4.6466665015731934e-07

which nicely demonstrates the O(h) global error.

The values and errors for the trapezoidal method are in my re-implementation

0.1  -> val = 0.09901920896465118,  err =  9.307974552161258e-06
0.05 -> val = 0.09901220700730341,  err =  2.3060172043981586e-06
Rich -> val = 0.09900987302152083,  err = -2.7968578189541127e-08

which is in accord with an error O(h^2), as the second error is about 1/4 of the first.