Infinite while loop for Euler's method in order to solve ode

7 次查看（过去 30 天）
Left Terry
Left Terry 2024-12-29
评论： Walter Roberson 2025-1-3
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

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2024-12-29
编辑：John D'Errico 2024-12-29
Ran in:
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)
yexact = 
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)
ans = 
2
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 个评论
显示 5更早的评论
Walter Roberson
Walter Roberson 2025-1-3
@Left Terry
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 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2016a

Community Treasure Hunt

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

Start Hunting!

Translated by