Using DifferentialEquations to solve a gravity problem results in runaway behaviour

General Usage
1

Hi,

As a test and for learning, I’m using DifferentialEquations to solve a simple Newtonian gravity well motion problem. For the given initial conditions, I’m expecting circular motion, but I find that the orbit radius decays continuously. There is no drag force included in the system. I’m wondering whether I need to set some special solver fields, but not sure about it. I’d like to “learn to fish” here, so any tips on how to debug something like this would be greatly appreciated.

using ModelingToolkit, DifferentialEquations, Plots

function orbit(du,u,p,t)
    rx,ry,rz,vx,vy,vz = u
    G,M = p
    du[1] = vx
    du[2] = vy
    du[3] = vz
    r = sqrt(rx^2+ry^2+rz^2)
    du[4] = -G*M*rx/(r^3)
    du[5] = -G*M*ry/(r^3)
    du[6] = -G*M*rz/(r^3)
end

f = ODEFunction(orbit)
tspan = [0.0,240 * 3600.0]
R = 7300000.0
G = 6.674e-11
M = 5.972e24
V = sqrt(G*M/R) # So we get circular motion
prob = ODEProblem(f, [R,0.0,0.0,0.0,V,0.0], tspan, [G,M])

rad(u) = sqrt(u[1]^2+u[2]^2+u[3]^2)
ϕ(u) = acos(u[3]/rad(u))
θ(u) = atan(u[2],u[1])
sol = solve(prob)
plot(sol.t, rad.(sol.u), title="Distance from earth's centre")

Printing rad(sol.u[1]), rad(sol.u[end]) results in

(7.3e6, 1.3642779459358605e6)

Thanks for any pointers at all!
-Srikumar

2

For this simple problem, I can quite easily hand roll a solver that works as expected. So this is puzzling to me.

3

Your problem has nothing to do with the physics of the two body problem.

You have encountered an issue that arises because of the use of a numerical integrator. Typically, depending if your numerical integrator is explicit or implicit, they modify the real behaviour of your physical system in different ways (making trajectories converge to zero or diverge) depending on where the eigenvalues of your original problem lie in the so-called region of stability.

Therefore, if you are in such a case and you want the trajectory to not diverge much from your true trajectory, you can use a symplectic integrator (which conserves energy) or lower drastically your tolerances. For example, if you use

sol = solve(prob, Tsit5(), abstol=1e-14, reltol=1e-14)

you will see that the r is almost constant throughout the integration.

1 Like
4

Indeed! This works as expected … and thanks for both those links too.

1 Like
5

And it looks like solve picks Tsit5 as I get exactly the same results when not specifying it. So the tolerances were the key.

6

Yes, it is a fast and reliable explicit method for a general use case. To see the vast amount of numerical integrators Julia has, and a little bit of insight on which to use for a particular problem refer to the documentation.

1 Like
7

For reference, I wrote a bit about this:

3 Likes