1

My task is to reproduce some known results from earlier research papers. To that end, I need to solve a nonlinear second-order system of ODEs. However, my code does not manage to solve the equations in my specified interval (i.e the list "result" is empty), thus yielding incorrect results. I have already tried several different solutions such as including t_eval, changing max_step, trying out different methods, scanning over the parameter space in one direction (i.e. over different w-s) but none seem to work. Attached is the code from which the problem originates

    import scipy.integrate as spint
    from math import exp
    from math import pi
    import numpy as np
    import matplotlib.pyplot as plt
    from math import sqrt

    #Model Parameters

    w = 0.7009561175558
    v = 0.05
    rmin = 1e-7
    rmaxx = 20

    #Initial Conditions
    initcond = [1e-24,1.0000000000000002,0.05450000000000014, 2.672546e-9]

    #For solve_ivp
    def eom(r, solss, w,v):
        m,a,phi,psi = solss
        V = phi**2 - 2*phi**4/v**2 + phi**6/v**4
        dVdr = 2*phi - 8*phi**3/v**2 + 6*phi**5/v**4
        rho = w**2*phi**2/a + (1-2*m/r)*psi**2 + V
        pres = w**2*phi**2/a + (1-2*m/r)*psi**2 - V
        dydr = [4*pi*r**2*rho, 2*r*a/(1 - 2*m/r)*(m/r**3 + 4*pi*pres), psi, -r*psi/(1 - 2*m/r)*(m/r**3 + 4*pi*pres) + (3*m*psi + r*(-2 + 4*pi*r**2*rho)*psi)/(r*(r - 2*m)) + (1/2*dVdr - w**2*phi/a)/(1-2*m/r)]
        return dydr

    result = []
    scansize = 1e-8
    scanpoints = 1000

    #Scan over parameter space to try out different w-s for which a smooth solution might exist.
    for i in range(scanpoints):
         try:
             w = w + i*scansize
             solutions = spint.solve_ivp(eom, t_span = [rmin,rmaxx], y0 = initcond, method = "Radau", args = (w, v), rtol = 1e-5)
             if solutions.status == 0:
                 result.append(solutions)
         except:
             pass
    # This would solve the ODE for a single w, for which I know that well-behaved solution exists.
    #solutions = spint.solve_ivp(eom, t_span = [rmin,rmaxx], y0 = initcond, method = "Radau", args = (w, v), max_step = 1e-7)

Any kind of help would be very much appreciated! :)

HiggsBoson
  • 11
  • 1

1 Answers1

1

Trying to run you code here first thing to note is that your try, except, pass is a very bad idea, as you cannot see error messages. It actually seems that the args are not properly passed to the function. This is discussed here. Following this post the solution for you would be

solutions = spint.solve_ivp( lambda t, y: eom(t, y, w, v), [rmin,rmaxx], initcond )

Finally, I think that your loop update might be wrong. You probably want a w = w0 + i * scanwidth and defining w0above or w += scanwidth

mikuszefski
  • 3,943
  • 1
  • 25
  • 38
  • Hey, thanks for your comment! I removed the `try, except` to explicitly see the errors, and also modified the `args` part as per your suggestion. Sadly, the error still persists, in that up to some timesteps the equation will be solved, but the "solution" diverges, causing the algorithm to terminate the iterations. Regardless, I do agree with your comment about the loop part, that was indeed wrong and has been fixed as suggested, thanks! – HiggsBoson Sep 25 '20 at 08:37
  • 1
    args is supported since v1.4 as described in your link – Matthew Flamm Oct 14 '20 at 23:37
  • @MatthewFlamm thanks for the info. I'll update....probably...sometime soon...I guess... – mikuszefski Oct 15 '20 at 10:22