考虑常微分方程

欧拉曾给过一个算法，这个算法是所有数值求解常微分方程的算法中最简单最直观的。即

这个算法可以想象精度非常差。这点可以通过考虑一类特殊情况非常明显地看到。假设f仅是x的函数，这时方程可以直接积出来，

其实就是要对f函数做个数值积分，也就是要求下图中f曲线下的面积。

但是如果用euler法，我们求的是图中阴影部分矩形的面积。这个数值积分方法甚至连梯形法都不如！梯形法的精度至少还是步长h的平方的量级，这里的矩形法的精度很明显只是h的量级。

我们试着用euler法求解一维谐振子的运动

相关程序如下：

x0 = 1; p0 = 2; h = 0.05;

tlist = 0 : h : 50;

xlist = zeros(1, length(tlist));

plist = zeros(1, length(tlist));

xlist(1) = x0; plist(1) = p0 ;

for s = 2 : length(tlist)

   xlist(s) = xlist(s-1) + h * plist(s-1);

   plist(s) = plist(s-1) - h* xlist(s-1);

end

h1 = figure;

plot(tlist, xlist, tlist, plist)

xlabel('t'); ylabel('x/p')

h2 = figure;

plot(xlist, plist)

xlabel('x'); ylabel('p')

输出结果如下。误差明显，振幅随时间迅速增大。

Euler法很差，但是其改进的余地很大，改进也很容易。比如我们有Heun方法：

很容易证明，euler法下，y_{n+1}对h的泰勒展开的前两项是对的（也就是到h的一次项），而heun法下，其对h的泰勒展开的前三项是对的（即到h的二次项）。下面我们用heun法重新求解谐振子的运动：

x0 = 1; p0 = 2; h = 0.1;

tlist = 0 : h : 50;

xlist = zeros(1, length(tlist));

plist = zeros(1, length(tlist));

xlist(1) = x0;plist(1) = p0 ;

for s = 2 : length(tlist)

   tempx = xlist(s-1) + h * plist(s-1);

   tempp = plist(s-1) - h* xlist(s-1);

   xlist(s) = xlist(s-1) + 0.5 * h * (plist(s-1) + tempp);

   plist(s) = plist(s-1) - 0.5* h * (xlist(s-1) + tempx);

end

h1 = figure;

plot(tlist, xlist, tlist, plist)

xlabel('t');ylabel('x/p')

h2 = figure;

plot(xlist, plist)

xlabel('x');ylabel('p')

结果见如下二图。可见，虽然步长h扩大了一倍，误差却小得多，x和p的振幅未见显著变化，系统在相空间的轨迹闭合为一个圆。

欧拉法和霍因法都是系统的Runge-Kutta法的特殊例子，并且是最简单的例子。以后我们会介绍精度更高的RK法。这类方法基本的思想是，重现尽可能多的taylor级数系数。

作业1：证明下面的方法也具有平方精度

作业2： 以谐振子作为例子其实是存在一个潜在问题的。这里我们有一个两分量的方程组，而不是一个单分量的方程。对单分量方程发展的算法并不能任意推广到多分量情况，不过就这个例子不存在这个问题，试证明之。