Python中的Odeint、拍摄方法和边界条件

Question

我一直在使用 odeint 和边界条件。基本上我要做的是解决这个图中给出的微分方程1

在我的代码中，R=R，ph = Phi，al = alpha，a = a，m = m，l = l 和 om = omega。我试图实现的初始条件是 R(0)=O(r^l); Phi(0)=O(r^{l-1}) 如果 l/=0 并且 Phi(0)=O(r) 如果 l=0； a(0) = 1 和 a(inf)=1/alpha(inf) （另外我需要 R(inf)=0）。我尝试应用拍摄方法来找到与我的边界条件最匹配的 alpha 初始条件。我还需要找到与 R 在无穷远处的边界条件最匹配的欧米茄。我写的代码如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import time

start = time.clock()
def system_DE(IC,p,r):
    l = p[0]
    m = p[1]
    om = p[2]

    R = IC[0]
    ph = IC[1]
    a = IC[2]
    al = IC[3]

    dR_dr = ph
    da_dr = a*((2*l+1)*r/2*(om**2*a**2*R**2/al**2+ph**2+l*(l+1)*a**2*R**2/r**2+m**2*a**2*R**2)-(a**2-1)/(2*r))
    dal_dr = al*(da_dr/a-l*(l+1)*(2*l+1)*a**2*R**2/r-(2*l+1)*m**2*a**2*r*R**2+(a**2-1)/r)
    dph_dr = -2*ph/r-dal_dr*ph/al+da_dr*ph/a-om**2*a**2*R/al**2+l*(l+1)*a**2*R/r**2+m**2*a**2*R

    return [dR_dr,da_dr,dal_dr,dph_dr]

def init(u,p,r):
    if p==0:
        return np.array([1,r,1,u])
    else:
        return np.array([r**l,l*r**(l-1),1,u])

l = 0
m = 1
ep = 0.3
n_om = 10
omega = np.linspace(m-ep,m+ep,n_om)
r = np.linspace(0.0001, 100, 1000)


niter = 100
u = 0
tol = 0.1
ustep = 0.01

p = np.zeros(3)
p[0] = l
p[1] = m
for j in range(len(omega)):
    p[2] = omega[j]
    for i in range(niter):
        u += ustep
        Y = odeint(system_DE(init(u,p[0],r[0]),p,r), init(u,p[0],r[0]), r)
        print Y[-1,2]
        print Y[-1,3]
        if abs(Y[len(Y)-1,2]-1/Y[len(Y)-1,3]) < tol:
            print(i,'times iterations')
            print("a'(inf)) = ", Y[len(Y)-1,2])
            print('y"(0) =',u)
            break
    if abs(Y[len(Y)-1,0]) < tol:
        print(j,'times iterations in omega')
        print("R'(inf)) = ", Y[len(Y)-1,0])
        break

但是，当我运行它时，我得到：错误：函数及其雅可比必须是可调用函数。

有人可以帮助我理解我的错误是什么吗？

问候，

路易斯·帕迪拉。

Answer 1

我修复了我的代码，现在它给了我一些结果。但是，当我运行它时，我收到了一些我不知道如何解决的警告。有人可以帮我解决吗？基本上我的代码是这样的：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import time

def system_DE(IC,r,l,m,om):    
    R = IC[0]
    ph = IC[1]
    a = IC[2]
    al = IC[3]

    dR_dr = ph
    da_dr = a*((2*l+1)*r/2*(om**2*a**2*R**2/al**2+ph**2+l*(l+1)*a**2*R**2/r**2+m**2*a**2*R**2)-(a**2-1)/(2*r))
    dal_dr = al*(da_dr/a-l*(l+1)*(2*l+1)*a**2*R**2/r-(2*l+1)*m**2*a**2*r*R**2+(a**2-1)/r)
    dph_dr = -2*ph/r-dal_dr*ph/al+da_dr*ph/a-om**2*a**2*R/al**2+l*(l+1)*a**2*R/r**2+m**2*a**2*R

    return [dR_dr,dph_dr,da_dr,dal_dr]

def init(u,p,r):
    if p==0:
        return np.array([1.,r,1.,u])
    else:
        return np.array([r**p,l*r**(p-1),1,u])


l = 0.
m = 1.
ep = 0.2
n_om = 30
omega = np.linspace(m-ep,m+ep,n_om)
r = np.linspace(0.001, 100, 1000)


niter = 1000
tol = 0.01
ustep = 0.01

for j in range(len(omega)):
    print('trying with $omega =$',omega[j])
    p = (l,m,omega[j])
    u = 0.001
    for i in range(niter):
        u += ustep
        ini = init(u,p[0],r[0])
        Y = odeint(system_DE, ini,r,p,mxstep=500000)
        if abs(Y[len(Y)-1,2]-1/Y[len(Y)-1,3]) < tol:
            break
    if abs(Y[len(Y)-1,0]) < tol and abs(Y[len(Y)-1,2]-1/Y[len(Y)-1,3]) < tol:
        print(j,'times iterations in omega')
        print(i,'times iterations')
        print("R'(inf)) = ", Y[len(Y)-1,0])        
        print("alpha(0)) = ", Y[0,3])
        print("\omega",omega[j])
        break

plt.subplot(2,1,1)
plt.plot(r,Y[:,0],'r',label = '$R$')
plt.plot(r,Y[:,1],'b',label = '$d R /dr$')
plt.xlim([0,10])
plt.legend()
plt.subplot(2,1,2)
plt.plot(r,Y[:,2],'r',label = 'a')
plt.plot(r,Y[:,3],'b', label = '$alpha$')
plt.xlim([0,10])
plt.legend()
plt.show()

但是当我运行它时，我得到了这个：

lsoda--  warning..internal t (=r1) and h (=r2) are
   such that in the machine, t + h = t on the next step  
   (h = step size). solver will continue anyway
  in above,  r1 =  0.1243782486482D+01   r2 =  0.8727680448722D-16

我该如何解决这个问题？

问候，

路易斯·帕迪拉。

Answer 2

首先，odeint 的第一个参数是您的派生函数 system_DE。只传递它的名字，没有括号或参数。 Odeint 在内部调用它并提供参数。