MATLAB Answers

Versions of Euler Methods

2 views (last 30 days)
Hmm!
Hmm! on 28 Apr 2021
Commented: Jan on 30 Apr 2021
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;

Accepted Answer

Jan
Jan on 29 Apr 2021
Edited: Jan on 29 Apr 2021
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;
  4 Comments
Jan
Jan on 30 Apr 2021
@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?

Sign in to comment.

More Answers (1)

James Tursa
James Tursa on 29 Apr 2021
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)));
  4 Comments
Jan
Jan on 30 Apr 2021
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?

Sign in to comment.

Products


Release

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!