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

1、CENTRAL SOUTH UNIVERSITY数值分析实验报告Euler方法与改进的Euler方法的应用一、问题背景在工程和科学技术的实际问题中，常需求解微分方程，但常微分方程中往往只有少数较简单和典型的常微分方程（例如线性常系数常微分方程等）可求出其解析解,对于变系数常微分方程的解析求解就比较困难，而一般的非线性常微分方程的求解困难就更不用说了。大多数情况下，常微分方程只能用近似方法求解。这种近似解法可分为两大类：一类是近似解析法，如级数解法、逐次逼近法等;另一类是数值解法，它给出方程在一些离散点上的近似值。二、数学模型在具体求解微分方程时，需具备某种定解条件，微分方程和定解条件合在一起组

2、成定解问题。定解条件有两种：一种是给出积分曲线在初始点的状态，称为初始条件，相应的定解问题称为初值问题。另一类是给出积分曲线首尾两端的状态，称为边界条件，相应的定解问题称为边值问题。在本文中主要讨论的是给定初值条件的简单Euler方法和改进的Euler方法来求解常微分方程。三、算法及流程Euler方法是最简单的一种显式单步法。对于方程考虑用差商代替导数进行计算，取离散化点列，则得到方程的近似式即得到简单Euler方法。具体计算时由出发，根据初值，逐步递推二得到系列离散数值。简单Euler方法计算量小，然而精度却不高，因而我们可以构造梯形公式其中。这是一个二阶方法，比Euler方法精度高。但是上

3、述公式右边有，因而是隐式差分方程，可以用迭代方法计算。初值可以由Euler公式提供，一般而言迭代一两次即可。在计算中的迭代公式为容易看出这实际上是一种预估-校正方法。改进的Euler方法又称为Henu方法。MATLAB实现过程：（1）简单Euler方法函数functiont,x=Euler(fun,t0,tt,x0,N) % MyEuler用前向差分的欧拉方法解微分方程 % fun表示f（t,x） % t0,tt表示自变量的初值和终值 % x0表示函数在x0处的值，其可以为向量形式 % N表示自变量在t0,tt 上取得点数 h=(tt-t0)/N; % 步长h t=t0+0:N'*h;

4、 % 时间t x(1,:) = x0' % 赋初值 for k = 1:N f=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;%表示出离散的自变量x x(1,:)=x0' for i=1:N f1=h*feval(fun,t(i),x(i,:); f1=

5、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保存。四、计算结果与分析计算常微分方程的初值问题编写函数文件，打开Editor编辑器，输入以下语句并以文件名GjEuler.m保存。function f=GjEuler(t,x) f=1/t*(x2+x);再编写主函数文件，打开Editor编辑器，输入下列语句并以文件名GjEuler_main.m保存文件。function GjEuler_main() % 比较改进Euler法、简单Euler法及

6、文芬方程符号解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*(x2+x)','x(1)=-2','t'); % 符号计算for k=1:16 % 时间离散化 st(k)=t(k); % 取数值方法的离散的时间 sx(k)=subs(sh,st(k); end plot(t,x,'*',tgj,xg

7、j,'+',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, -38/2

8、37/5，-6694439885567781/4503599627370496,-43-14/923/15,-3-33468975972395-46/31669915380311/4503599627370496,-10/79/5,-5972253212408169/4503599627370496,-6244734289756621/4503599627370496,-18/1329/15,-14531/15,-285-297-62/4711/5,7/3,-5527499547085519/4503599627370496, -5738095816108559/450359962737049

9、6, -14/1137/15,-5455688972551323/4503599627370496,-5654358103897161/4503599627370496,-74/5913/5,41/15,-5338702448959659/4503599627370496,-275-82/6743/15,-529-272-86/713, -13521-6/5其中时间和改进的Euler方法两列的结果序列就是我们要的结果。为了分析问题，下图给出了分析解和Euler方法的结果进行比较。其中“*”线代表的是Euler法，“+”线代表的是改进Euler法，“”线代表的是分析解。可以看出，改进的Euler

10、方法比原方法要精确许多，在本例中的区间内，几乎改进后的方法和分析方法重合，而Euler方法则有比较大的误差。运行结果中还有下列数据：ans = 1, 0, 017/15, -16/285, 176/7267519/15, -6016/87975, 7/5, -2800435813076915/40532396646334464, 6797504630227/253327479039590423/15,-578292948083527/8725724278030336,5/3,-984409/5, -3425501531362731/58546795155816448, 29/15,-2651331/15, 11/5, -36899857630325880151/382805968326492167/3, -2247899765246235/49539595901075456, 68659194007205/4953959590107545637/15,-11380713/5, -3833349331259339/94575592174780416, 56376328559315/47287796087390208 3, -39224285523514680241/5629499534213120上面结果给出的是简单E