科学计算选讲.doc

科学计算选讲 两种EULER法求解常微分方程初值的比较 学院：精密仪器与光电子工程学院 专业：仪器科学与技术 姓名：崔建军 学号：1011202025 一、微分方程的初值计算方法及其原理分析 对于微分方程的解析表达式 （1） 则它的初值表达式： （2） 依据：定理 （存在与唯一性定理）如果方程（1）的右端函数在闭矩形域 上满足如下条件： （1）在上连续； （2）在上关于变量满足利普希茨（Lipschitz）条件，即存在常数，使 对于上任何一对点和有不等式： , 则初值问题（2）在区间上存在唯一解 , 其中． 根据存在与唯一性定理，只要关于满足Lipschitz条件 ， 即可保证其解存在并唯一。 不过解析法只能求解极少数典型的常微分方程，如线性常系数微分方程等，对于变系数常微分方程非线性常微分方程解析法计算变得十分困难，实际应用中建立的微分方程物理模型主要采用数值分析法求解。目前采用具体的数值算法也比较多，下面主要介绍改进的EULER法EULER算法 首先，基本思想是迭代。EULER法取切线的端点作为下一步的起点进行计算，代替y(x)，即在节点处用差商近似代替导数，当f(x,y(x))中的x取[xn,xn+1]的左端点xn，即f(xn,y(xn))。通过泰勒级数展开方法，能够获得证明。具体过程如下： 在处泰勒展开有 , （4） 略去余项，得 . 用近似值代替，把上式右端所得结果记为，得 . 或者也可在区间上对微分方程（1）进行积分，得 , （5） 利用左矩形公式得 . 用代替，把上式右端所得结果记为，即得公式（3）. 2.1 前向EULER法 如果将y(xn)的近似值记为yn，即yn≈y(xn)，yn+1≈y(xn+1)，则得到前向EULER公式：。 2.2 后向EULER法 当f(x,y(x))中的x取[xn,xn+1]的右端点xn+1，即f(xn+1,y(xn+1))，可得到后向欧拉公式：。 2.3 改进的EULER法 如果用梯形公式计算（5）式右端的积分，得 . 对于上式右端，用代替，把上式右端所得结果记为，即得梯形公式 . （6） 梯形公式（6）是关于的隐式方程，因此称梯形法为隐式方法，而公式（3）是的显式形式，因此称欧拉法为显式方法。 或者先用Euler格式求得一个初步的近似值，称之为预测值，预测值的精度可能很差，再用梯形公式将它校正一次，得，这个结果称之为校正值．而这样建立的预测—校正系统通常称为改进的Euler格式． 预测 校正 这一计算格式亦可表示为： , 三、截断误差与精度分析 在实际计算中，即使初始值是准确的，但所得的数值解往往与初始值问题的准确解有一定的误差。设 , 表示在节点处准确解与数值解之差，称之为整体截断误差，它不仅与这步计算有关，而且与这些步计算有关。假定为准确的，即在的前提下估计误差 , 这种误差称为局部截断误差。如果局部截断误差这里是等距节点的步长，则称所用的数值解的方法的阶数是，或称该方法有阶精确度.显然，步长越小，阶数越高，则局部截断误差越小，计算结果的精度也越高。 对于EULER法，通过逐次替代，最后求出所要求的解，并达到一定的精度。单步，显式，一阶求导精度，截断误差阶当步数增多时，误差会因积累而越来越大。采用区间两端的函数值的平均值作为直线方程的斜率。改进欧拉法的精度为二阶。 4.1 微分方程的函数： function z=myf4(t,u) z=u-2*t/u; 4.2 准确解的函数文件： function fz=fz(x) fz=(2*t + 1)^(1/2); 通过如下语句得到上式。 dsolve(Du = u-2*t/u;,u(0) = 1,s) 4.3 欧拉格式： x=zeros(1,11); y=zeros(1,11); x(1)=0; y(1)=0; h=0.1; for n=1:10 x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(@myf4,x(n),y(n)); end E=[x,y] y1=fz(x) plot(x,y,:or,x,y1,--*b); title(Euler格式与准确解比较图); 4.4 改进的Euler梯形格式： x=zeros(1,11); y=zeros(1,11); x(1