数值分析知识内容 (37).pdf

《数值分析知识内容 (37).pdf》由会员分享，可在线阅读，更多相关《数值分析知识内容 (37).pdf（7页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、7.2 一阶初值问题的欧拉方法 将一阶初值问题转化为单步法(One-step Methods)，常用的离散化方法有：有限差分离散化方法、数值积分离散化方法、Taylor 展开离散化方法.首先介绍欧拉方法（Eulers Method）.7.2.1 欧拉方法 设)(),(11kkkkxyyxyy，给定初值)(00 xyy 及步长h，下面分别利用三种离散方法建立欧拉公式.1、利用一阶向前差分公式)(,()(1kkkkkxyxfxyhyy，得到欧拉公式 1,1,0),(1nkyxhfyykkkk.（8.4）2、利用 Taylor 展开方法，将)(xy在节点kx作 Taylor 展开，即有),()(,(

2、)()(2)()()(22hOxyxhfxyxyhxyhxyhxykkkkkkk （8.5）略去高阶项，则有欧拉公式1,1,0),(1nkyxhfyykkkk.3、数值积分方法.利用左矩形数值积分公式，对方程(8.3)两端在1,kkxx上积分)(,()(,()()(111kkkkxxkkxxxyxfdxxyxfxyxykk，则有欧拉公式1,1,0),(1nkyxhfyykkkk.当给定初值)(00 xyy 及步长h时，按公式（8.4）逐步递推，求出ky),2,1(nk值，称此方法为欧拉方法.欧拉方法（即折线法）的几何意义（见图 8-2）：在0,1的每个区间1,kkxx上，用过点),(kkyx且

3、斜率等于函数值),(kkyxf的折线，来逼近方程（8.3）的解)(xyy.00.10.20.30.40.50.60.70.80.9111.522.533.544.55Euler ApproximationExact Solution 图 8-2 欧拉方法几何意义 图 8-3 欧拉方法局部截断误差 通过几何直观图 8-3 来考察欧拉方法的精度.假设)(kkxyy，即顶点kP落在)(xyy 上，那么按欧拉方法做出的折线1kkPP就是)(xyy 过点kP的切线.从图 8-3 可知，这样求出的顶点1kP显著地偏离了原曲线，可见欧拉方法是相当粗糙的.为了分析数值计算公式的精度，引入局部截断误差的概念.如

4、果假设kx点的函数值是准确的，即)(kkxyy，则称111)(kkkyxyT为局部截断误差(Local Truncation Error).若某方法的局部截断误差为)(11pkhOT，则称该方法具有p阶精度.由 Taylor 展开公式（8.5）可知，)(2)(211kkkxyhyxy，即局部截断误差为)(21hOTk，所以欧拉方法具有 1 阶精度，是一阶方法.7.2.2 后退的欧拉方法 利用右矩形数值积分公式，对方程(8.3)两端在1,kkxx上积分)(,()(,()()(11111kkkkxxkkxxxyxfdxxyxfxyxykk，则有后退的欧拉公式 1,1,0),(111nkyxhfyy

5、kkkk.（8.6）与欧拉公式相比，后退的欧拉公式(Backward Euler One-step Method)有明显的区别，前者是关于1ky的一个直接的计算公式，称为显式的（Explicit）；后者公式的右端含有未知的1ky，它实际上是关于1ky的一个函数方程，称为隐式的（Implicit）.Tk+1 Pk+1 Pk y(x)O y x【注】（1）显式与隐式各有优点.在计算时显式算法比隐式方便，但考虑到数值稳定性等因素时，隐式算法比显式算法优越.（2）隐式公式通常用迭代法求解，而迭代过程的实质是逐步显式化.（3）后退的欧拉公式与欧拉公式的误差相似，也是一阶方法.这是因为局部截断误差)(,(

6、)()()(111111kkkkkkkxyxhfxyxyyxyT)()()()(2)()(22 kkkkkkxyhxyhxyxyhxyhxy )(22kxyh.7.2.3 梯形方法 在数值积分离散化的方法中，若采用梯形求积公式，则得到 )(12-)(,()(,(2)()(3111kkkkkkkxyhxyxfxyxfhxyxy.若以)(),(11kkkkxyyxyy，并略去高阶项，则得到梯形公式（隐式方程）),(),(2111kkkkkkyxfyxfhyy.（8.7）它的局部截断误差为)(31hOTk，所以梯形方法具有 2 阶精度，是二阶方法.显然梯形方法是隐式的(Implicit Trapez

7、oidal Method)，同后退的欧拉方法一样，可以用迭代法求解.利用欧拉方法提供初值，则梯形方法的迭代公式(Trapezoidal With Iteration)为),()0(1kkkkyxhfyy 1,2,),(),(2)1(11)(1iyxfyxfhyyikkkkkik.（8.8）迭代过程的实质是逐步显式化，这样就可以求出1ky的近似解序列)(1iky.分析迭代过程（8.8）的收敛性.若),(yxf满足李氏条件，212211),(),(yyLyxfyxf，则用(8.7)式减去(8.8)式，得到)1(11)1(1111)(112),(),(2ikkikkkkikkyyLhyxfyxfhy

8、y.从而递推得到)0(11)(11)2(kkiikkyyLhyy.如果选取h充分小，使2Lh，则当i时，就有迭代序列)(1iky收敛到1ky.【注】这里指出了迭代序列收敛到公式（8.7）的解，并不能由此断定迭代序列是否收敛到一阶初值问题（8.3）的解.7.2.4 改进的欧拉方法 梯形方法虽然提高了精度，但计算复杂，每迭代一次都要计算函数值，若干次的迭代使得计算量很大，而且很难预测.实际计算时，为简化计算，通常只迭代一两次就转入下一步的计算.具体地，先用欧拉公式求得预测值)0(1ky，这个预测值的精度较差，然后用梯形公式（8.7）将它校正一次得校正值1ky，这样建立的预测校正系统（Predict

9、or-Corrector Method）通常称为改进的欧拉方法(Modified Euler Method)，改进的欧拉公式为 1,1,0),(),(2),()0(111)0(1nkyxfyxfhyyyxhfyykkkkkkkkkk.（8.9）例 1 分别用 Euler 法，改进的 Euler 法解初值问题.1)0(1,0,3yxxyy 解 取5n，则5,4,3,2,1,0,2.0,2.051kkkhxhk.1）Euler 法计算公式为.4,1,0,1),3(2.001kyyxyykkkk.2）改进的 Euler 法计算公式为.4,1,0,1,)(3 1.0),3(2.00)0(111)0(1

10、kyyyxxyyyxyykkkkkkkkkk 计算结果见表 8-1.表 8-1 Euler 法与改进的 Euler 法计算结果对比 kx ky(Euler 法)ky(改进 Euler 法)(kxy 0 1 1 1 0.2 1.2 1.28 1.28561 0.4 1.5600 1.7536 1.76730 0.6 2.1120 2.463392 2.48848 0.8 2.8944 3.4613382 3.50216 1 3.9533 4.8108327 4.87313 利用 MATLAB 解析求解例 1 微分方程，输入 y=dsolve(Dy=y+3*x,y(0)=1,x)得到解析解 y=-

11、3*x-3+4*exp(x).由解析解xexy433计算出准确值)(kxy，与近似值ky一起列在表 8-1 中，两者比较可知，改进的 Euler 法比 Euler 法精度高.改进的 Euler 法具有二阶精度，这是因为),()(),(kkkkkyxfxyxyy,有)()()()(),(),(),(2211hOxyhxyhOyxf hyxfyxfkkkkkkkk,则改进的 Euler 法)()(2)(),(),(232111hOxyhxyhyyxfyxfhyykkkkkkkkk.比较 Taylor 展式 )(2)()()()(21kkkkkxyhxyhxyhxyxy，得到)()(311hOyxy

12、kk.Euler公式的MATLAB程序：Euler.m function x,y=Euler(fun,a,b,n,y0)%Euler公式，其中，%fun为一阶微分方程的函数；%a,b为求解区间的左右端点；%n为等分区间数；%y0为初始条件.x=zeros(1,n+1);y=zeros(1,n+1);h=(b-a)/n;x(1)=a;y(1)=y0;for k=1:n x(k+1)=x(k)+h;y(k+1)=y(k)+h*feval(fun,x(k),y(k);end 后退Euler公式的MATLAB程序：Euler_back.m function x,y=Euler_back(fun,a,b

13、,n,y0)%后退的Euler公式，其中，%fun为一阶微分方程的函数；%a,b为求解区间的左右端点；%n为等分区间数；%y0为初始条件.x=zeros(1,n+1);y=zeros(1,n+1);h=(b-a)/n;x(1)=a;y(1)=y0;for k=1:n x(k+1)=x(k)+h;%用迭代法求y(k+1)z0=y(k)+h*feval(fun,x(k),y(k);for i=1:5 z1=y(k)+h*feval(fun,x(k+1),z0);if abs(z1-z0)1e-3 break;end z0=z1;end y(k+1)=z1;end 梯形公式的MATLAB程序：Eul

14、er_trap.m function x,y=Euler_trap(fun,a,b,n,y0)%梯形公式，其中，%fun为一阶微分方程的函数；%a,b为求解区间的左右端点；%n为等分区间数；%y0为初始条件.x=zeros(1,n+1);y=zeros(1,n+1);h=(b-a)/n;x(1)=a;y(1)=y0;for k=1:n x(k+1)=x(k)+h;%用迭代法求y(k+1)z0=y(k)+h*feval(fun,x(k),y(k);for i=1:5 z1=y(k)+h/2*(feval(fun,x(k),y(k)+feval(fun,x(k+1),z0);if abs(z1-z

15、0)1e-3 break;end z0=z1;end y(k+1)=z1;end 改进的Euler公式MATLAB程序：Euler_correct.m function x,y=Euler_correct(fun,a,b,n,y0)%改进的Euler公式，其中，%fun为一阶微分方程的函数；%a,b为求解区间的左右端点；%n为等分区间数；%y0为初始条件.x=zeros(1,n+1);y=zeros(1,n+1);h=(b-a)/n;x(1)=a;y(1)=y0;for k=1:n x(k+1)=x(k)+h;y0=y(k)+h*feval(fun,x(k),y(k);y(k+1)=y(k)+

16、h/2*(feval(fun,x(k),y(k)+feval(fun,x(k+1),y0);end （1）对于例 1，利用 Euler 法函数 Euler.m 递推求解 fun=inline(y+3*x,x,y)a=0;b=1;n=5;y0=1;x,y=Euler(fun,a,b,n,y0)计算得 x=0 0.2000 0.4000 0.6000 0.8000 1.0000 y=1.0000 1.2000 1.5600 2.1120 2.8944 3.9533（2）对于例 1，利用后退 Euler 法函数 Euler_back.m 递推求解 fun=inline(y+3*x,x,y)a=0;b

17、=1;n=5;y0=1;x,y=Euler_back(fun,a,b,n,y0)计算得 x=0 0.2000 0.4000 0.6000 0.8000 1.0000 y=1.0000 1.3999 2.0498 3.0122 4.3651 6.2063（3）对于例 1，利用梯形法函数 Euler_trap.m 递推求解 fun=inline(y+3*x,x,y)a=0;b=1;n=5;y0=1;x,y=Euler_trap(fun,a,b,n,y0)计算得 x=0 0.2000 0.4000 0.6000 0.8000 1.0000 y=1.0000 1.2888 1.7751 2.5029

18、3.5257 4.9092（4）对于例 1，利用改进的 Euler 法函数 Euler_correct.m 递推求解 fun=inline(y+3*x,x,y)a=0;b=1;n=5;y0=1;x,y=Euler_correct(fun,a,b,n,y0)计算得 x=0 0.2000 0.4000 0.6000 0.8000 1.0000 y=1.0000 1.2800 1.7536 2.4634 3.4613 4.8108 【注】改进的 Euler 公式与梯形公式精度高于 Euler 公式与后退的 Euler 公式.前者采用两个节点上斜率（即节点函数值),(yxf）的平均值为斜率；而后者只用一个节点上的斜率.由此启发我们用多个节点斜率的加权平均为斜率，可能构造出精度更高的数值方法.