6 views (last 30 days)

Show older comments

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;

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;

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?

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)));

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?

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

Start Hunting!