使用odeint避免发散的解决方案？拍摄方法
我正在试着用Python解一个方程。基本上我想要做的就是解这个方程：
(1/x^2)*d(Gam*dL/dx)/dx)+(a^2*x^2/Gam-(m^2))*L=0这是Schwarzschild时空中大质量标量场的克莱因-戈登方程。假设我们知道m和Gam=x^2-2*x。我所知道的初始/边界条件是L(2+epsilon)=1和L(infty)=0。请注意，方程的渐近行为是
L(x-->infty)-->Exp[(m^2-a^2)*x]/x and Exp[-(m^2-a^2)*x]/x 然后，如果是a^2>m^2，我们就会有振荡解，而如果是a^2 < m^2，我们就会有散度解和衰减解。
我感兴趣的是衰减解，然而，当我试图解决上面的方程时，将其转换为一阶微分方程组，并使用打靶法，以便找到能给出我感兴趣的行为的"a“，我总是有一个发散的解。我认为这是因为odeint总是在寻找发散的渐近解。有没有办法避免或告诉odeint我对衰减解决方案感兴趣？如果没有，你知道我可以解决这个问题的方法吗？也许用另一种方法解我的微分方程组？如果是，采用哪种方法？
基本上，我要做的是为"a“添加一个新的方程系统。
(d^2a/dx^2=0, da/dx(2+epsilon)=0,a(2+epsilon)=a_0)为了让"a“成为一个常量。然后我考虑"a_0“的不同值，并询问我的边界条件是否满足。
耽误您时间，实在对不起。致以敬意，
Luis P。
回答 1
Stack Overflow用户
发布于 2018-06-17 12:39:06
考虑到渐近行为，我将无穷大处的值合并进来，这意味着我将有一个域和它的导数之间的关系。如果代码对你有帮助，我会帮你贴出来：
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]
#Here I ask if my asymptotic behavior is fulfilled or not. This should be basically my value 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：bc返回应该有形状(11，)，但实际上有(8，)。ValueError：bc返回应具有形状(11，)，但实际具有(8，)。
https://stackoverflow.com/questions/50885131
复制相似问题
腾讯云开发者
