Forward Euler method in MATLAB-CodePudding

I'm trying to write a MATLAB code for the Forward Euler method, but I don't think the output is completely correct. This is the question:

Write a Matalb function that implements the Forward Euler Scheme to solve generic systems of p ODEs. Your function must conform to the following interface function [tHist ,uHist] = ForwardEuler(f,tspan ,u0,nsteps) with the following inputs:

•a function handle f for the function f of (1);

•a 2-vector tspan, with components t0 and T;

•a p ×1 vector u0, containing the initial condition y0;

•the number of steps nteps to be taken;

and the following outputs:

•an (Nh 1) ×1 vector tHist, storing the approximation times tn

•an (Nh 1) ×p matrix uHist, whose nth row stores the approximation un ≈y(tn) ∈R^p.

This is my code:

function [tHist ,uHist] = ForwardEuler(f,tspan ,u0,nsteps)
%Inputs:
% f:  a function handle for the function f 
% tspan: a 2-vector tspan, with components t0 and T, the initial and final components respectively;
% u0 : a p ×1 vector containing the initial condition y0;
% nsteps : he number of steps to be taken;
% Outputs:
% tHist: a (Nh   1) ×1 vector  storing the approximation times tn
% uHist: an (Nh   1) ×p matrix uHist, whose nth row stores the approximation un ≈y(tn) ∈ R^p.
T = [tspan(1) :nsteps: tspan(2)];
u = zeros(1,nsteps);
u(n) = u0;
h = (tspan(2) - tspan(1))/nsteps;
for i= 1: nsteps 
    u(i 1) = u(i)   h*f(T(i),u(i));
    tHist = [u(i 1)*T(i)];
    uHist = [u(i 1)];
end
end

CodePudding user response：

A minimal straightforward implementation could look like

function [tHist ,uHist] = ForwardEuler(f,tspan ,u0,nsteps)
    tHist = linspace(tspan(1), tspan(2), nsteps 1);
    h = tHist(2)-tHist(1);
    uHist = zeros(nsteps 1,length(u0));
    uHist(1,:)=u0;
    for i = 1:nsteps
        uHist(i 1,:) = uHist(i,:)   h*f(tHist(i),uHist(i,:));
    end%for
end%function

Test the routine with the circle equation

[T,U] = ForwardEuler(@(t,u)[u(2), -u(1)], [0,6.5], [1,0], 60);
plot(U(:,1),U(:,2));

which produces the expected outward spiral.

CodePudding user response：

The core algorithm seems to be alright only the initial values had to be corrected.

function [T ,u] = ForwardEuler(f,tspan ,u0,nsteps)
u = zeros(1,nsteps 1);%corrected to nsteps 1
u(1) = u0;%corrected n to 1
h = (tspan(2) - tspan(1))/nsteps;
T = [tspan(1) :h: tspan(2)];%corrected nsteps to h
for i= 1: nsteps 
    u(i 1) = u(i)   h*f(T(i),u(i));
end
end

tHist and uHist still not done but there is not much work left to do.

Another solution is to store for each step the T and u vectors in matrices.

function [tHist ,uHist] = ForwardEulerWithHist(f,tspan ,u0,nsteps)
tHist=zeros(nsteps 1,nsteps);
uHist=zeros(nsteps 1,nsteps);
   for j=1:nsteps
    [T,u] = ForwardEuler(f,tspan,u0,j);
    tHist(1:(j 1),j) = T;
    uHist(1:(j 1),j) = u;
  end
end

Then you may compare how much the iterations smooths to the best iteration.Example compare iteration at step 20 with iteration at step 10:

[tHist,uHist]=ForwardEulerWithHist(@(x,y)x y.^2,[0,0.5],0,20);
plot(tHist(:,20),uHist(:,20)), hold on;plot(tHist(1:11,10),uHist(1:11,10))