2

Been trying to solve the newtonian two-body problem using RK45 from scipy however keep running into the TypeError:'Required step size is less than spacing between numbers.' I've tried different values of t_eval than the one below but nothing seems to work.

from scipy import optimize
from numpy import linalg as LA
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import numpy as np
from scipy.integrate import solve_ivp

AU=1.5e11
a=AU
e=0.5
mss=2E30
ms = 2E30
me = 5.98E24
mv=4.867E24
yr=3.15e7
h=100
mu1=ms*me/(ms+me)
mu2=ms*me/(ms+me)
G=6.67E11
step=24

vi=np.sqrt(G*ms*(2/(a*(1-e))-1/a))
#sun=sphere(pos=vec(0,0,0),radius=0.1*AU,color=color.yellow)
#earth=sphere(pos=vec(1*AU,0,0),radius=0.1*AU)

sunpos=np.array([-903482.12391302, -6896293.6960525, 0.  ])
earthpos=np.array([a*(1-e),0,0])

earthv=np.array([0,vi,0])
sunv=np.array([0,0,0])





def accelerations2(t,pos):
    norme=sum( (pos[0:3]-pos[3:6])**2 )**0.5
    gravit = G*(pos[0:3]-pos[3:6])/norme**3
    sunaa = me*gravit
    earthaa = -ms*gravit
    tota=earthaa+sunaa
    return [*earthaa,*sunaa]

def ode45(f,t,y,h):
        """Calculate next step of an initial value problem (IVP) of an ODE with a RHS described
        by the RHS function with an order 4 approx. and an order 5 approx.
        Parameters:
        t: float. Current time.
        y: float. Current step (position).
        h: float. Step-length.
        Returns:
        q: float. Order 2 approx.
        w: float. Order 3 approx.
        """

        s1 = f(t, y[0],y[1])
        s2 = f(t + h/4.0, y[0] + h*s1[0]/4.0,y[1] + h*s1[1]/4.0)
        s3 = f(t + 3.0*h/8.0, y[0] + 3.0*h*s1[0]/32.0 + 9.0*h*s2[0]/32.0,y[1] + 3.0*h*s1[1]/32.0 + 9.0*h*s2[1]/32.0)
        s4 = f(t + 12.0*h/13.0, y[0] + 1932.0*h*s1[0]/2197.0 - 7200.0*h*s2[0]/2197.0 + 7296.0*h*s3[0]/2197.0,y[1] + 1932.0*h*s1[1]/2197.0 - 7200.0*h*s2[1]/2197.0 + 7296.0*h*s3[1]/2197.0)
        s5 = f(t + h, y[0] + 439.0*h*s1[0]/216.0 - 8.0*h*s2[0] + 3680.0*h*s3[0]/513.0 - 845.0*h*s4[0]/4104.0,y[1] + 439.0*h*s1[1]/216.0 - 8.0*h*s2[1] + 3680.0*h*s3[1]/513.0 - 845.0*h*s4[1]/4104.0)
        s6 = f(t + h/2.0, y[0] - 8.0*h*s1[0]/27.0 + 2*h*s2[0] - 3544.0*h*s3[0]/2565 + 1859.0*h*s4[0]/4104.0 - 11.0*h*s5[0]/40.0,y[1] - 8.0*h*s1[1]/27.0 + 2*h*s2[1] - 3544.0*h*s3[1]/2565 + 1859.0*h*s4[1]/4104.0 - 11.0*h*s5[1]/40.0)
        w1 = y[0] + h*(25.0*s1[0]/216.0 + 1408.0*s3[0]/2565.0 + 2197.0*s4[0]/4104.0 - s5[0]/5.0)
        w2 = y[1] + h*(25.0*s1[1]/216.0 + 1408.0*s3[1]/2565.0 + 2197.0*s4[1]/4104.0 - s5[1]/5.0)
        q1 = y[0] + h*(16.0*s1[0]/135.0 + 6656.0*s3[0]/12825.0 + 28561.0*s4[0]/56430.0 - 9.0*s5[0]/50.0 + 2.0*s6[0]/55.0)
        q2 = y[1] + h*(16.0*s1[1]/135.0 + 6656.0*s3[1]/12825.0 + 28561.0*s4[1]/56430.0 - 9.0*s5[1]/50.0 + 2.0*s6[1]/55.0)

        return w1,w2, q1,q2
t=0
T=10**5
poss=[-903482.12391302, -6896293.6960525, 0. ,a*(1-e),0,0 ]
sol = solve_ivp(accelerations2, [0, 10**5], poss,t_eval=np.linspace(0,10**5,1))
print(sol)

Not sure what the error even means because I've tried many different t_evl and nothing seems to work.

Linus
  • 147
  • 4
  • 15

1 Answers1

