为什么solve_ivp对我的二阶非线性ODE系统不起作用?

问题描述

我的任务是从早期的研究论文中复制一些已知的结果。为此,我需要解决ODE的非线性二阶系统。但是,我的代码无法在指定的时间间隔内解决方程式(即列表“结果”为空),从而产生错误的结果。我已经尝试了几种不同的解决方案,例如包括t_eval,更改max_step,尝试不同的方法,在一个方向上(即在不同的w-s上)扫描参数空间,但似乎没有任何效果。 附带的是导致问题的代码

    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,max_step = 1e-7) 

任何帮助将不胜感激! :)

解决方法

尝试在此处运行代码首先要注意的是 try,except,pass是一个非常糟糕的主意,因为您看不到错误消息。 实际上似乎args没有正确地传递给函数。 here对此进行了讨论。在这篇文章之后,您的解决方案将是

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

最后,我认为您的循环更新可能是错误的。您可能需要w = w0 + i * scanwidth并在w0上方定义w += scanwidth