yexact =
Infinite while loop for Euler's method in order to solve ode
显示 更早的评论
I am trying to solve an ode with Euler's method using while loop but i get into an infinite loop for my ymax value. If i change it to a lower value i don't have this problem. What is happening?
clc, clear all, close all
%--- Solving the diff. eqn. y'(t) = y(t) - c*y^2(t) with c = 0.5 using
%--- Euler's method in while loop
c = 0.5;
y0 = 0.1;
f = @(y) y - c*y.^2;
h = 0.1;
ymax = 1/c; % ymax is obtained by setting y' = 0 and solving for y, here is ymax = 2
t = 0;
y = 0.1;
[t,y]
while y < ymax % here is ymax = 2 and i get infinite loop. If, for example, i write y < 1 there is no problem, but i have to use ymax 2
y = y + h*f(y);
t = t + h;
[t,y]
end
采纳的回答
I think you misunderstand what is happening. Your ODE is
y' = y - y^2/2
ymax is found when y' = 0, and so you find that ymax = 2. All good so far. But, WHEN does that happen? I'll solve the ode here:
syms y(t)
c = 0.5;
y0 = 1;
ODE = diff(y) == y - c*y^2;
yexact = dsolve(ODE,y(0) == 0.1)
fplot(yexact,[0,10])
Now, when does y(t) exceed 2? Never, of course. But does it EVER reach 2 EXACTLY? NO. Ony at t==infinity. And infinity is a long way away. But your test in the while triggers only at y==2 (or greater, which can never happen.)
limit(yexact,t,inf)
while y < ymax % here is ymax = 2 and i get infinite loop. If, for example, i write y < 1 there is no problem, but i have to use ymax 2
And of course, if you change ymax to some lower value, it stops. Should that surprise you in the least bit? NO! Of course not.
7 个评论
OMG, you are so right and i feel so dumb right now. I changed it and, of course, it works. (I added some new lines trying to put the results in an array and plotting them).
clc, clear all, close all
%--- Solving the diff. eqn. y'(t) = y(t) - c*y^2(t) with c = 0.5 using
%--- Euler's method in while loop
c = 0.5;
f = @(y) y - c*y.^2;
h = 0.1;
ymax = 1/c;
t(1) = 0; y(1) = 0.1;
while y < 1.999999999999999
y(end+1) = y(end) + h*f(y(end));
t(end+1) = t(end) + h;
end
plot(t,y)
@John D'Errico Hello. When i run your code on my computer, the plot pops out but the axes limits are both from 0 to 1 instead being the ones we see below. What is wrong?
syms y(t)
c = 0.5;
y0 = 1;
ODE = diff(y) == y - c*y^2;
yexact = dsolve(ODE,y(0) == 0.1)
fplot(yexact,[0,10])
Walter Roberson
2024-12-30
When I run John's code, the plot limits are 0 to 10.
Left Terry
2024-12-30
Walter Roberson
2024-12-30
You are encountering a bug specific to R2016a. The FunctionLine object returned by the plotting is using the wrong limits. The work-around is to use
xlim([0 10]);
ylim([0 2.2]);
Left Terry
2024-12-31
Walter Roberson
2025-1-3
I tested the code in R2016a. When I did not use xlim(), the output was restricted to [0 1] when I fplot with upper bound exceeding 1 . When I used xlim(), the plot worked properly.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
