I tried to solve this system of ODE and it took so long to give me results.

2013 3 13
1 个回答
回答已采纳
1 次查看（30 天）

0 个投票

function [xprime] = CRungeKutta1(t,x)
r=sqrt(x(1)^2+x(2)^2);
while (2<r & r<=2.01)
A=(16.3096-7.1548*r)/r;
L=-log(r-2);
B=2*(0.4056*L^2+1.49681*L-1.9108)/(r*(L^2+6.0425*L+6.32549));
end
while (2.01<r & r<2.5)
A=-4.3833+17.7176*r^-1+14.8204*r^-2-92.4471*r^-3-46.3151*r^-4+232.2304*r^-5;
B=-3.1918+12.3641*r^-1+11.4615*r^-2-65.2926*r^-3-36.4909*r^-4+154.8074*r^-5;
end
if (r>=2.5)
A=5*r^-3-8*r^-5+25*r^-6-35*r^-8+125*r^-9-102*r^-10+12.5*r^-11+430*r^-12;
B=(1/3)*(16*r^-5+10*r^-8-36*r^-10-25*r^-11-36*r^-12);
end
e=x(1).*x(2).*(B-A)/r.^2;
xprime=[x(2)+e.*x(1)-0.5*B.*x(2); e.*x(2)-0.5*B.*x(1)];
end
I type on command window
X0=[-10,0.2]
tspan=[0:0.1:200]
[t,x]=ode45(@CRungeKutta1,tspan,X0)
I tried to solve this system of ODE but it takes so long to do it. Anybody can tell me why? Please helps . Thanks

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

 采纳的回答

ChristianW
ChristianW 2013-3-13

0 个投票

Pretty sure your function hang up in the while loops. I guess it should be IF instead of WHILE.
while 2<r
This is an endless loop because you are not changing r in the loop.

5 个评论
显示 3更早的评论 隐藏 3更早的评论

Phong Pham
Phong Pham 2013-3-13
编辑：Phong Pham 2013-3-13
Can I use all If instead of while?
I replace all the while by If. Now they show the error message " Undefined B"
ChristianW
ChristianW 2013-3-13
编辑：ChristianW 2013-3-13
if (2<r && r<=2.01)
A=(16.3096-7.1548*r)/r;
L=-log(r-2);
B=2*(0.4056*L^2+1.49681*L-1.9108)/(r*(L^2+6.0425*L+6.32549));
elseif (2.01<r && r<2.5)
A=-4.3833+17.7176*r^-1+14.8204*r^-2-92.4471*r^-3-46.3151*r^-4+232.2304*r^-5;
B=-3.1918+12.3641*r^-1+11.4615*r^-2-65.2926*r^-3-36.4909*r^-4+154.8074*r^-5;
elseif (r>=2.5)
A=5*r^-3-8*r^-5+25*r^-6-35*r^-8+125*r^-9-102*r^-10+12.5*r^-11+430*r^-12;
B=(1/3)*(16*r^-5+10*r^-8-36*r^-10-25*r^-11-36*r^-12);
else % THIS IS FOR r<=2
A = ...
B = ...
end
You need to define A and B if the radius is becoming r<=2.
Phong Pham
Phong Pham 2013-3-13
编辑：Phong Pham 2013-3-13
function [xprime] = CRungeKutta1(t,x)
r=sqrt(x(1)^2+x(2)^2);
if r<=2
r=2
A=(16.3096-7.1548*r)/r;
L=-log(r-2);
B=2*(0.4056*L^2+1.49681*L-1.9108)/(r*(L^2+6.0425*L+6.32549));
elseif (2<=r && r<=2.01)
A=(16.3096-7.1548*r)/r;
L=-log(r-2);
B=2*(0.4056*L^2+1.49681*L-1.9108)/(r*(L^2+6.0425*L+6.32549));
elseif (2.01<r && r<2.5)
A=-4.3833+17.7176*r^-1+14.8204*r^-2-92.4471*r^-3-46.3151*r^-4+232.2304*r^-5;
B=-3.1918+12.3641*r^-1+11.4615*r^-2-65.2926*r^-3-36.4909*r^-4+154.8074*r^-5;
elseif (r>=2.5)
A=5*r^-3-8*r^-5+25*r^-6-35*r^-8+125*r^-9-102*r^-10+12.5*r^-11+430*r^-12;
B=(1/3)*(16*r^-5+10*r^-8-36*r^-10-25*r^-11-36*r^-12);
end
e=x(1).*x(2).*(B-A)/r.^2;
xprime=[x(2)+e.*x(1)-0.5*B.*x(2); e.*x(2)-0.5*B.*x(1)];
end
Undefined function or variable "B" in the "e" expression.
ChristianW
ChristianW 2013-3-13
编辑：ChristianW 2013-3-13
If your radius is r<=2 then you put r=2. Then L become inf because of log(0). Then B becomes NaN.
I don't know your system, but if the radius r can't be smaller then 2, then you may need a event function. For an example, check out:
ballode
open ballode
odeexamples

请先登录，再进行评论。

更多回答（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