How to solve ode45 error?

Hi!

I have to do a Restricted three body analisys, which launch a rocket to Moon,

The program is:

%constants G = 6.67384e-11; Mt = 5.9722e24; Ml = 7.349e22; m= 100; d = 384400000; w =2.59e-6; R = 6371000; fi0 = 41.61; v0 = 11000;

theta0 = 50;

%condiciones iniciales x0=zeros(1,4); x0(1) = R; %posicion inicial x x0(2) = fi0; %velocidad inicial x x0(3) = m*v0*cos(theta0-fi0); %posicion inicial y x0(4) = m*R*v0*sin(theta0-fi0); %velocidad inicial y

tf=3000; tspan=[0 tf];

fg=@(t,x)[ x(3)/m; x(4)/m*x(1)*x(1); (x(4)*x(4)/m*x(1)*x(1)*x(1))-(G*Mt*m/x(1)*x(1))-(G*Ml*m*(x(1)-d*cos(x(2)-w*t)))/(sqrt(x(1)*x(1)+d*d-2*x(1)*d*cos(x(2)-w*t)))^3 ; -G*Ml*m*x(1)*d*sin(x(2)-w*t)/(sqrt(x(1)*x(1)+d*d-2*x(1)*d*cos(x(2)-w*t)))^3 ]; [t,x]=ode45(fg,tspan,x0);

figure(1);

plot(x(:,1),x(:,2),'r')

I have the problem when run the send me a warning:

Warning: Failure at t=4.830103e-18. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.232595e-32) at time t. > In ode45 at 309

Can anybody help me?

Thanks.

Ben

Answer 1

[Star Strider]

I didn’t run your code, but that Warning usually appears when you have a discontinuity, such as a singularity, in your ODE function at that time.

The first step in finding the problem is to plot the output of ode45 to see what your function is doing. Look at the integrated function and the derivatives.

Consider starting your ‘tspan’ interval at some time other than zero (experiment with 0.001 for example) and see if that helps. Also look at all your initial conditions to be certain they are not creating a singularity or discontinuity near t=0.

[Star Strider]

I see several problems.

First, you need to vectorise your code if you intend element-by-element operations. This means replacing ‘*’ by ‘.*’, ‘/’ by ‘./’ and ‘^’ by ‘.^’. See the documentation for Array vs. Matrix Operations for details on vectorisation.

Second, you need to attend to parentheses (). For example, in the second row, you have:

x(4)/m*x(1)*x(1);

that is the same as:

x(4)*x(1)*x(1)/m;

if you want all the terms to the right of the division operator / to be in the denominator, it should be:

x(4)./(m.*x(1).*x(1));

(I added the vectorisation operators).

Answer 2

[Ali]

I think that could hapen, and i tried to change everething before asking in the post, but i have always the same problem....

I don't know, but the program is running a problem with polar arguments, can be the problem related with that? Or maybe the ODE is depening with time vector and may this can cause problems?