为Scipy的“odeint”边界条件指定不同的时间点？

Question

我试图数值求解一个由两个非线性常微分方程组成的系统。我正在使用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))

Answer 1

odeint求解initial value problem。你描述的问题是两点boundary value problem。为此，可以使用^{}

您还可以查看^{}和{a5}，尽管看起来{}已经很久没有更新了。在

例如，下面是如何使用scipy.integrate.solve_bvp。我改变了参数，这样溶液就不会衰减得这么快，频率也会更低。当b=0.25时，衰减速度足够快，对于ω（0）=0且|θ（0）|为1阶的所有解，θ（100）≈0。在

函数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修改为

^{pr2}$

另一种解决方案如下图所示：

在这个解中，摆锤几乎在倒立位置（result.y[0, 0] = 3.141592653578858）开始几乎。它开始下降非常缓慢，逐渐下降得更快，在t=100时达到直线下降位置。在

平凡解θ（t）∠0和ω（t）∠0也满足边界条件。在