I’m trying to simulate a one-body problem: Two-Body Numerical Solution in an Inertial Frame — Orbital Mechanics & Astrodynamics
\ddot{r} = \frac{-GM}{|r|^3}r = \frac{\mu}{|r|^3}r
Which can be split into:
\ddot{r_x} =\frac{\mu}{|r|^3}r_x
\ddot{r_y} =\frac{\mu}{|r|^3}r_y
If I manually create a state space model of it, it works as expected:
# μ is std gravitational parameter
μ = 1.0
# r is our orbital radius
r = 100
#From orbital velocity v = sqrt(μ/r)
# ṙy = sqrt(μ/r) = sqrt(1/100) = 0.1
# ṙ = A*r
# R is state vector [x,y]
# ṙx rx ṙy ry
R = [0,100,0.1,0.0]
A = [0 -μ/r^3 0 0
1 0 0 0
0 0 0 -μ/r^3
0 0 1 0
]
function fn!(dR,R,p,t)
dR .= A*R
end
Output is:
julia> sol = solve(prob)
retcode: Success
Interpolation: 3rd order Hermite
t: 23-element Vector{Float64}:
0.0
0.014072221109095914
0.15479443220005507
1.5620165431096464
15.634237652205558
96.97240029695304
271.1932065940258
522.6516504395706
843.4302866085757
1245.9148660620922
1708.1541865456325
2178.9375254988713
2747.703974500907
3370.4647853707534
3975.243341613904
4677.755465434197
5390.680246830623
6112.981695144672
6899.484518437533
7641.912879733987
8470.982854946626
9234.576078883347
10000.0
u: 23-element Vector{Vector{Float64}}:
[0.0, 100.0, 0.1, 0.0]
[-1.4072221108631458e-6, 99.99999999009863, 0.09999999999009863, 0.0014072221108631462]
[-1.547944315818753e-5, 99.9999988019342, 0.09999999880193419, 0.015479443158187531]
[-0.00015620159079168211, 99.99987800524076, 0.09999987800524075, 0.15620159079168203]
[-0.0015633600748802023, 99.98777877959019, 0.09998777877959017, 1.5633600748802023]
[-0.00968204893840455, 99.53018601465037, 0.09953018601465036, 9.68204893840455]
[-0.026788122517086115, 96.34519433183367, 0.09634519433183367, 26.788122517086112]
[-0.04991795510010098, 86.64985620491788, 0.08664985620491787, 49.91795510010099]
[-0.07469283623861765, 66.49044686598229, 0.06649044686598231, 74.69283623861764]
[-0.0947688643780251, 31.919595942781413, 0.03191959594278143, 94.76886437802513]
[-0.0990580914338696, -13.692780937011964, -0.013692780937011987, 99.05809143386973]
[-0.08207098337516477, -57.13450045574537, -0.057134500455745474, 82.07098337516476]
[-0.038377611821215296, -92.34266651347288, -0.09234266651347289, 38.37761182121527]
[0.022689372568765884, -97.39217126710521, -0.09739217126710525, -22.689372568765894]
[0.07404065795574148, -67.21642583667328, -0.06721642583667317, -74.0406579557415]
[0.09994106260058923, -3.4589246742302793, -0.003458924674229971, -99.94106260058923]
[0.07786227017821544, 62.751711755656665, 0.06275171175565686, -77.86227017821577]
[0.0169311147981259, 98.55858780280901, 0.09855858780280946, -16.931114798125872]
[-0.057813359031274894, 81.59899573665678, 0.0815989957366571, 57.8133590312753]
[-0.09776738729401874, 21.035901878435276, 0.021035901878434956, 97.76738729401917]
[-0.0815569509416117, -57.87876095288729, -0.05787876095288781, 81.55695094161223]
[-0.01888607969783396, -98.20921579651987, -0.09820921579652088, 18.88607969783399]
[0.0544281983351936, -83.90207159436801, -0.08390207159436888, -54.428198335194296]
The plot shows a cosine and sine wave, with expected period:
When I tried to use SecondOrderODEProblem:
function fn2!(ddR,dR,R,p,t)
ddR[1] = R[1] -μ/r^3 # r̈x = rx - μ/(|r|^3)
ddR[2] = R[2] -μ/r^3 # r̈y = ry - μ/(|r|^3)
end
# ṙx ṙy rx ry
prob2 = SecondOrderODEProblem(fn2!, [0,0.1], [100,0.0], tspan)
It explodes:
julia> sol2 = solve(prob2)
┌ Warning: At t=703.5196841226399, dt was forced below floating point epsilon 1.1368683772161603e-13, and step error estimate = NaN. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase ~/.julia/packages/SciMLBase/dHarz/src/integrator_interface.jl:600
retcode: Unstable
Interpolation: 3rd order Hermite
t: 628-element Vector{Float64}:
0.0
1.407222124981799e-5
0.00015479443374799787
0.0015620165587297966
0.015634237808547786
0.0969731329706833
0.27126821080651753
0.5236150230721681
0.8470009810336272
1.2552488401436728
⋮
703.5196841226399
u: 628-element Vector{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}:
([0.0, 0.1], [100.0, 0.0])
([0.0014072221109560224, 0.09999999999582916], [100.00000000990137, 1.4072221249292297e-6])
([0.015479443281823324, 0.10000000104327142], [100.00000119806583, 1.5479443424637106e-5])
([0.15620171783026238, 0.10000012043279412], [100.00012199481007, 0.00015620171817233144])
([1.5634874571189656, 0.1000122060836606], [100.01222171841309, 0.0015634873505366552])
([9.712518927880645, 0.10047046088097153], [100.47055800145617, 0.009712514319425774])
([27.460741331814937, 0.10370166538149508], [103.7019399518917, 0.02746070458702246])
([54.78719834109313, 0.11402417522989955], [114.0247229616412, 0.0547870586417341])
([95.19709414632987, 0.1380660198069198], [138.06697139720106, 0.09519671442858303])
([161.18528000798202, 0.18968419120315827], [189.6858021591164, 0.16118438476180436])
([276.45096347045137, 0.29397881575861484], [293.9815783284611, 0.2764490264191585])
([495.2571730914945, 0.5052472348066727], [505.2521833259307, 0.4952531255221926])
([937.7026295469472, 0.943010504890484], [943.0198734866742, 0.9376942087251551])
)
⋮
([1.5567977511949736e307, 1.556782198785263e304], [1.5567977511949736e307, 1.556782198785263e304])
Is there anything missing?
