计算科学 - 使用 odeint 避免不同的解决方案？拍摄方法 - 吾爱随笔录

我正在尝试在 Python 中求解一个方程。基本上我想做的是解方程：

\frac{1}{x^{2}} \frac{d}{d x} (G a m \frac{d L}{d x}) + L (\frac{a^{2} x^{2}}{G a m} - m^{2}) = 0

$\frac{1}{x^2}\frac{d}{dx}\left(Gam \frac{dL}{dx}\right)+L\left(\frac{a^2x^2}{Gam}-m^2\right)=0$

这是史瓦西时空中大规模标量场的克莱因-戈登方程。我们知道 $m$ 和 $Gam=x^2-2x$ .

我知道的初始/边界条件是 $L\rvert_{2+\epsilon}=1$ 和 $L\rvert_{\infty}=0$ .

请注意，方程的渐近行为是

L (x \to \infty) \to \frac{e^{\pm (m^{2} - a^{2}) x}}{x}

$L(x\to\infty)\to \frac{e^{\pm(m^2-a^2)x}}{x}$

那么，如果 $a^2>m^2$ 我们将有振荡解，而如果 $a^2 < m^2$ 我们将有一个发散和衰减的解决方案。我感兴趣的是衰变解决方案；然而，当我试图解决上述方程时，将其转换为一阶微分方程系统并使用射击方法来找到 $a$ 这可以给我我感兴趣的行为，我总是有一个不同的解决方案（对于 $0<a^2<m^2$ ）。我想它正在发生，因为odeint总是在寻找不同的渐近解。

有没有办法避免或告诉odeint我对衰变解决方案感兴趣？如果没有，你知道我可以解决这个问题的方法吗？也许使用另一种方法来解决我的微分方程组？如果是，是什么方法？

基本上我正在做的是添加一个新的方程组 $a$

\frac{d^{2} a}{d x^{2}} = 0, \frac{d a}{d x} (2 + ϵ) = 0, a (2 + ϵ) = a_{0})

$\frac{d^2a}{dx^2}=0,\\ \frac{da}{dx}(2+\epsilon)=0,\\ a(2+\epsilon)=a_0)$

为了有 $a$ 作为常数。然后我正在考虑不同的值 $a_0$ 并询问我的边界条件是否满足。

编辑 1。

你好，让我试着解释一下我的问题

考虑到渐近行为，我在无穷大处合并值，这意味着我将在场与其导数之间建立关系。如果有帮助，我会为您发布代码：

from IPython import get_ipython
get_ipython().magic('reset -sf')
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from math import *
from scipy.integrate import ode

这些是 Schwarzschild 的初始条件。该字段在重新缩放下是不变的，然后我可以使用 $L(2+\epsilon)=1$

def init_sch(u_sch):
    om = u_sch[0]
    return np.array([1,0,om,0]) #conditions near the horizon, [L_c,dL/dx,a,da/dx]

这些是我们的方程组

def F_sch(IC,r,rho_c,m,lam,l,j=0,mu=0):
    L = IC[0]
    ph = IC[1]
    om = IC[2]
    b = IC[3]

    Gam_sch=r**2.-2.*r

    dR_dr = ph
    dph_dr = (1./Gam_sch)*(2.*(1.-r)*ph+L*(l*(l+1.))-om**2.*r**4.*L/Gam_sch+(m**2.+lam*L**2.)*r**2.*L)
    dom_dr = b
    db_dr = 0.
    return [dR_dr,dph_dr,dom_dr,db_dr]

然后我尝试不同的“om”值并询问是否满足我的边界条件。p_sch 是我模型的参数。一般来说，我想做的事情要复杂一些，一般来说，我需要更多的参数，而不是在大规模的情况下。但是我需要从最简单的开始，这就是我在这里要问的

p_sch = (1,1,0,0) #[rho_c,m,lam,l], lam and l are for a more complicated case 
ep = 0.2
ep_r = 0.01
r_end = 500
n_r = 500000
n_omega = 1000
omega = np.linspace(p_sch[1]-ep,p_sch[1],n_omega)
r = np.linspace(2+ep_r,r_end,n_r)
tol = 0.01
a = 0

for j in range(len(omega)): 
    print('trying with $omega =$',omega[j])
    omeg = [omega[j]]
    ini = init_sch(omeg)
    Y = odeint(F_sch,ini,r,p_sch,mxstep=50000000)
    print Y[-1,0]
    #This is my boundary condition that I am implementing. Basically this should be my conditions at "infinity
    if abs(Y[-1,0]*((p_sch[1]**2.-Y[-1,2]**2.)**(1/2.)+1./(r[-1]))+Y[-1,1]) < tol: 
        print(j,'times iterations in omega')
        print("R'(inf)) = ", Y[-1,0])        
        print("\omega",omega[j])
        omega_1 = [omega[j]] 
        a = 10
        break           
    if a > 1:
        break

基本上我在这里要做的是求解给出不同初始条件的方程组，并找到应该接近我的边界条件的“a=”（或代码中的“om”）的值。我需要这个，因为在此之后我可以将这样的初始客人提供给割线方法并尝试为“a”获得最佳价值。但是，总是在运行此代码时，我有不同的解决方案，当然，这是一种我不感兴趣的行为。我正在尝试相同但考虑到 scipy.integrate.solve_vbp，但是当我运行以下代码时：

from IPython import get_ipython
get_ipython().magic('reset -sf')
import numpy as np
import matplotlib.pyplot as plt
from math import *
from scipy.integrate import solve_bvp

def bc(ya,yb,p_sch):
    m = p_sch[1]
    om = p_sch[4]
    tol_s = p_sch[5]
    r_end = p_sch[6]

    return np.array([ya[0]-1,yb[0]-tol_s,ya[1],yb[1]+((m**2-yb[2]**2)**(1/2)+1/r_end)*yb[0],ya[2]-om,yb[2]-om,ya[3],yb[3]])

def fun(r,y,p_sch):
    rho_c = p_sch[0]
    m = p_sch[1]
    lam = p_sch[2]
    l = p_sch[3]

    L = y[0]
    ph = y[1]
    om = y[2]
    b = y[3]

    Gam_sch=r**2.-2.*r

    dR_dr = ph
    dph_dr = (1./Gam_sch)*(2.*(1.-r)*ph+L*(l*(l+1.))-om**2.*r**4.*L/Gam_sch+(m**2.+lam*L**2.)*r**2.*L)
    dom_dr = b
    db_dr = 0.*y[3]
    return np.vstack((dR_dr,dph_dr,dom_dr,db_dr))

eps_r=0.01
r_end = 500
n_r = 50000
r = np.linspace(2+eps_r,r_end,n_r)
y = np.zeros((4,r.size))
y[0]=1
tol_s = 0.0001
p_sch= (1,1,0,0,0.8,tol_s,r_end)


sol = solve_bvp(fun,bc, r, y, p_sch)

我得到这个错误：ValueError：bcreturn 应该有形状（11，），但实际上有（8，）。ValueError：bcreturn 应该具有形状 (11,)，但实际上具有 (8,)。