Versions of Euler Methods

Please, pardon me for my ignorance I am new to MATLAB. I want to compare three versions of Euler methods with the step sizes of n = [ 2^4, 2^5, 2^6, 2^7, 2^8, 2^9] so that I can compute the error to the reference solution(yref) where the error is the max( abs (yref - yi) ) so that I can print out a table for each of the three columns each ie, time step size, corresponding error and the error rate. Unfortunately, the code that I have wrriten works for only the step size of n = 2^4 (I am not even certain if I am correct). So I provide my codes below for anyone to possible help me out with the correct codes. Thanks in advance.

Here's my forward Euler method

function y = forwardEuler(func,t,y)
% solve the ODE y'=f(t,y)
% input: func is the name of the f function, 
% t is a vector of the t points: [t1,t2,...,tn]
% y1 is the initial condition
% output: the vector y 

% initialize y to be the same size at t, and let N be the total number of
% y points we want to find

T=0.5;
n_steps = 2^4;
t=linspace(0,T,1+n_steps);

y = 0 * t;
N = length(y);

% Set initial condition
y(1)=100;

% use FE to find y_i+1 
for i=1:N-1
    y(i+1) = y(i) + ( t(i+1) - t(i) ) * func(t(i),y(i));
end

Here's method A code

function y = methodA(func,t,y)
T=0.5;
n_steps = 2^4;
t=linspace(0,T,1+n_steps);

y = 0*t;
N = length(y);
y(1)=100; % Set initial condition

for i=1:N-1 % use FE to find y_i+1   
h=t(i+1)-t(i);  %step size
 %y(i+1) = y(i)+ h*func( t(i) + h/2, (y(i)+y(i+1))/2 ); earlier codes
 
 y(i+1) = y(i)+ h*func( t(i)      , y(i)            );
y(i+1) = y(i)+ h*func( t(i) + h/2, (y(i)+y(i+1))/2 ); %James code
end

method B code is:

function y = methodB(func,t,y)
T=0.5;
n_steps = 2^4;
t=linspace(0,T,1+n_steps);

y = 0*t;
N = length(y);
y(1)=100; % Set initial condition

for i=1:N-1 % use FE to find y_i+1
    
h=t(i+1)-t(i); %step size

    %y(i+1) = y(i)+ h*func(t(i) + h/2, (h/2)*func(t(i),y(i))); eaelier codes
y(i+1) = y(i)+ h*func(t(i) + h/2, y(i) + (h/2)*func(t(i),y(i))); %James code
end

callingAllcodes

func = @(t,y) 10*cos(t)-2*y;
y0=100;
T=0.5;
[t,y] = ode45(func,[0 T],y0);
plot(t,y,'-k')
ylim([0,100])
hold on;
n_steps=2^4; %n_steps = [2^4,2^5,2^6,2^7,2^8,2^9]; I would like to do all these steps at
%once
t=linspace(0,T,1+n_steps);
y = forwardEuler(func, t, y0);
max_FE=max(abs(y))
plot(t,y,'r-o')
 
y = methodA(func, t, y0);
max_BE=max(abs(y))
plot(t,y,'b-x')
y = methodB(func, t, y0);
max_FE=max(abs(y))
plot(t,y,'g-o')
% 
legend('ode45', 'forward Euler', 'methodA', 'methodB', 'location','se')
hold off;

Answer 1

[Jan]

You provide the time span and y0 as inputs already:

y = forwardEuler(func, t, y0)

Inside the functions, you redefine these inputs:

function y = forwardEuler(func,t,y)
T=0.5;                      % == t(1)
n_steps = 2^4;              % == length(t) - 1
t=linspace(0,T,1+n_steps);  % == t
y = 0 * t;                  % Overwrites y0 called "y" in inputs
N = length(y);
y(1)=100;                 
...    

Prefer to use the inputs, because this allows to modify them in the caller:

