Problem with ode45 code ?

2 次查看（过去 30 天）
Jenny Briggs
Jenny Briggs 2016-2-12
评论： Brian Russell 2020-3-25
Hi, i'm really new to Matlab. I am using ode45 to solve this differential equation:
y1prime = y2
y2prime = -y1 - (1/t)*(y2)
I am plotting it over the closed interval 1, 50.
It has the intial conditions:
y1(0)=0
y2(0)=0
y2prime(0)=0.
However i keep getting error messages about my m-script. I am using a code from quite an old book, so i'm not sure if that has anything to do with it. I would appreciate if anybody could have a quick look to see where i've gone wrong. Thank you
function yprime = pend(t,y)
yprime = [y(2); -y(1)-((y(2))/t)];
tspan = [1 50];
yazero = [0;0]; ybzero = [0;0]; yczero = [0;0];
[ta, ya] = ode45(@pend, tspan, yazero);
[tb, yb] = ode45(@pend, tspan, ybzero);
[tc, yc] = ode45(@pend, tspan, yczero);
[y1,y2] = meshgrid(-5:.5:5,-3:.5:3);
Dy1Dt = y2; Dy2Dt = -y(1)-((y(2))/t);
quiver(y1,y2,Dy1Dt,Dy2Dt)
hold on
plot(ya(:,1),ya(:,2),yb(:,1),yb(:,2),yc(:,1),yc(:,2))
axis equal, axis([-5 5 -3 3])
xlabel y_1(t), ylabel y_2(t), hold off
end
  5 个评论
显示 3更早的评论
Brian Russell
Brian Russell 2020-3-25
I see your problem. The old book you refer to is "MATLB Guide" by Higham and Higham, which I think is in general a brilliant book. However, they often leave out statements in their code which they assume their reader will know to put in. I also struggled with this example until I realized that their function is very short. In fact, it is only two lines:
function yprime = pend(t,y)
yprime = [y(2); -y(1)-((y(2))/t)];
There should be an end after these two lines. This short function should be saved in a file called pend.m. The rest of your code should be put in another program, say pend_run.m. This will fix your problem.
Brian Russell
Brian Russell 2020-3-25
I meant "MATLAB Guide", of course!

请先登录，再进行评论。

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

回答（1 个）

John D'Errico
John D'Errico 2016-2-12
So, just suppose I decide to try to solve your problem, assuming that you really intended to solve it over the interval [1,50]. That would take no more complicated code than this:
yp = @(y,t) [y(2);-y(1) - y(2)./t);
[tout,yout] = ode45(yp,[1,50],[0 0]);
The result would be rather boring though, since you have postulated a differential equation system of the form:
y1' = 0
y2' = 0
Note that since y1(1) = y2(1) = 0 at the beginning, they will STAY at zero forever. And ever. And ever.
Kind of boring.

类别

在 Help Center 和 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