Forward Euler solution plotting for dy/dt=y^2-y^3

Hi,

I am trying to solve the flame propagation model dy/dt=y^2-y^3 with y(0)= 1/100 and 0<t<200, using the forward and backward euler method with step size 0.01. But it has been giving me errors. How should I go about this? Please help I need this for my project

Thank you

Here are my codes for the Forward Euler

h=0.01;
y(0)=2
for n=1:N
t(n+1)=n*h
opts = odeset('RelTol',1.e-4);
y(n+1)= y(n)+h*(y.^2-y.^3);
end
plot(t,y)

Answer 1

[Davide Masiello]

There are a couple of issues with your code

1) Indexes in MatLab start at 1, not 0, so y(0) is not valid syntax and must be replaced with y(1).

2) First index, then raise to the power, i.e. y^2(n) becomes y(n)^2

See example below (since you have not specified the value of delta, I arbitrarily replaced it with 0.1)

N    = 100;                       % number of steps
t    = linspace(0,200,N);
h    = t(2)-t(1);               % step size

y    = zeros(1,N);              % Initialization of solution (speeds up code)
y(1) = 1/100;                   % Initial condition

for n = 1:N-1
    y(n+1) = y(n)+h*(y(n)^2-y(n)^3); % FWD Euler solved for y(n+1)
end

figure(1)
plot(t,y)

Now, the backward euler method is a bit more complicated because it's an implicit method.

You can use a root finder algorithm like fzero and do

N    = 100;                       % number of steps
t    = linspace(0,200,N);
h    = t(2)-t(1);               % step size

y    = zeros(1,N);              % Initialization of solution (speeds up code)
y(1) = 1/100;                   % Initial condition

for n = 1:N-1
    f       = @(x) x-y(n)-h*(x.^2-x.^3);
    y(n+1)  = fzero(f,y(n)); % BWD Euler solved for y(n+1)
end

figure(2)
plot(t,y)

[David]

Thanks for your help but the graph does not look right. The given differential y^2-y^3 is coming from the frame propagation model. Since I am working on solving stiff ode i am now trying to solve the given ode which is stiff with both Backward and forward euler to show the differences in stability isn the two solutions, Pls help if you have any clue on how to go about it. Thank you

[Torsten]

Change

h       = 0.01;         % step size
N       = 6;            % number of steps
y       = zeros(1,N);   % Initialization of solution (speeds up code)
y(1)    = 0.1;          % Initial condition

to

h       = 0.01;         % step size
N       = 100;          % number of steps
y       = zeros(1,N);   % Initialization of solution (speeds up code)
y(1)    = 2;            % Initial condition

in Davide's code.

syms y(t)
eqn = diff(y,t) ==y^2-y^3;
cond = y(0)==2;
sol = dsolve(eqn,cond);
y = matlabFunction(sol);
t = 0:0.01:1;
plot(t,y(t))

[David]

Sorry I made a mistake there by the initial condition the right pne is y(0)=1/100 and 0<t<200 Thank you

[Torsten]

And does it work with the new setting

h       = 0.01;         % step size
N       = 200*100;      % number of steps
y       = zeros(1,N);   % Initialization of solution (speeds up code)
y(1)    = 0.01;         % Initial condition

?

[David]

N=100

[Torsten]

If you set N = 100 with step size 1/100, you arrive at t=N*h = 1, not t=200.

[David]

Okay Lets just consider 200 then

[Davide Masiello]

I have edited my answer, please check.

[David]

Thanks very much

[Davide Masiello]

My pleasure. If the answer helped, please consider accepting it.

[David]

Good day sir, how do i graphically show the stability regions of the two methods(Backward and Forward Euler) on the same model?

[Torsten]

Look at the regions plotted in pink under

https://en.wikipedia.org/wiki/Euler_method

https://en.wikipedia.org/wiki/Backward_Euler_method

[David]

Thank you for the information. So they always look the same?

[Torsten]

They are not dependent on the model you solve - stability regions only depend on the discretization method for the ODE

y' = f(t,y)

[David]

Thank you very much, this helped me so much