Python ODE 系统四阶方程组

Question

我正在为以下颂歌系统的解决方案编写一个Python脚本。我使用 scipy

solve_bvp

库。

这是我的实现：

from scipy.integrate import solve_bvp
import numpy as np


def ode(r, y):
    Re=100
    return np.vstack((
#       y[0]    #  f        
        y[1],   # [f]' = [f']
        2/r**y[0]-2/3*Re/r**2*y[3]*y[0],   #[f']'=f''= ...
#       y[3]    #  g        
        y[4],   # [g]'=g'
        y[5],   # [g']'=g''
        y[6],   # [g'']'=g'''
        12/r**2*y[4]-24/r**3*y[4]+2*Re/(5*r**3)*(
            5*y[0]*(r*y[1]-2*y[0])+
            2*r*y[3]*y[6]+
            (r*y[4]-4*y[3])*y[5]
            -18/r*y[4]*y[5]
            +48/r**2*y[3]**2
            )
    ))


def bc(ya, yb):
    return np.array([
                     ya[0],     yb[0]-1,
                     ya[1],     yb[1],
                     ya[3],     yb[3], 
                     ya[4],     yb[4], 
                     ya[5],     yb[5],
                     ya[6],     yb[6]
                     ])

x_flow = np.linspace(0.7, 1, 10)
y_flow = np.ones((6, x_flow.shape[0]))

res_flow = solve_bvp(ode, bc, x_flow, y_flow, verbose=2)

目前求解器返回错误：

  y[6],   # [g'']'=g'''

IndexError: index 6 is out of bounds for axis 0 with size 6

使用调试器，在使用

np.asarray()

包装自定义代码时，在 _bvp.py 的第 647 行抛出错误。

Answer 1

问题出在索引上。您正在创建一个仅包含 6 个元素的 y 向量（因此索引可以从 0 到 5），但随后您调用

y[6]

，因此是第 7 个元素。然而，只需创建一个包含 7 个元素的 y_flow，例如

y_flow = np.ones((7, x_flow.shape[0]))

不起作用，因为 ODE 仅返回 6 个值，因此要么

ode

函数中 y 的索引有错误，要么缺少方程。我添加了一个虚拟方程，取消了 ODE 中第一个

y[0]

的注释，并在边界条件中出现错误，因为您传递的是 12 个值而不是 7 个值。

我认为您需要在将系统放入代码之前更好地重新构建系统，以便与索引保持一致。