Odeint, метод стрельбы и граничные условия в Python - PullRequest
Новая
       23

Odeint, метод стрельбы и граничные условия в Python

0 голосов
/ 01 мая 2018

Я работал с 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). Я попытался применить метод стрельбы, чтобы найти начальные условия для альфы, которые лучше всего соответствуют моим граничным условиям. Мне также нужно найти омегу, которая лучше всего соответствует граничным условиям для 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

Однако, когда я запускаю его, я получаю: error: функция и ее якобиан должны быть вызываемыми функциями.

Может ли кто-нибудь помочь мне понять, в чем моя ошибка?

С уважением,

Луис Падилья.

Ответы [ 2 ]

0 голосов
/ 03 мая 2018

Я исправил свой код, и теперь он дает мне некоторые результаты. Однако, когда я запускаю его, я получаю некоторые предупреждения, которые я не знаю, как решить это. Может ли кто-нибудь помочь мне решить это? В основном мой код таков: 

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

Как я мог решить проблему?

С уважением,

Луис Падилья.

0 голосов
/ 01 мая 2018

Для начала, первым аргументом odeint является ваша производная функция system_DE. Просто передайте его имя, без скобок и аргументов. Odeint с внутренним вызовом и предоставлением аргументов.

...