Question

1. Implement the Explicit Euler Scheme for Initial Value Problems of the form:
        y'(t) = F(t, y(t)) ,  t0 ≤ t ≤ tend
        y(t0) = y0

   The function F(t,y) should be coded in a function subprogram FCN(...).
   Input data: t0, y0, tend, Nsteps.  Thus the time-step will be h=(tend-t0)/Nsteps.
   Your code should print out the input data and then the pairs:
                tn          Yn
   At the end, it should print out the final n, tn, Yn
   (appropriately labelled, of cource).

2. Solve, on paper, the simple (integration) problem:
                y' = 2t ,   0 ≤ t ≤ 1
                y(0) = −1

3. To debug your code, run it on the problem above. 
   Compare the numerical solution Yn with the exact solution yEXACT(tn),
   i.e. modify your output to print out:
        tn      Yn      yEXACTn         ERRn
   where ERRn = |Yn - yEXACT(tn)|, and keep track of the maximum error.
   At the end of the run, print out the above values (at time tend) 
   and the maximum overall error ERRmax.  Test with N=10 and N=100.
   Turn off printing of tn  Yn ... and test with N = 1000, 10000 and larger.

   Once the code is debugged, only FCN(...) and input data need to be changed 
   to solve other IVPs.

4. Now solve the IVP:
        y' = −t/y ,   0 ≤ t ≤ 1
        y(0) = 1

   Find the exact solution at t = 1 (by hand),
   and compare the numerical and exact values at t = 1.  
   Try small (N=10) and larger (N=100, 1000, 10000, ... ) number of time-steps.
   At which time does the worst error occur in this problem ?
   Plot the exact solution.  Do you see why it occurs there?

Answer 1

We have developed a MATLAB code for explicit Euler scheme for initial value problems.

MATLAB code:

function TT_Euler=Euler_explicit(t0,y0,tend,Nsteps)

fprintf('\n\nStarting value of t is, t0 = %g\n',t0);
fprintf('End value of t is, tn = %g\n',tend);
fprintf('Initial value of y is, y(%g) = %g\n',t0,y0);
fprintf('Number of steps for the given problem is, N = %g\n\n',Nsteps);
  
h= (tend-t0)/Nsteps;
tn=t0:h:tend;
Yn=y0;
A_Euler(1,:)=[tn(1) Yn Yn abs(Yn-Yn)];
  
for i=1:numel(tn)-1
[Fa,Fe]=Fty(tn(i),Yn);
yn1=Yn+h*Fa;
A_Euler(i+1,:)=[tn(i+1) yn1 Fe abs(Fe-yn1)];
Yn=yn1;
end
  
VarNames = {'tn','Yn','yEXACTn','ERRn'};
TT_Euler=table(A_Euler(:,1),A_Euler(:,2),A_Euler(:,3),A_Euler(:,4), 'VariableNames',VarNames);
  
if Nsteps<=100
fprintf('\n\nDetails of Euler method with h=%g\n\n',h)
disp(TT_Euler);
end
  
temp1=A_Euler(:,4)==max(A_Euler(:,4));
fprintf('\n\nValue of n=%g\nValue of tn=%g\nValue of Yn=%g\n',Nsteps,tn(end),Yn);
fprintf('Worst error occurs at t=%g and the worst error is = %g\n\n',A_Euler(temp1,1),max(A_Euler(:,4)));
  
plot(A_Euler(:,1),A_Euler(:,2),A_Euler(:,1),A_Euler(:,3),'--','LineWidth',1.5);
xlabel('t','fontsize',14)
ylabel('y(t)','fontsize',14)
legend('Numerical Solution','Exact Solution');
title('Numerical Solution vs Exact Solution','fontsize',14);
  
  
function [F_approx,F_exact]=Fty(t,y)
F_approx=-t./y; % 2*t; %
F_exact=sqrt(1-t.^2); % t.^2-1; %
end

end