我正在尝试对包含两个非线性 ODE 的系统进行数值求解。我正在使用 Scipy 的 odeint
函数。 odeint
需要一个指定初始条件的参数 y0
。然而,它似乎假设 y0
的初始条件在同一时间点开始(即两个条件都在 t=0)。在我的例子中,我想指定为不同时间指定的两个不同的边界条件(即 omega(t=0) = 0,theta(t=100) = 0)。我似乎不知道该怎么做,非常感谢任何帮助!
下面是一些示例代码:
from scipy.integrate import odeint
def pend(y, t, b, c):
theta, omega = y
dydt = [omega, -b*omega - c*np.sin(theta)]
return dydt
b = 0.25
c = 5.0
t = np.linspace(0, 100, 101)
# I want to make these initial conditions specified at different times
y0 = [0, 0]
sol = odeint(pend, y0, t, args=(b, c))
最佳答案
odeint
解决了 initial value problem .您描述的问题是两点boundary value problem .为此,您可以使用 scipy.integrate.solve_bvp
你也可以看看 scikits.bvp1lg
和 scikits.bvp_solver
,虽然看起来 bvp_solver
已经很久没有更新了。
例如,这里是您可以如何使用 scipy.integrate.solve_bvp
。我更改了参数,因此溶液不会衰减得那么快并且频率较低。 b = 0.25 时,衰减足够快,对于所有 ω(0) = 0 和 |θ(0)| 的解,θ(100) ≈ 0数量级为 1。
函数 bc
将在 t=0 和 t=100 时传递 [θ(t), ω(t)] 的值。它必须返回两个值,它们是边界条件的“残差”。这只是意味着它必须计算必须为 0 的值。在您的情况下,只需返回 y0[1]
(即 ω(0))和 y1[0]
(即 θ(100))。 (如果 t=0 处的边界条件为 ω(0) = 1
,则 bc
返回值的第一个元素将为 y0[ 1] - 1
。)
import numpy as np
from scipy.integrate import solve_bvp, odeint
import matplotlib.pyplot as plt
def pend(t, y, b, c):
theta, omega = y
dydt = [omega, -b*omega - c*np.sin(theta)]
return dydt
def bc(y0, y1, b, c):
# Values at t=0:
theta0, omega0 = y0
# Values at t=100:
theta1, omega1 = y1
# These return values are what we want to be 0:
return [omega0, theta1]
b = 0.02
c = 0.08
t = np.linspace(0, 100, 201)
# Use the solution to the initial value problem as the initial guess
# for the BVP solver. (This is probably not necessary! Other, simpler
# guesses might also work.)
ystart = odeint(pend, [1, 0], t, args=(b, c,), tfirst=True)
result = solve_bvp(lambda t, y: pend(t, y, b=b, c=c),
lambda y0, y1: bc(y0, y1, b=b, c=c),
t, ystart.T)
plt.figure(figsize=(6.5, 3.5))
plt.plot(result.x, result.y[0], label=r'$\theta(t)$')
plt.plot(result.x, result.y[1], '--', label=r'$\omega(t)$')
plt.xlabel('t')
plt.grid()
plt.legend(framealpha=1, shadow=True)
plt.tight_layout()
plt.show()
这是结果图,您可以在其中看到 ω(0) = 0 和 θ(100) = 0。
请注意,边值问题的解不是唯一的。如果我们修改创建 ystart
为
ystart = odeint(pend, [np.pi, 0], t, args=(b, c,), tfirst=True)
找到了一个不同的解决方案,如下图所示:
在此解决方案中,钟摆开始时几乎处于倒置位置 (result.y[0, 0] = 3.141592653578858
)。它开始下降非常缓慢;它逐渐下降得更快,并在 t = 100 时到达直线下降位置。
平凡解θ(t) ≡ 0 和 ω(t) ≡ 0 也满足边界条件。
关于python - 为 Scipy 的 "odeint"的边界条件指定不同的时间点?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/51733696/