2

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 my code:

function [t ,u] = Euler(f,tspan ,u0,n)

T = [tspan(1) :n: tspan(2)];
u = zeros(1,n);
u(n) = u0;
h = (tspan(2) - tspan(1))/n;
for i= 1: n 
    u(i+1) = u(i) + h*f(T(i),u(i));
    t = [u(i+1)*T(i)];
    u = [u(i+1)];
end
end
Nancy
  • 33
  • 8
  • 1
    Please note that `u(i)` will always be a scalar (element `i` of the data of `u` in column-first order). If you want to address the vector, you need matrix or 2D array notation, `u(i,:)`. // See https://stackoverflow.com/questions/27657685/what-is-the-fastest-way-of-appending on how to implement appending in Matlab. // Contemplate on the difference of `u` and `uHist`. The same for `T` and `tHist`. – Lutz Lehmann Jan 01 '22 at 08:56
  • Thank you so much! For uHist and tHist I really can't figure out what else I should do. From the. question all I could figure out was that the approximation should be u(i+1) * T(i). As for uHist, I don't get where the p should come from, and don't know what to do. Is it, for each ODE I give, the columns are the approximation at that row, each for every function 1:p ? – Nancy Jan 01 '22 at 17:29
  • 1
    p is the dimension of the ODE system, in its first-order formulation. As I understand it, your code should handle not only scalar first-order equations. The arrays u and T are already what uHist and tHist should contain. I do not see where you get the idea of the product `u(i+1) * T(i)` from. – Lutz Lehmann Jan 01 '22 at 18:43
  • "an (Nh + 1) ×1 vector tHist, storing the approximation times tn " this is where I got it from. I can't figure out how else to solve it – Nancy Jan 01 '22 at 20:59
  • 1
    Ok, one can indeed read that wrong. It actually means to store the times `t_n` at which the approximations are computed, the "approximation times", not a multiplication operation. – Lutz Lehmann Jan 01 '22 at 23:25
  • Ooh ok, thank you soo much!! I did read it wrong. So is tHist = T(i) in the for loop and uHist = u(i+1)? – Nancy Jan 02 '22 at 09:48
  • 1
    Not exactly. It is tHist=T and uHist=u, the history variables are to contain the lists of the times and approximations computed during the algorithm. It would make more sense to have a t and u variable that are just the current value and state vector, and use them after each loop to update the history arrays. – Lutz Lehmann Jan 02 '22 at 10:33
  • You need to use the variable names in the output in the function declaration. So either uppercase `T` everywhere or lowercase `t`. In the construct `a:h:b`, `h` is the step size. If you want to prescribe the step number use `linspace` or compute the correct value of `h=(b-a)/N` and then the array as `t=a+[0:N]*h` (to force the exact number of steps). – Lutz Lehmann Jan 06 '22 at 11:21

1 Answers1

1

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.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51