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)
采纳的回答
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)
15 个评论
David
2022-11-15
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
2022-11-15
Torsten
2022-11-15
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
2022-11-15
David
2022-11-15
Davide Masiello
2022-11-15
I have edited my answer, please check.
David
2022-11-15
Davide Masiello
2022-11-15
My pleasure. If the answer helped, please consider accepting it.
Torsten
2022-11-23
Look at the regions plotted in pink under
David
2022-11-23
Torsten
2022-11-23
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
2022-11-23
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
