Turning points in a squared first order ODE using Scipy.integrate-CodePudding

I am trying to code a radius function based on the Schwarzchild solution to a black hole given the expression:

(dr/dtau)^2= Emu^2- Veff^2

As it is a square the sign in front of the root will depend on the turning points that I have manually found and labeled tp1 and tp2. However, even though I am changing the functions sign depending on its position it behaves relatively well until it hits these turning points.

Here is the code I have so far: (P.S: I hope this is the correct formatting and way to present a question although i have been a reader for a few years this is actually my first post).

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode



tp1 = 1.08329*10**11
tp2 = 4.13115*10**11


#arbitrary initial radius
r_start = 3e11
#constants:
M = 4*10**6*(1.9891*10**30) # SMBH mass in kg
G = 6.67408*10**(-11) # Gravitational constant in N kg^-2 m^2
c = 299792458


Emu2 = 0.88*10**17
Lmu = 10**19
def odes(tau,rs):
    
    Vef = (1-(((2*G*M)/c**2)/rs))*((c**2) (Lmu/rs)**2)

    sign = (Emu2)-Vef
    signcount = 1
    if sign <= 0: 
        if rs <= tp1: 
            rs = tp1 5 
            Vef = (1-(((2*G*M)/c**2)/rs))*((c**2) (Lmu/rs)**2)
            sign = (Emu2)-Vef
            drdTau = np.sqrt(Emu2 - Vef)
            signcount = 1
        if rs >= tp2:
            rs = tp2-5
            Vef = (1-(((2*G*M)/c**2)/rs))*((c**2) (Lmu/rs)**2)
            sign = (Emu2)-Vef
            drdTau = (-1)*np.sqrt((Emu2) - Vef)
            signcount = 2
        return [drdTau]
    if sign > 0 :
        if signcount == 1: 
            Vef = (1-(((2*G*M)/c**2)/rs))*((c**2) (Lmu/rs)**2)
            sign = (Emu2)-Vef
            drdTau = np.sqrt(Emu2 - Vef)
        if signcount == 2: 
            Vef = (1-(((2*G*M)/c**2)/rs))*((c**2) (Lmu/rs)**2)
            sign = (Emu2)-Vef
            drdTau = (-1)*np.sqrt(Emu2 - Vef)
        return [drdTau]
            
if __name__ == '__main__': 
    
    r0 = [r_start]
    tspan = 5*3.154e7
    
    # timestep
    dtau = 1000
    
    # total number of steps
    n_steps = int(np.ceil(tspan/dtau))

    #Initialise arrays
    t = np.zeros((n_steps,1))
    rs = np.zeros((n_steps,1))
    ts = np.zeros((n_steps,1))
    step = 1


    r0 = [r_start]
    t0 = [0]
    t[0] = np.array(t0)
    rs[0] = np.array(r0)
# initiate solver
    solver = ode(odes)
    solver.set_integrator('DOP83')
    solver.set_initial_value(r0,0)
        
    #propagate orbit
    while solver.successful() and step<n_steps:
        solver.integrate(solver.t dtau)
        ts[step] = solver.t
        rs[step] = solver.y
        step  = 1

plt.plot(ts,rs,'s',color='#0066FF')
    

   
    # axes labels 
plt.xlabel('$x$') 
plt.ylabel('$y$')
plt.legend('pos')
    
    # check for and set axes limits
max_yval = np.amax(rs)
max_xval = np.amax(ts)
plt.xlim(0,max_xval)
plt.ylim(tp1 - 300000,max_yval)
    
plt.show()
print(rs)

CodePudding user response：

I have solved it by changing the function it is integrating once it reaches tp1 and tp2 I am hoping to add a root solver to automate the turning point solving but for now a rough tratment is what I have. I'll put up the code for anyone working on this area that has come across this issue once i have refined it, but essentially instead of modifying the function I added a second one called odes2with a negative sign and solved for that once the turning points were reached.

 solver = ode(odes)
 solver.set_integrator('lsoda')
 solver.set_initial_value(r0,0)
    
#propagate orbit
while solver.successful() and step<n_steps:
    solver.integrate(solver.t dtau)
    ts[step] = solver.t
    rs[step] = solver.y
    if rs[step]>tp2:
        rs[step] = tp2 
        solver = ode(odes2)
        solver.set_initial_value(rs[step],ts[step])
    if rs[step]<tp1:
        rs[step] = tp1 
        solver = ode(odes)
        solver.set_initial_value(rs[step],ts[step])   
    step  = 1

CodePudding user response：

Not perfect, but using the second order dynamic really avoids a lot of overhead

def Veff(rs): return (1-(((2*G*M)/c**2)/rs))*((c**2) (Lmu/rs)**2)

def grad_Veff(rs):
    rp=rs*(1 1e-5)-1e-5; rm = rs*(1 1e-5) 1e-5
    return (Veff(rp)-Veff(rm))/(rp-rm);

def ode(t,y): return [y[1], -grad_Veff(y[0])]

t = np.arange(0, tspan, dtau)

vstart = (Emu2 - Veff(rstart))**0.5

res = solve_ivp(ode, (0,tspan), [rstart, vstart], t_eval=t, atol=1e2, rtol=1e-8)

This should integrate through the turning points, to find the turning points search for sign changes or roots in the second component. One could also define an event for this so that the solver automatically registers the turning points.