Advanced Search

Question: Trying to make a simulated orbit accurate for thousands of orbits

Question:Trying to make a simulated orbit accurate for thousands of orbits

Aragorn 20 Maple
ode dsolve numeric differential-equations duplicate-question
December 10 2020
1

I have to create a simulation of the orbit of Hyperion around Saturn, which should be reliable for up to thousands of orbits. In practice I don't seem to be able to make the model stay consistent above a few hundred orbits using dsolve for the ODEs, but is there any way to do this so that the orbit stays on the same path?

I have uploaded my code here for your perusal.


 

``

> 

with(plots)
> 

G := 0.667408e-10
> 

MH := 0.55855e19
> 

MSat := 0.56832e27
> 

M := MSat+MH
> 

a := 0.1501e10
> 

ecc := .232
> 

beta := .89
> 

Digits := 10
> 

TH := sqrt(4*3.141592653589793^2*a^3/(G*M))
> 

omegaH := (2*3.141592653589793)/TH
> 

Eqns := diff(xH(t), t) = vxH(t), diff(yH(t), t) = vyH(t), diff(vxH(t), t) = -G*M*xH(t)/(xH(t)^2+yH(t)^2)^(3/2), diff(vyH(t), t) = -G*M*yH(t)/(xH(t)^2+yH(t)^2)^(3/2), diff(theta(t), t) = omega(t), diff(omega(t), t) = -G*MSat*beta^2*(xH(t)*sin(theta(t))-yH(t)*cos(theta(t)))*(xH(t)*cos(theta(t))+yH(t)*sin(theta(t)))/(xH(t)^2+yH(t)^2)^2.5
> 

ICs := xH(0) = a*(1-ecc), yH(0) = 0., vxH(0) = 0., vyH(0) = sqrt(G*M*(1+ecc)/(a*(1-ecc))), theta(0) = 0., omega(0) = omegaH
> 

soln := dsolve({Eqns, ICs}, numeric, method = rkf45, abserr = 0.1e-9, relerr = 0.1e-9, maxfun = 0)
> 

T := 100*TH
> 

odeplot(soln, [xH(t)/a, yH(t)/a], 0 .. T, scaling = constrained, labels = ["x/a", "y/a"], numpoints = 2000)``
> 

odeplot(soln, [theta(t), omega(t)], 0 .. T, labels = ["theta", "omega"], numpoints = 2000)

``


 

Download Project.mw
Please Wait...