由于您有二阶 ODE，您必须提供 2 个初始条件，例如：

$$y'(0) = y'_0 \qquad y(0) = y_0\tag{1}$$

您可以做的一件事是转换原始二阶方程： 到一阶 ODE 系统 初始条件：

$$\left(A+By-Cy'^{-1/3}\right)y''=-\left(By'^2+C\right) \tag{2}$$ $$\left\{\begin{array}{ll}% v' =& -\frac{Bv^2+C}{A+By-Cv^{-1/3}}\\ y'=&v \end{array}\right. \tag{3}$$ $$v(0) = y'_0 \qquad y(0) = y_0$$

系统可以放入一个通用向量方程： 是 for and and。因此： $(2)$$\vec{x} = [x_1,x_2]^{T}$$x_1=v$$x_2=y$$\vec{f} = \left[-\frac{Bv^2+C}{A+By-Cv^{-1/3}},\,v\right]^T$

$$\frac{d \vec{x}}{dt} = \vec{f}(\vec{x})$$

Yu 可以使用任何支持刚性方程的 ODE 求解器（隐式数值方法）通过 matlab 轻松求解。

您还必须指定要使用的方案。

假设您要使用简单的 Backward-Euler (BE) 方案。对于 BE 方案，您的控制方程变为

$\left(A+By_{n+1}-C(y'_{n+1})^{-1/3}\right)y''_{n+1}=-\left(B(y'_{n+1})^2+C\right) \tag{1}$ with $y'_{n+1}=\frac{y_{n+1}-y_n}{\Delta t} \tag{2}$

$y''_{n+1}=\frac{y'_{n+1}-y'_n}{\Delta t} \tag{3}$

使用 (3) 中的 (2)， $y''_{n+1}=\frac{1}{\Delta t^2}[y_{n+1}-y_n] - \frac{y'_n}{\Delta t} \tag{4}$

现在，使用 (1) 中的 (2) 和 (4)，您将得到如下形式的非线性方程

$f(y_n,y'_n, y_{n+1})=0 \tag{5}$ 可以使用迭代方案求解，例如 Newton-Raphson 方案。得到后，您可以分别从 (2) 和 (4) 计算速度和加速度。$y_{n+1}$$y_{n+1}$

您可以对其他方案遵循相同的过程。对于不同的方案，等式 (1)、(2)、(3) 和 (4) 将有所不同。