时间戳

首页 > 解决方案 > 为什么 solve_ivp 不适用于我的二阶非线性 ODE 系统？

python - 为什么 solve_ivp 不适用于我的二阶非线性 ODE 系统？

问题描述

我的任务是重现早期研究论文中的一些已知结果。为此，我需要求解 ODE 的非线性二阶系统。但是，我的代码无法在我指定的时间间隔内求解方程（即列表“结果”为空），从而产生不正确的结果。我已经尝试了几种不同的解决方案，例如包括 t_eval、更改 max_step、尝试不同的方法、在一个方向（即在不同的 ws 上）扫描参数空间，但似乎没有一个有效。附上问题出处的代码

    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)

任何形式的帮助将不胜感激！:)

标签： pythonmathscipyode

解决方案

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

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

最后，我认为您的循环更新可能是错误的。你可能想要一个w = w0 + i * scanwidth和定义w0上面或w += scanwidth

推荐阅读