Euler方法与改进的Euler方法的应用

CENTRALSOUTHUNIVERSITY数值分析实验报告Euler方法与改进的Euler方法的应用一、问题背景在工程和科学技术的实际问题中，常需求解微分方程，但常微分方程中往往只有少数较简单和典型的常微分方程（例如线性常系数常微分方程等）可求出其解析解,对于变系数常微分方程的解析求解就比较困难，而一般的非线性常微分方程的求解困难就更不用说了。大多数情况下，常微分方程只能用近似方法求解。这种近似解法可分为两大类：一类是近似解析法，如级数解法、逐次逼近法等;另一类是数值解法，它给出方程在一些离散点上的近似值。二、数学模型在具体求解微分方程时，需具备某种定解条件，微分方程和定解条件合在一起组成定解问题。定解条件有两种：一种是给出积分曲线在初始点的状态，称为初始条件，相应的定解问题称为初值问题。另一类是给出积分曲线首尾两端的状态，称为边界条件，相应的定解问题称为边值问题。在本文中主要讨论的是给定初值条件的简单Euler方法和改进的Euler方法来求解常微分方程。三、算法及流程Euler方法是最简单的一种显式单步法。对于方程yxfdxdy,考虑用差商代替导数进行计算，取离散化点列nhxxn0，Ln,2,1,0则得到方程的近似式nnnnxyxfhxyxy,1即nnnnyxhfyy,1得到简单Euler方法。具体计算时由0x出发，根据初值，逐步递推二得到系列离散数值。简单Euler方法计算量小，然而精度却不高，因而我们可以构造梯形公式0111,,2yytfytfhyynnnnnn其中Nabh。这是一个二阶方法，比Euler方法精度高。但是上述公式右边有1ny，因而是隐式差分方程，可以用迭代方法计算1ny。初值可以由Euler公式提供，一般而言迭代一两次即可。在计算中的迭代公式为Lkyxfyxfhyyyxhfyynknnnnnknnnn,2,1,0,,2,111110容易看出这实际上是一种预估-校正方法。改进的Euler方法又称为Henu方法。MATLAB实现过程：（1）简单Euler方法函数function[t,x]=Euler(fun,t0,tt,x0,N)%MyEuler用前向差分的欧拉方法解微分方程%fun表示f（t,x）%t0,tt表示自变量的初值和终值%x0表示函数在x0处的值，其可以为向量形式%N表示自变量在[t0,tt]上取得点数h=(tt-t0)/N;%步长ht=t0+[0:N]'*h;%时间tx(1,:)=x0';%赋初值fork=1:Nf=feval(fun,t(k),x(k,:));f=f';x(k+1,:)=x(k,:)+h*f;end将文件以文件名Euler.m保存。（2）改进的Euler方法函数function[t,x]=GjEuler(fun,t0,tt,x0,N)%改进的Euler方法解微分方程h=(tt-t0)/N;%计算所取的两离散点之间的距离t=t0+[0:N]'*h;%表示出离散的自变量xx(1,:)=x0';fori=1:Nf1=h*feval(fun,t(i),x(i,:));f1=f1';f2=h*feval(fun,t(i+1),x(i,:)+f1);f2=f2';x(i+1,:)=x(i,:)+1/2*(f1+f2);end将文件以文件名GjEuler.m保存。四、计算结果与分析计算常微分方程的初值问题21,3,112xtxxtdtdx编写函数文件，打开Editor编辑器，输入以下语句并以文件名GjEuler.m保存。functionf=GjEuler(t,x)f=1/t*(x^2+x);再编写主函数文件，打开Editor编辑器，输入下列语句并以文件名GjEuler_main.m保存文件。functionGjEuler_main()%比较改进Euler法、简单Euler法及文芬方程符号解[t,x]=Euler('GjEuler_fun',1,3,-2,15);%Euler方法¨[tgj,xgj]=GjEuler('GjEuler_fun',1,3,-2,15);%改进的Euler方法的解sh=dsolve('Dx=1/t*(x^2+x)','x(1)=-2','t');%符号计算fork=1:16%时间离散化st(k)=t(k);%取数值方法的离散的时间sx(k)=subs(sh,st(k));endplot(t,x,'*',tgj,xgj,'+',st,sx)[t,x,xgj,sx']%t时间%x简单Euler方法计算结果%xgj改进Euler方法的计算结果%sx符号计算结果[t,sx'-x,sx'-xgj]%两种方法的计算误差%plot(t,sx'-x,'+',t,sx'-xgj,'*')运行程序，在MATLAB命令窗口输入：GjEuler_main运行结果如下所示：ans=[1,-2,-2,-2][17/15,-26/15,-6854/3825,-34/19][19/15,-6058/3825,-3726675967822917/2251799813685248,-38/23][7/5，-6694439885567781/4503599627370496,-438605242064379/281474976710656,-14/9][23/15,-399017934858279/281474976710656,-3346897597239529/2251799813685248,-46/31][5/3,-3076228012319145/2251799813685248,-6443669915380311/4503599627370496,-10/7][9/5,-5972253212408169/4503599627370496,-6244734289756621/4503599627370496,-18/13][29/15,-1456996799034161/1125899906842624,-3041377995996661/2251799813685248,-58/43][31/15,-2854895172348175/2251799813685248,-2974161832326041/2251799813685248,-62/47][11/5,-5611129535974289/4503599627370496,-2917483895703271/2251799813685248,-22/17][7/3,-5527499547085519/4503599627370496,-5738095816108559/4503599627370496,-14/11][37/15,-5455688972551323/4503599627370496,-5654358103897161/4503599627370496,-74/59][13/5,-5393344808589217/4503599627370496,-2790627213541703/2251799813685248,-26/21][41/15,-5338702448959659/4503599627370496,-2758440471997471/2251799813685248,-82/67][43/15,-5290411912308893/4503599627370496,-2729881239997317/2251799813685248,-86/71][3,-2623711206582523/2251799813685248,-1352184591147197/1125899906842624,-6/5]其中时间和改进的Euler方法两列的结果序列就是我们要的结果。为了分析问题，下图给出了分析解和Euler方法的结果进行比较。其中“*”线代表的是Euler法，“+”线代表的是改进Euler法，“—”线代表的是分析解。可以看出，改进的Euler方法比原方法要精确许多，在本例中的区间内，几乎改进后的方法和分析方法重合，而Euler方法则有比较大的误差。运行结果中还有下列数据：ans=[1,0,0][17/15,-16/285,176/72675][19/15,-6016/87975,45154339887667/51791395714760704][7/5,-2800435813076915/40532396646334464,6797504630227/2533274790395904][23/15,-578292948083527/8725724278030336,171034084903991/69805794224242688][5/3,-984402050618465/15762598695796736,69693133957217/31525197391593472][9/5,-3425501531362731/58546795155816448,116752474167145/58546795155816448][29/15,-2651332238403269/48413695994232832,174864634112039/96827391988465664][31/15,-5431515348121151/105834591243206656,174017670838551/105834591243206656][11/5,-3689989690587999/76561193665298432,57630325880151/38280596832649216][7/3,-2247899765246235/49539595901075456,68659194007205/49539595901075456][37/15,-11380723044888647/265712378014859264,340755704515795/265712378014859264][13/5,-3833349331259339/94575592174780416,56376328559315/47287796087390208][41/15,-11602105364083519/301741175033823232,167926901640221/150870587516911616][43/15,-11690322179931253/319755573543305216,166784062878179/159877786771652608][3,-392242849198873/11258999068426240,5523514680241/5629499534213120]上面结果给出的是简单Euler、改进Euler与分析解的误差比较结果，将其画出图像得到下面的图形：“+”线表示的是Euler方法误差，“*”线表示的是改进方法的误差。当将步长减小时，可以得到精度比原来更高，误差也在减小，下图给出的是将步长减小到35步时的对比图。可以明显看到，在同样区间上，Euler法和改进Euler法都比原来更为接近分析解，足以说明精度在提高。