function y = forwardEuler(func, t, y0)
N = length(t);
y = zeros(1, N);
y(1) = y0;
for i = 1:N-1
    y(i+1) = y(i) + (t(i+1) - t(i)) * func(t(i), y(i));
end
end

Now you can call this method using a loop:

axes('NextPlot', 'add');  % As: hold on
for N = 2 .^ (2:9)
   t = linspace(0, T, N + 1);
   y = forwardEuler(func, t, y0);
   plot(t, y);
end

Moving the definition of the initial position to the caller helps to avoid problems like this:

% forwardEuler:
y(1)=100;

% methodB:
y(1)=1000;

[Hmm!]

Jan, I am getting errors with your codes. Could you re-run your code to check it? Would be glad if you could post a workable codes for me to run.

[Jan]

figure
axes('NextPlot', 'add');  % As: hold on
y0 = 100;
T  = 0.5;
func = @(t,y) 10 * cos(t) - 2 * y;
for N = 2 .^ (2:9)
   t = linspace(0, T, N + 1);
   y = forwardEuler(func, t, y0);
   plot(t, y);
end

function y = forwardEuler(func, t, y0)
N = length(t);
y = zeros(1, N);
y(1) = y0;
for i = 1:N-1
    y(i+1) = y(i) + (t(i+1) - t(i)) * func(t(i), y(i));
end
end

For me it is working. If you get error messages, please post a copy of them. It is easier to fix a problem, if it is known, what the problem is.

[Hmm!]

Jan and James, I am getting more confused with the codes. For you to get clearer understanding of my codes, here's what I intend to do attached in this picture.

[Jan]

@Hmm!: The question asks for 0 <= t <= 1. So why do you implement T=0.5 and t = linspace(0, T, N + 1) ? Why 0.5 and not 1.0?

Except for this detail, I have posted working code for the forwardEuler method. One of the problem of your code is, that you do not use the input arguments.

James has explained, where the errors of your implementations in the other two methods are. XYou can re-use my code and insert just one line for the other methods.

So what exactly does still confuse you?

Answer 2

[James Tursa]

methodA has a fundamental flaw:

y(i+1) = y(i)+ h*func( t(i) + h/2, (y(i)+y(i+1))/2 );

You can't use y(i+1) on the right hand side because it isn't known yet. What your code currently does is use 0 for y(i+1) on the rhs because that is what it is initialized to, but this is obviously incorrect. You need to fix this method. I am guessing that this method was supposed to take an Euler step first to get a preliminary y(i+1) value, and then use that in the averaging formula. E.g., something like this:

y(i+1) = y(i)+ h*func( t(i)      , y(i)            );
y(i+1) = y(i)+ h*func( t(i) + h/2, (y(i)+y(i+1))/2 );

methodB also has a fundamental flaw:

y(i+1) = y(i)+ h*func(t(i) + h/2, (h/2)*func(t(i),y(i)));

The 2nd input argument to func( ) is supposed to be a y value, but you are feeding it a delta y value based on the derivative. That 2nd argument should look like this instead, where the current y(i) is added to the delta y:

y(i+1) = y(i)+ h*func(t(i) + h/2, y(i) + (h/2)*func(t(i),y(i)));

[Hmm!]

James Tursa, this is how the two methods look like for your reference.

[Hmm!]

James, when I changed methodA and methodB to yours, I get errors with it. eg. for methoB (and similar for methodA) I get something like this

>> methodB

Not enough input arguments.

Error in methodB (line 29)

y(i+1) = y(i)+ h*func(t(i) + h/2, y(i) + (h/2)*func(t(i),y(i)));

[James Tursa]

Can you post your complete current code?

[Jan]

methodB still starts at 1000, while the other 2 methods start at 100.

Calling the function with inputs is useless, if you overwrite the inputs inside the functions:

forwardEuler(func, t, y0)
%                  ^  ^^  Then use these values

Did you read my answer?