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

 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

Answer 1

[ChristianW]

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.

[Phong Pham]

Can I use all If instead of while?

I replace all the while by If. Now they show the error message " Undefined B"

[ChristianW]

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]

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]

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

[Phong Pham]

I just got it. Thanks