为 Scipy 的“odeint"的边界条件指定不同的时间点? [英] Specify different timepoints for boundary conditions of Scipy's "odeint"?

问题描述

我正在尝试对两个非线性 ODE 的系统进行数值求解.我正在使用 Scipy 的 odeint 函数.odeint 需要一个参数 y0 来指定初始条件.然而，它似乎假设 y0 的初始条件在同一时间点开始(即两个条件都在 t=0).就我而言，我想指定为不同时间指定的两个不同边界条件(即 omega(t=0) = 0, theta(t=100) = 0).我似乎无法弄清楚如何做到这一点，非常感谢您的帮助！

下面的一些示例代码:

from scipy.integrate import odeintdef 挂起(y, t, b, c):θ, 欧米茄 = ydydt = [omega, -b*omega - c*np.sin(theta)]返回 dydtb = 0.25c = 5.0t = np.linspace(0, 100, 101)# 我要使这些初始条件在不同的时间指定y0 = [0, 0]sol = odeint(pend, y0, t, args=(b, c))

解决方案

odeint 解决了

注意边值问题的解不是唯一的.如果我们修改创建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 也满足边界条件.

I am trying to numerically solve a system of two nonlinear ODEs. I am using Scipy's odeint function. odeint requires an argument y0 which specify the initial conditions. However, it seems to assume that the initial conditions of y0 begin at the same timepoint (i.e both conditions are at t=0). In my case, I want to specify two different boundary conditions that are specified for different times (i.e. omega(t=0) = 0, theta(t=100) = 0). I can't seem to figure out how to do this and would greatly appreciate any help!

Some example code below:

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 solves the initial value problem. The problem that you describe is a two-point boundary value problem. For that, you can use scipy.integrate.solve_bvp

You could also take a look at scikits.bvp1lg and scikits.bvp_solver, although it looks like bvp_solver hasn't been updated in a long time.

For example, here is how you could use scipy.integrate.solve_bvp. I changed the parameters so the solution does not decay so fast and has a lower frequency. With b = 0.25, the decay is rapid enough that θ(100) ≈ 0 for all solutions where ω(0) = 0 and |θ(0)| is on the order of 1.

The function bc will be passed the values of [θ(t), ω(t)] at t=0 and t=100. It must return two values that are the "residuals" of the boundary conditions. That just means it must compute values that must be 0. In your case, just return y0[1] (which is ω(0)) and y1[0] (which is θ(100)). (If the boundary condition at t=0 had been, say ω(0) = 1, the first element of the return value of bc would be 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()

Here's the plot of the result, where you can see that ω(0) = 0 and θ(100) = 0.

Note that the solution to the boundary value problem is not unique. If we modify the creation ystart to

ystart = odeint(pend, [np.pi, 0], t, args=(b, c,), tfirst=True)

a different solution is found, as seen in the following plot:

In this solution, the pendulum starts out almost at the inverted position (result.y[0, 0] = 3.141592653578858). It starts to fall very slowly; gradually it falls faster, and it reaches the straight down position at t = 100.

The trivial solution θ(t) ≡ 0 and ω(t) ≡ 0 also satisfies the boundary conditions.

这篇关于为 Scipy 的“odeint"的边界条件指定不同的时间点?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！