orbit equation with ode45

Question

Douglas Alves 2014-10-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/160650-orbit-equation-with-ode45

回答： Torsten 2014-10-30

Hi, I am trying to solve an equation with ode45 but it does not actually give me what it's expect to. The equation is the classical mechanics one for orbits. a(phi)'' = a(phi) + mi*gamma/l^2 , a(phi) = 1/r(phi)
I have included F = -gm1m2/r^2 already. I guess my initial values might be wrong but I do not see another alternative. ( I dont want to use the analitical solution for this). I put the code down below the first part is just a adjustment of scales plus variables. Could anyone help me see what's wrong? What I get is a vector y with y(1) = -y(2) which doesnt make sense. 
Since a = 1/r  and I get y1 and y2 ... and I have used y1=da/dphi and y2=a then, what I got is y1 and y2 and it's simple to convert them to give me radius and velocity but they dont return reasonable values. 
function SunsOrbit
% Simple 2 body problem with Jupiter and the Sun 
% CIRCULAR ORBIT AND JUPITER's VELOCITY CONSTANT WERE SUPPOSED
%
%Scales
Scale.d=7.785e11;
Scale.M=1.989e30;
Scale.mi=Scale.M;
Scale.V=1e3;
Scale.G=6.67384e-11;
Scale.gamma=Scale.G*Scale.M^2;
Scale.P=Scale.M*Scale.V;
Scale.l=Scale.d*Scale.P;
%Setting up Constants
K.d=7.785e11/Scale.d ; %m distance Jupiter-Sun
K.Mj=1.898e27/Scale.M; %kg
K.Msun= 1.989e30/Scale.M; %kg
K.mi=K.Mj*K.Msun/(K.Mj+K.Msun);
K.Vj=1.707e4/Scale.V; %m/s Jupiter's velocity
K.Pj=K.Mj*K.Vj; % Jupiter's linear momentum
K.l=K.d*K.Pj; % Angular Momentum Modulus
K.G=6.67384e-11/Scale.G; % m3 kg-1 s-2
K.gamma = K.G*K.Mj*K.Msun;
%
Inits=[1e-5,1/K.d];
ThetaRange=linspace(-pi,pi,500);
[r,y]=ode45(@orbit,ThetaRange,Inits,[],K); % y(1) is 1/r(theta)
r=1./y(:,2);
Vj=-(r.^2).*y(:,1);
Ueff=-K.gamma./r + (K.l.^2)./2*K.mi.*r.^2;
[u,i]=pol2cart(Theta,r);
plot(r*Scale.d,Vj*Scale.V,'*')
axis image;
% Orbit's equation
function da=orbit(Theta,y,K)
da=[-y(2)+K.mi*K.gamma/K.l^2;-y(1)]; % u'' = -u + mi*gamma/l 
                                    % 
                                    % TRY TO SET IT ALL UP SO THAT THE
                                    % ENERGY IS CONSTANT ALSO USE ENERGY AS
                                    % A FUNCTION OF POTENTIAL ENERGY AND
                                    % KINETIC ENERGY