常微分方程(组)的数值解

                     

贡献者: addis

预备知识 天体运动的简单数值计算

   数值解微分方程或微分方程组时,一般需要先把微分方程(组)化成一阶微分方程组。例如在 “弹簧振子受迫运动的简单数值计算” 中列出的微分方程为

\begin{equation} m y'' = \alpha y' - ky + f(t)~. \end{equation}
新增变量,令 $v = y'$,则可变为一阶微分方程组
\begin{equation} \begin{cases} v' = [-\alpha u - ky + f(t)]/m\\ y' = v \end{cases}~. \end{equation}
方程组中,$t$ 是自变量,$y$ 和 $u$ 是 $t$ 的函数。给出某时刻的 $y(t), u(t), t$ 就可以通过方程组求出 $u'(t)$ 和 $y'(t)$ 的数值。

   又例如在 “天体运动的简单数值计算” 中,列出的二阶微分方程组为

\begin{equation} \left\{\begin{aligned} x'' &= -\frac{GM}{(x^2+y^2)^{3/2}}x\\ y'' &= -\frac{GM}{(x^2+y^2)^{3/2}}y \end{aligned}\right. ~.\end{equation}

   新增变量,令 $v_x = x'$,$v_y = y'$,上式也可变为一阶微分方程组

\begin{equation} \begin{cases} x' = v_x\\ v'_x = -GMx/(x^2 + y^2)^{3/2}\\ y' = v_y\\ v'_y = -GMy/(x^2 + y^2)^{3/2} \end{cases}~. \end{equation}
该式中同样 $t$ 是自变量,其他都是 $t$ 的函数,给出某时刻的 $x, y, v_x, v_y, t$ 就可以由该方程组求出各函数的一阶导数。

   在以上两个例子中,我们使用了一种较为原始的方法(微分近似)。这种方法相当于把某时刻 $t$ 的各函数值代入一阶微分方程组,得到 $t$ 时刻各函数的一阶导数,再通过微分近似由这些一阶导数来计算其 $[t, t + \Delta t]$ 时间内的增量,得到 $t +\Delta t$ 时刻的各函数值,再代入一阶方程组得到 $t +\Delta t$ 时刻各函数的一阶导数,如此一直循环,得到各函数每隔 $\Delta t$ 时间的值。这种方法叫做欧拉法。对一阶常微分方程 $y'(t) = f(y, t)$,令 $h = \Delta t$,$t_n = t_0 + nh$,$y_n = y(t_n)$ 欧拉法可以表示为

\begin{equation} y_{n+1} = y_n + h f(y_n, t_n)~. \end{equation}

   以下介绍三种更精确的算法,在步长($\Delta t$)相同时它们的精确度递增。注意虽然大多数问题中自变量是时间,但这些方法适用于大多数单个自变量的微分方程(组)。


致读者: 小时百科一直以来坚持所有内容免费无广告,这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。 因此,我们请求广大读者热心打赏 ,使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 20 元,我们一周就能脱离亏损, 并在接下来的一年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获取知识, 我们在此表示感谢。

                     

友情链接: 超理论坛 | ©小时科技 保留一切权利