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

问题描述

我正在尝试数值求解两个非线性 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))

标签： pythonscipyode

odeint解决了初值问题。您描述的问题是两点边值问题。为此，您可以使用scipy.integrate.solve_bvp

你也可以看看scikits.bvp1lgand 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”的边界条件指定不同的时间点？

问题描述

解决方案

推荐阅读