orbit equation with ode45

24 views (last 30 days)
Douglas Alves
Douglas Alves on 30 Oct 2014
Answered: Torsten on 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

Answers (1)

Torsten
Torsten on 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.

Categories

Find more on Reference Applications in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!