为Scipy的"odeint“的边界条件指定不同的时间点？
我正在尝试数值求解一个由两个非线性常微分方程组成的系统。我使用的是Scipy的odeint函数。odeint需要指定初始条件的参数y0。然而，它似乎假设y0的初始条件开始于相同的时间点(即两个条件都在t=0)。在我的例子中，我想指定两个不同的边界条件，这两个条件是为不同的时间指定的(即ω(t=0)= 0，θ(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))回答 1
Stack Overflow用户
发布于 2018-08-08 04:19:36
odeint解决了问题。您描述的问题是两点。为此，您可以使用scipy.integrate.solve_bvp
你也可以看看scikits.bvp1lg和scikits.bvp_solver，尽管看起来bvp_solver已经很长时间没有更新了。
例如，下面是使用scipy.integrate.solve_bvp的方法。我更改了参数，以便解决方案不会衰减得太快，并且具有较低的频率。当b= 0.25时，衰减足够快，使得θ(100)对所有解≈0，其中ω(0) =0并且|θ(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也满足边界条件。
https://stackoverflow.com/questions/51733696
复制相似问题
腾讯云开发者
Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有
深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 
粤公网安备44030502008569号
腾讯云计算（北京）有限责任公司 京ICP证150476号 | 京ICP备11018762号