改进euler方法 c语言,科学网—计算方法：Euler法及其改进 - 张江敏的博文

考虑常微分方程

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

这个算法可以想象精度非常差。这点可以通过考虑一类特殊情况非常明显地看到。假设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： 以谐振子作为例子其实是存在一个潜在问题的。这里我们有一个两分量的方程组，而不是一个单分量的方程。对单分量方程发展的算法并不能任意推广到多分量情况，不过就这个例子不存在这个问题，试证明之。

转载本文请联系原作者获取授权，同时请注明本文来自张江敏科学网博客。

链接地址：http://blog.sciencenet.cn/blog-100379-1094683.html

上一篇：彩虹的数学与物理

下一篇：计算方法：Euler法，Heun法，RK法的精度