我尝试使用 odeint
在Python中以y'' + ay' + by + c = 0
(二阶微分方程(的形式求解方程。据我了解,odeint
仅在y(0(= y1,y'(0(= y2的形式的初始条件下工作。而我的条件是边界:y'(0(= 0,y'(pi/4(= 0。是否有任何方法可以在这种情况下使用odeint
?这是我拥有的代码:
from scipy.integrate import odeint
import numpy as np
def lagrange(x, teta):
y = x[0]
dy = x[1]
xdot = [[],[]]
xdot[0] = dy
xdot[1] = (y + dy**2*y*(y**2 + dy**2)**(-1.5) + y*(y**2+dy**2)**(-0.5))/((y**2+dy**2)**(-0.5) - dy**2*(y**2+dy**2)**(-1.5))
return xdot
phi = np.linspace(0,np.pi/4,100)
U = odeint(lagrange, [u1_0, u2_0], phi)
官方文档指出,它解决了一阶ode-s的僵硬或非建立系统的初始值问题。因此,它无法解决边界价值问题。
您很可能需要实现有限差异方法,有限元方法或拍摄方法来解决问题。
我会亲自推荐有限的差异方法,因为它基本上包括添加3个对角矩阵(其中2个是tridiagonal(,设置解决方案向量并使用线性求解器。
。如果您需要设置一些帮助,我不久前就编写了一个扩散 - 添加求解器,您可以向其移植并扩展以包含反应术语(应该是b*Identity matrix(。
根据您的价值,您也可能需要实现人工扩散(如果Peclet号码太大(。