# Solution to the first MATLAB assignment

First we set some parameters that will be used over and over:
tStart=0;
tEnd=20;
y0=1/2;

## First, let's look with a small number of points

We compare the three different methods with a big stepsize .
N=40;
h=(tEnd-tStart)/N;
t=linspace(tStart,tEnd,N+1);
yEuler=zeros(1,N+1);
yEuler(1)=y0;
yEulerImproved=yEuler;
yRunge=yEuler;
yX=yExact(t);
for k=1:N
yEuler(k+1)=eulerStep(@f,t(k),yEuler(k),h);
yEulerImproved(k+1)=eulerImprovedStep(@f,t(k),yEulerImproved(k),h);
yRunge(k+1)=RKStep(@f,t(k),yRunge(k),h);
end
plot(t,yEuler,t,yEulerImproved,t,yRunge,t,yX,'k--')
xlabel('$t$')
ylabel('$y$')
legend('Euler','Improved','Runke-Kutta','Exact','Location',"eastoutside")

## Plot the errors of the three methods

plot(t,yX-yEuler,t,yX-yEulerImproved,t,yX-yRunge)
xlabel('$t$')
ylabel('Error')
legend('Euler','Improved Euler','Runge-Kutta','location','northwest')
errorEuler=abs(yX(end)-yEuler(end));
errorImproved=abs(yX(end)-yEulerImproved(end));
errorRunge=abs(yX(end)-yRunge(end));
fprintf('With 40 points, the final-time error is:\n')
With 40 points, the final-time error is:
fprintf('Euler method: %0.2d',errorEuler)
Euler method: 2.70e-01
fprintf('Improved Euler method: %0.2d',errorImproved)
Improved Euler method: 6.34e-03
fprintf('Runge-Kutta method: %0.2d',errorRunge)
Runge-Kutta method: 3.27e-05

## Experimentally find h small enough to make error below

We just keep cuttingh in half (doubling N) until we get below the threshold.
while errorEuler>0.01
N=2*N;
h=(tEnd-tStart)/N;
t=linspace(tStart,tEnd,N+1);
y=zeros(1,N+1);
y(1)=y0;
for k=1:N
y(k+1)=eulerStep(@f,t(k),y(k),h);
end
errorEuler=abs(yX(end)-y(end));
fprintf('N = %i, error=%0.3f\n',N,errorEuler)
end
N = 80, error=0.189 N = 160, error=0.110 N = 320, error=0.058 N = 640, error=0.030 N = 1280, error=0.015 N = 2560, error=0.007
plot(t,y,t,yExact(t),'--')
xlabel('$t$');ylabel('$y$')
legend('Euler','Exact')
That's good enough. Somewhere between 1280 and 2560 points, Euler's method returns an error below 0.01.

## For the improved Euler method, I found points to be enough

Here I notice that was good enough in the initial experiment, so I just reduced it a bit at a time by hand until it failed.
N=36;
h=(tEnd-tStart)/N;
t=linspace(tStart,tEnd,N+1);
y=zeros(1,N+1);
y(1)=y0;
yX=yExact(t);
for k=1:N
y(k+1)=eulerImprovedStep(@f,t(k),y(k),h);
end
plot(t,y,t,yExact(t),'--')
xlabel('$t$');ylabel('$y$')
legend('Euler','Exact')
error=abs(yX(end)-y(end));
fprintf('N = %i, error=%0.3f\n',N,error)
N = 36, error=0.010

## Runge Kutta can get away with points

Clearly this is not enough points to generate a nice smooth curve but the numbers it gives are very accurate.
N=18;
h=(tEnd-tStart)/N;
t=linspace(tStart,tEnd,N+1);
y=zeros(1,N+1);
y(1)=y0;
yX=yExact(t);
for k=1:N
y(k+1)=RKStep(@f,t(k),y(k),h);
end
plot(t,y,t,yExact(t),'--')
xlabel('$t$');ylabel('$y$')
legend('Euler','Exact')
error=abs(yX(end)-y(end));
fprintf('N = %i, error=%0.3f\n',N,error)
N = 18, error=0.006

