orbit equation with ode45

19 vues (au cours des 30 derniers jours)
Douglas Alves
Douglas Alves le 30 Oct 2014
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

Réponses (1)

Torsten
Torsten le 30 Oct 2014
y(2) contains the solution of the ODE
-u'' = -u + K.mi*K.gamma/K.l^2
I guess you either want to solve
u'' = -u + mi*gamma/l
as written as a comment behind your orbit-function, or
u'' = u + mi*gamma/l
as written in your introductionary text.
Best wishes
Torsten.

Catégories

En savoir plus sur Reference Applications dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by