scipy solve_ivp初始状态字段的参数传入顺序逻辑是什么

问题背景 #

范德波尔（van der Pol）方程表达式为：
$$\frac{d^2x}{dt2} - \mu(1-x^2)\frac{dx}{dt} + x = 0$$
这是二阶常微分方程，无法直接传入Scipy的solve_ivp求解，必须先降阶为一阶方程组，降阶结果为：

• $\frac{dx}{dt} = v$
• $\frac{dv}{dt} = \mu(1-x^2)v - x$
  参考实现代码如下：

import sympy as sp
import matplotlib.pyplot as plt
import numpy as np 
from scipy.integrate import solve_ivp

def vdp_derivatives(t, y):
    x = y[0]
    v = y[1]
    return [v, mu*(1 - x*x)*v - x] 

t = np.linspace(0, 10, 500)
mu = 1
sol = solve_ivp(fun=vdp_derivatives, t_span=[t[0], t[-1]], y0=[1,0], t_eval=t) 

y0顺序的核心逻辑 #

solve_ivp本身不会强制规定y0里元素的排列顺序，顺序完全由你在导数计算函数里定义的状态向量索引规则决定，要求三者严格对齐：从y中取变量的顺序、返回导数值的顺序、y0传参的顺序。
拿上面的范德波尔方程代码举例：

• 导数函数里明确写了x = y[0]、v = y[1]，相当于你自己定了规则：状态数组y的第0位存位移$x$，第1位存速度$v$（也就是$dx/dt$）
• 函数返回值是[v, mu*(1 - x*x)*v - x]，也是按照「x的一阶导、v的一阶导」的顺序排列，和前面的状态索引一一对应
• 因此传y0的时候必须遵守同一个规则：[x初始值, v初始值]，示例里的[1,0]就代表初始时刻x=1、速度v=0。
  调换x和v的初始值顺序代码不会报错，是因为语法层面完全合法，但实际上你给方程设置了完全不同的初始条件（变成初始时刻x=0、v=1），解出来的曲线偏移是必然结果，本质是初始条件给错了，不是代码逻辑有语法问题。

适用于所有高阶常微分方程的通用方法 #

不管是几阶的常微分方程，降阶求解时按下面的步骤做就不会出顺序错误：

1. 先完成降阶替换：对n阶常微分方程，依次定义状态变量$y_0=x, y_1=\frac{dx}{dt}, y_2=\frac{d^2x}{dt2},...,y_{n-1}=\frac{d^{n-1}x}{dt{n-1}}$，把原方程拆成n个一阶微分方程
2. 写导数函数时，严格按照上面定义的状态变量顺序，从输入的y数组按索引取对应变量，再按完全相同的顺序返回每个状态变量的一阶导数值
3. 构造y0时，和你定义的状态变量顺序一一对应，依次传入每个状态变量在初始时刻的取值即可。

说白了状态向量的顺序是自定义的，没有强制要求，只要取变量、返回导数、传初始值三个环节的顺序完全匹配，求解结果就是正确的。

内容的提问来源于stack exchange，提问作者Gary Bollenbach