欧拉法求解微分方程

by Conmajia
2014
原文写于2014年，用C# 完成，现在改为JavaScript. 欧拉和中心差分逼近，是最朴素的想法，但代数精度太低，而龙格库塔的稳定性又是个问题. 总之只能用来计算普通的东西，高大上问题一般还是使用广义 $\alpha$，Wilson-$\Theta$ 之类的，利用系数调节人工阻尼，让低频的部分少受算法影响，高频阻尼增大，增加微分方程的稳定性，又不影响算法的精度.

当然，现在大把的lib、app，你完全可以用py、mtlb之类的搞定，不需要搞懂任何原理.

Math.NET计划是用于.NET的数学运算库，包括了数值计算库Iridium.NET，非常经典（且古老）.

工程中有时候需要解微分方程，比如这种：

$$y'=y-\cfrac{2x}{y}$$

$y(0)=1$. 两边积分，得到：

$$\begin{align*} \int_{x_n}^{x_{n+1}}y'\mathrm{d} x&=\int_{x_n}^{x_{n+1}}\left(y-\cfrac{2x}{y}\right)\mathrm{d} x \\ y_{n+1}-y_n&=\int_{x_n}^{x_{n+1}}\left(y-\cfrac{2x}{y}\right)\mathrm{d}x \\ y&=\sqrt{1+2x} \end{align*}$$

工程上一般不需要精确解，求个数值就好，通常采用差分近似代替微分，最简单最朴素的办法是：

$$\tag{1} y_{n+1}=y_n+\Delta x f(x_n,y_n)\quad n=0,1,2,\cdots$$

推导很简单，对 $x_n$，有：

$$y'_n=f(x_n,y_n).$$

对 $y'(x_n)$ 有：

$$\tag{2} y'_n\approx \cfrac{y_{n+1}-y_n}{\Delta x}.$$

(2)的左边叫微商，右边叫差商，$\Delta x=\left|x_{n+1}-x_n\right|$.

代回(2)，有：

$$\begin{align*} y'_n=f(x_n,y_n) \Rightarrow\cfrac{y_{n+1}-y_n}{\Delta x} &\approx f(x_n,y_n) \\ y_{n+1}-y_n &\approx \Delta x f(x_n,y_n) \\ y_{n+1} &\approx y(x_n)+\Delta x f(x_n,y_n). \end{align*}$$

于是，(1)成了：

$$\tag{3} y_{n+1}=y_n+\Delta x\left(y_n-\cfrac{2x_n}{y_n}\right).$$

用myFn函数表示余项 $y_n-\dfrac{2x_n}{y_n}$：

function myFn(x, y) {
  return y - 2 * x / y;
}

Euler可以这样实现：

function calculate(x0, y0, delta, xn) {
  var yn;
  while(x0 < xn) {
    yn = y0 + delta * myFn(x0, y0);
    y0 = yn;
    x0 = x0 + delta;
  }
  return yn;
}

现在可以开始试验. 理论上：

$$y=\sqrt{1+2x}\Rightarrow y(1)=\sqrt{3}\approx 1.7321$$

$\sqrt{3}$ 是方程的精确解，程序的目标是通过计算，得到尽量接近精确解的结果.

取 $x_0=0$，$y_0=1$，$\Delta x=0.1$，程序计算结果为：

ans = 1.784771

误差 $\varepsilon_1=3.04\%$. 现在减小 $\Delta x$，比如 $\Delta x=0.0001$，计算结果为：

ans = 1.732112

$\varepsilon_2=0.0007\%$，比较接近真值了.

尽管这个方法可以求到比较精确的解，但是 $\Delta x$ 太小的话，显然while会执行很多次 $N_\mathrm{while}\prop\dfrac{1}{\Delta x}$，效率低下.

引入定积分的梯形公式：

$$\int_{x_n}^{x_{n+1}}f(x,y)\mathrm{d}x\approx\cfrac{\Delta x}{2}\left[f(x_n,y_n)+f(x_{n+1},y_{n+1})\right]$$

(3)可以写成：

$$\tag{4} y_{n+1}\approx y_n+\cfrac{\Delta x}{2}\left[f(x_n,y_n)+f(x_{n+1},y_{n+1})\right].$$

(3)和(4)的特点是，前者速度快，精度低；后者速度慢，精度高. 所以对于粗算，可以使用：

$$y'_{n+1}=y_n+\Delta x f(x_n,y_n)$$

精算，可以用：

$$y_{n+1}=y_n+\cfrac{\Delta x}{2}\left[f(x_n,y_n)+f(x_{n+1},y'_{n+1})\right].$$

结合起来，先通过粗算得到 $y'_{n+1}$ 的近似值，再进行精算，得到高精度的最终解. 这是一个比较均衡的折衷，既保证了计算结果的准确度，又保证了计算效率. 下面是改进后的代码：

function calculate(x0, y0, delta, xn) {
  var yp, yc;
  while(x0 < xn) {
    yp = y0 + delta * myFn(x0, y0);
    yc = y0 + delta * myFn(x0 + delta, yp);
    y0 = 1 / 2 * (yp + yc);
    x0 = x0 + delta;
  }
  return y0;
}

取 $x_0=0$，$y_0=1$，$\Delta x=0.1$：

ans = 1.737867

$\varepsilon_3=0.33\%$.

和 $\varepsilon_1$ 相比，在 $\Delta x$ 相同——意味着执行时间相近——的情况下，精度提高了10倍.

手工求微分方程基本是这么个思路. 当然，更多的时候，做个伸手党其实很不错的，大把的现成玩意儿可以用，我干嘛还要费那劲呢？

The End. $\Box$