使用 ODEINT 的二阶耦合 ODE

问题描述

我是用 python 解决耦合 ODE 的新手，我想知道我的方法是否正确，目前这段代码输出的图形与预期的输出完全不同。这些是我试图解决的方程：

enter image description here

这是我正在使用的代码（对于函数 f_gr、f_sc_phi 和 f_gTheta，您可以输入任何常量值）

import Radial as rd
import ScatteringAzimuthal as sa
import PolarComponent as pc
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt


#gamma for Now set to 1
g_mm = 1

def f(u,t):
    #y1 = thetadot ::  y2 = phidot  ::  y3 = cdot
    rho,theta,y1,phi,y2,c,y3 = u
    p = [y1,(pc.f_gTheta(theta,524.1+rho)/(c*np.cos(phi))-(g_mm*y1)+(2*y1*y2*np.tan(phi))-(2*y3*y1/c)),((sa.f_sc_phi(theta,524.1+rho/c))-(g_mm*y2)-(2*y3*y2/c)-(np.sin(phi)*np.cos(phi)*y2**2)),y3,(rd.f_gr(theta,524.1+rho)-(g_mm*y3)+(c*y2**2)+(c*(y1**2)*(np.cos(phi)**2))),phi]
    return p

time = np.linspace(0,10,100)
z2 = odeint(f,[0.1,np.pi/2,0.1,0.1],time)

rhoPl = z2[:,0]
thetaPl = z2[:,1]
phiPl = z2[:,3]

'''
plt.plot(rhoPl,time)
plt.plot(thetaPl,time)
plt.plot(phiPl,time)
plt.show()

'''

x = rhoPl*np.sin(thetaPl)*np.cos(phiPl)
y = rhoPl*np.sin(thetaPl)*np.sin(phiPl)
z = rhoPl*np.cos(thetaPl)

plt.plot(x,time)
plt.plot(y,time)
plt.plot(z,time)
plt.show()

当我将时间从 0.1 更改为 5 时，出现错误：

ODEintWarning：在此调用上完成了过多的工作（可能是错误的 Dfun 类型）。使用 full_output = 1 运行以获取定量信息。

关于如何改进此代码的任何想法，或者我的方法是否完全不正确？

Radial.py 的代码

import numpy as np
from scipy.special import spherical_jn
from scipy.special import spherical_yn
import sympy as sp
import matplotlib.pyplot as plt
R_r = 5.6*10**(-5) 
l = 720
n_w = 1.326
#k = 524.5/R_r
X_r = 524.5



# R is constant r is changing  
def f_gr(theta,x):
    f = ((sp.sin(theta))**(2*l-2))*(1+(sp.cos(theta))**2)
    b = (spherical_jn(l,n_w*x)*spherical_jn(l,n_w*x,True))+(spherical_yn(l,n_w*x)*spherical_yn(l,True))
    c = (spherical_jn(l,n_w*X_r)*spherical_jn(l,n_w*X_r,n_w*X_r)*spherical_yn(l,True))
    n = b/c
    f = f*n
    return f

ScatteringAzimuthal.py 代码

from scipy.special import spherical_jn,spherical_yn
import numpy as np
import matplotlib.pyplot as plt

l = 720
n_w = 1.326
n_p = 1.572
X_r = 524.5
R_r = 5.6*10**(-5)
R_p = 7.5*10**(-7)
k = X_r/R_r


def f_sc_phi(theta,x):
    f = (2/3)*(n_w**2)*((X_r**3)/x)*((R_P**3)/(R_r**3))*(((n_P**2)-(n_w**2))/((n_P**2)+(2*(n_w**2))))
    g = np.sin(theta)**(2*l-3)
    numerator = (l*(1+np.sin(theta))- np.cos(2*theta))\
                *((spherical_jn(l,n_w*x))+(spherical_yn(l,n_w*x)))
    denominator = ((spherical_jn(l,True))\
                   +(spherical_yn(l,True)))
    m = numerator/denominator
    final = f*g*m
    return final

PolarComponent.py 的代码

import numpy as np
from scipy.special import spherical_yn,spherical_jn
import matplotlib.pyplot as plt

l = 720
n_w = 1.326
X_r = 524.5 #this value is implemented in the ode file

#define dimensionless polar component
#X_r is radius,x is variable 
def f_gTheta(theta,x):
    bessel1 = (spherical_jn(l,n_w*x)) + \
              (spherical_yn(l,n_w*x))
    bessel2 = ((spherical_yn(l,True)) + \
              (spherical_yn(l,True)))*n_w*x
    bessels = bessel1/bessel2
    rest = (np.sin(theta)**(2*l-3))*((l-1)*(1+(np.cos(theta)**2)) \
                                     -((np.sin(theta)**2)*np.cos(theta)))
    final = rest*bessels
    return final

解决方法

这是一个 link，我非常喜欢用来模拟二阶颂歌。它有一个优化扭曲，因为它正在拟合模型以匹配模拟。它有几个 odeint 和 gekko 的例子。

python scipy scipy