orbit equation with ode45
24 views (last 30 days)
Show older comments
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
0 Comments
Answers (1)
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.
0 Comments
See Also
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!