正面教材分享 http://blog.sciencenet.cn/u/wdlang 70%的以色列人是无神论者,不过他们都相信上帝给了他们那块土地。这个世界经不起思考

博文

计算方法:Euler法及其改进

已有 2691 次阅读 2018-1-13 20:32 |个人分类:计算方法|系统分类:教学心得

考虑常微分方程

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

这个算法可以想象精度非常差。这点可以通过考虑一类特殊情况非常明显地看到。假设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法的精度

2 刘全慧 徐令予

该博文允许注册用户评论 请点击登录 评论 (2 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2018-10-20 10:46

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部