5

The default values in solve_ivp are made for a "normal" situation where the scales of the variables are not too different from the range from 0.1 to 100. You could achieve these scales by rescaling the problem so that all lengths and related constants are in AU and all times and related constants are in days.

Or you can try to set the absolute tolerance to something reasonable like 1e-4*AU.

It also helps to use the correct first order system, as I told you recently in another question on this topic. In a mechanical system you get usually a second order ODE x''=a(x). Then the first order system to pass to the ODE solver is [x', v'] = [v, a(x)], which could be implemented as 

def firstorder(t,state):
    pos, vel = state.reshape(2,-1);
    return [*vel, *accelerations2(t,pos)]

Next it is always helpful to apply the acceleration of Earth to Earth and of the sun to the sun. That is, fix an order of the objects. At the moment the initialization has the sun first, while in the acceleration computation you treat the state as Earth first. Switch all to sun first

def accelerations2(t,pos):
    pos=pos.reshape(-1,3)
    # pos[0] = sun, pos[1] = earth
    norme=sum( (pos[1]-pos[0])**2 )**0.5
    gravit = G*(pos[1]-pos[0])/norme**3
    sunacc = me*gravit
    earthacc = -ms*gravit
    totacc=earthacc+sunacc
    return [*sunacc,*earthacc]

And then it never goes amiss to use the correctly reproduced natural constants like

 G = 6.67E-11

Then the solver call and print formatting as

state0=[*sunpos, *earthpos, *sunvel, *earthvel]
sol = solve_ivp(firstorder, [0, T], state0, first_step=1e+5, atol=1e-6*a)
print(sol.message)
for t, pos in zip(sol.t, sol.y[[0,1,3,4]].T): 
    print("%.6e"%t, ", ".join("%8.4g"%x for x in pos))

gives the short table

The solver successfully reached the end of the integration interval.

       t         x_sun       y_sun    x_earth    y_earth
0.000000e+00 -9.035e+05, -6.896e+06,  7.5e+10,        0
1.000000e+05 -9.031e+05, -6.896e+06, 7.488e+10, 5.163e+09

that is, for this step the solver only needs one internal step.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • Thank you it really helped. I don't fully understand what comes out of the solver. I understand the first column you printed are simply the time step but I don't see what the other four are supposed to represent? I would like to plot one of the planet's position like (earthx,earthy). – Linus Jan 11 '20 at 20:31
  • The table has t as first column, then x,y of the sun and then x,y of the earth, the z coordinates left out as they are zero. Note that 1e5 sec are a little more than a day, that is, about 1° of the full orbit circle. – Lutz Lehmann Jan 11 '20 at 20:57
  • Thank you again. I was just wondering how to use the firstorder function on it's own. If I try to just run firstorder([0,T],state0) I get 'AttributeError: 'list' object has no attribute 'reshape''. – Linus Feb 14 '20 at 20:00
  • The function assumes that it is passed numpy arrays, as will happen coming from `solve_ivp`. Also, it can only process one time point at a time, else you would need to do some dimension testing to do the right thing with `reshape`. While the time argument is not used, whatever you intended to happen by passing it a tuple will not work. – Lutz Lehmann Feb 14 '20 at 22:16
  • what is "a" in atol=1e-6*a? – Martin Alexandersson Jun 21 '21 at 16:16
  • 1
    @MartinAlexandersson : `atol` is the absolute tolerance, it should be a value that makes sense relative to the typical scale of the state vector components. Here `a` is the major half-axis, or simply `1 AU`, as defined in the constant block in the question. – Lutz Lehmann Jun 21 '21 at 17:33
  • Is the method: **firstorder** an example of "reduction of second order linear equations to canonical forms"? – Martin Alexandersson Aug 20 '21 at 06:15
  • 1
    @MartinAlexandersson : I'm not sure about "canonical", there should be a mass matrix involved and impulse variables. So it only applies if the mass matrix is the identity matrix. And the gravitation equation is expressively non-linear, so it does not fit from that angle too. Here it is just the transformation of a second order equation to a partitioned first-order system (conceptually) and then to a standard first-order system – Lutz Lehmann Aug 20 '21 at 06:23
  • When I tested the original example above, changing to G = 6.67E-11 made the original code work! I guess that this scales the problem into something that is solvable. – Martin Alexandersson Aug 20 '21 at 06:42
  • 1
    @MartinAlexandersson : Yes, reducing the right side of some equation `y'=F(y)` by a factor of `1e-22` tends to slow down any divergence considerably. Or seen the other way, taking a reasonable ODE system and multiplying the right side with `1e+22` will push it out of operating range for most standard solvers. // Still, the resulting system is not the physical simulation of the solar system it was intended to be. – Lutz Lehmann Aug 20 '21 at 08:12