## Show the Euler method is first order

N=1000;
h=(tEnd-tStart)/N;
t=linspace(tStart,tEnd,N+1);
y=zeros(1,N+1);
y(1)=y0;
for k=1:N
y(k+1)=eulerStep(@f,t(k),y(k),h);
end
error1 = abs(yExact(tEnd)-y(end));
fprintf('N = %i, error=%0.3f\n',N,error1)
N = 1000, error=0.019
N=2*N;
h=(tEnd-tStart)/N;
t=linspace(tStart,tEnd,N+1);
y=zeros(1,N+1);
y(1)=y0;
yX=yExact(t);
for k=1:N
y(k+1)=eulerStep(@f,t(k),y(k),h);
end
error2=abs(yExact(tEnd)-y(end));
fprintf('N = %i, error=%0.3f\n',N,error2)
N = 2000, error=0.010
errRatio=error1/error2;
fprintf('The error went down by a factor of %0.2f.\n',errRatio)
The error went down by a factor of 1.99.

## Show the Improved Euler method is second order

N=1000;
h=(tEnd-tStart)/N;
t=linspace(tStart,tEnd,N+1);
y=zeros(1,N+1);
y(1)=y0;
for k=1:N
y(k+1)=eulerImprovedStep(@f,t(k),y(k),h);
end
error1 = abs(yExact(tEnd)-y(end));
fprintf('N = %i, error=%0.d\n',N,error1)
N = 1000, error=9e-06
N=2*N;
h=(tEnd-tStart)/N;
t=linspace(tStart,tEnd,N+1);
y=zeros(1,N+1);
y(1)=y0;
yX=yExact(t);
for k=1:N
y(k+1)=eulerImprovedStep(@f,t(k),y(k),h);
end
error2=abs(yExact(tEnd)-y(end));
fprintf('N = %i, error=%0.3d\n',N,error2)
N = 2000, error=2.208e-06
errRatio=error1/error2;
fprintf('The error went down by a factor of %0.2f.\n',errRatio)
The error went down by a factor of 3.87.

## Show the Runge-Kutta method is fourth order

N=1000;
h=(tEnd-tStart)/N;
t=linspace(tStart,tEnd,N+1);
y=zeros(1,N+1);
y(1)=y0;
for k=1:N
y(k+1)=RKStep(@f,t(k),y(k),h);
end
error1 = abs(yExact(tEnd)-y(end));
fprintf('N = %i, error=%0.d\n',N,error1)
N = 1000, error=1e-10
N=2*N;
h=(tEnd-tStart)/N;
t=linspace(tStart,tEnd,N+1);
y=zeros(1,N+1);
y(1)=y0;
yX=yExact(t);
for k=1:N
y(k+1)=RKStep(@f,t(k),y(k),h);
end
error2=abs(yExact(tEnd)-y(end));
fprintf('N = %i, error=%0.3d\n',N,error2)
N = 2000, error=6.707e-12
errRatio=error1/error2;
fprintf('The error went down by a factor of %0.2f.\n',errRatio)
The error went down by a factor of 15.58.

## Functions used

function yprime= f(t,y)
yprime=(sin(t)*(y^3-y));
end
function y=yExact(t)
y=exp(cos(t))./sqrt(3*exp(2)+exp(2*cos(t)));
end
function yNew=eulerStep(f,t,y,h)
yNew=y+h*f(t,y);
end
function yNew=eulerImprovedStep(f,t,y,h)
k1=f(t,y);
k2=f(t+h,y+h*k1);
yNew=y+h*(k1+k2)/2;
end
function yNew=RKStep(f,t,y,h)
k1=f(t,y);
k2=f(t+h/2,y+h*k1/2);
k3=f(t+h/2,y+h*k2/2);
k4=f(t+h,y+h*k3);
yNew=y+h*(k1+2*k2+2*k3+k4)/6;
end