《《计算方法第八章》PPT课件》由会员分享，可在线阅读，更多相关《《计算方法第八章》PPT课件（58页珍藏版）》请在金锄头文库上搜索。
1、第八章第八章 常微分方程初值常微分方程初值问题数值解法问题数值解法本章主要研究常微分方程初值问题的数值求解：本章主要研究常微分方程初值问题的数值求解：通常，假设函数通常，假设函数 f 关于第二个变量满足李普希茨条件（关于第二个变量满足李普希茨条件（L条条件），即为存在常数件），即为存在常数 L 0，使得使得第一节第一节 一般概念一般概念1.1 欧拉法及其简单改进欧拉法及其简单改进方法：选择适当的节点，用差分近似微分，将方程离散化，方法：选择适当的节点，用差分近似微分，将方程离散化，从而求在这些节点上的函数值。从而求在这些节点上的函数值。欧拉欧拉方法方法例子：例子：下面我们分别取步长为下面我们分
2、别取步长为0.1与与0.01进行计算，进行计算，计算结果显示在下面的图中。计算结果显示在下面的图中。步长为步长为0.1的计算结果。的计算结果。步长为步长为0.01的计算结果的计算结果0.01 0.99005 0.99 0.1 0.90484 0.90438 0.2 0.81873 0.81791 0.3 0.74082 0.7397 0.4 0.67032 0.66897 0.41 0.66365 0.66228 0.59 0.55433 0.55268 0.6 0.54881 0.54716 0.9 0.40657 0.40473 0.91 0.40252 0.40068 0.99 0.37
3、158 0.36973 1 0.36788 0.36603 DOUBLE PRECISION h,y(0：100) OPEN(20,FILE=OUTPUT.DAT,STATUS=UNKNOWN)h=1.0/100y(0)=1.0do 10 i=1,100 y(i)=y(i-1)*(1.0-h)write(20,*) i*h,y(i)10 continueEND精度更高的欧拉公式：精度更高的欧拉公式：方法：选择计算导数值方法：选择计算导数值的精度更好的差分格式。的精度更好的差分格式。欧拉中点公式欧拉中点公式利用中点公式求解微分方程时，有一个问题，就是计算时需利用中点公式求解微分方程时，有一个问题
4、，就是计算时需要两个迭代初值！要两个迭代初值！对于这个问题，我们可以先用欧拉公式，通过给定的初值计对于这个问题，我们可以先用欧拉公式，通过给定的初值计算出第一个点的值，然后在利用这两个（第一和第二个点）算出第一个点的值，然后在利用这两个（第一和第二个点）的值进行计算，直到计算出全部节点上的值。的值进行计算，直到计算出全部节点上的值。下面，我们用中点公式求解刚才的例子，计算的步长取下面，我们用中点公式求解刚才的例子，计算的步长取0.01，可以看到，计算的精度比较高。，可以看到，计算的精度比较高。此时，计算公式为此时，计算公式为计计算算结结果果0.02 0.9802 0.98020.1 0.904
5、84 0.904840.2 0.81873 0.818740.3 0.74082 0.740840.4 0.67032 0.670350.5 0.60653 0.606560.6 0.54881 0.548850.7 0.49659 0.496630.8 0.44933 0.449380.9 0.40657 0.406631.2 欧拉方法的其他改进欧拉方法的其他改进微分方程数值解的关键在于对导数的处理，可以用差分来近似微分方程数值解的关键在于对导数的处理，可以用差分来近似导数，也可以通过积分，将导数项化掉。导数，也可以通过积分，将导数项化掉。对于方程：对于方程：首先，作出划分首先，作出划分设已
6、经求出第设已经求出第 n 个个节点的函数值节点的函数值 ，在区间，在区间 上对上对方程两边积分方程两边积分容易看出，要求第容易看出，要求第 n+1 个个节点的函数值，关键在于节点的函数值，关键在于选择适当的选择适当的积分公式计算积分积分公式计算积分！（1）如选择下矩形公式，则）如选择下矩形公式，则这正是前面的这正是前面的欧拉公式欧拉公式。（2）如选择上矩形公式，则）如选择上矩形公式，则这是所谓的这是所谓的后退欧拉公式后退欧拉公式。（3）如选择梯形公式，则）如选择梯形公式，则这是所谓的这是所谓的欧拉梯形公式欧拉梯形公式。直接利用已经求得的已知节点上的值计算未知节点上的函数直接利用已经求得的已知节
7、点上的值计算未知节点上的函数值的算法称为值的算法称为显式法显式法。例如：欧拉公式、欧拉中点公式例如：欧拉公式、欧拉中点公式计算未知节点上的函数值时，用到了未知节点上的函数值，计算未知节点上的函数值时，用到了未知节点上的函数值，这种算法称为这种算法称为隐式法。隐式法。例如：后退欧拉法、欧拉梯形公式例如：后退欧拉法、欧拉梯形公式显然，利用隐式法求微分方程的数值解是，显然，利用隐式法求微分方程的数值解是，需要从表达式中需要从表达式中反解未知节点上的函数值反解未知节点上的函数值。1.3 隐式法的具体计算：隐式法的具体计算：例如欧拉梯形公式例如欧拉梯形公式用迭代法计算用迭代法计算n+1步的值。步的值。（
8、1）简单迭代）简单迭代收敛的条件：收敛的条件：（2）牛顿迭代）牛顿迭代必须指出，在真正计算中，常用的是简单迭代，而且只迭代一必须指出，在真正计算中，常用的是简单迭代，而且只迭代一步，迭代初值步，迭代初值 称为预测，迭代称为校正，这样组成的称为预测，迭代称为校正，这样组成的一组计算公式称为一组计算公式称为预测校正公式预测校正公式。预测校正公式也称为改进的欧拉法，将上面的组合公式改预测校正公式也称为改进的欧拉法，将上面的组合公式改写为：写为：注意到注意到 ，将上式进一步改写为：，将上式进一步改写为：这是我们最终使用的计算格式。这是我们最终使用的计算格式。例子：例子：取步长为取步长为0.1计算，结果
9、如图。计算，结果如图。图：图： DOUBLE PRECISION h,y(0:10),ak1,ak2 OPEN(20,FILE=OUTPUT1.DAT,STATUS=UNKNOWN)h=1.0/10y(0)=1.0do 10 i=1,10 ak1=-h*y(i-1) ak2=-h*(y(i-1)+ak1) y(i)=y(i-1)+(ak1+ak2)/2.010 continue do 20 i=0,10 write(20,*) i*h,y(i),exp(-i*h)20 continueEND同理，对于后退欧拉公式同理，对于后退欧拉公式有预测校正公式有预测校正公式或改写为：或改写为：用此法解前面
10、的例子用此法解前面的例子步长步长0.1步长步长0.011.4 误差估计误差估计定义：定义：利用第利用第n个节点的函数值，通过数值方法计算第个节点的函数值，通过数值方法计算第n+1个节个节点的函数近似值，所引起的误差，称为点的函数近似值，所引起的误差，称为n+1个节点上的个节点上的局部局部截断误差截断误差。我们记我们记 为为n+1个节点上函数的精确值，个节点上函数的精确值， 为数值解，为数值解，则局部截断误差为：则局部截断误差为：如局部截断误差为如局部截断误差为 ，称为具有，称为具有p阶局部截断误差。阶局部截断误差。欧拉方法的误差分析：欧拉方法的误差分析：省略余项得公式：省略余项得公式：完全类似
11、的可以得到完全类似的可以得到后退欧拉公式的局部截断误差为：后退欧拉公式的局部截断误差为：欧拉中点公式的局部截断误差为：欧拉中点公式的局部截断误差为：欧拉梯形公式的局部截断误差为：欧拉梯形公式的局部截断误差为：定义：由初值引起的定义：由初值引起的第第 n 个节点的近似值与精确值之间的误差个节点的近似值与精确值之间的误差称为第称为第 n 个节点整体误差。个节点整体误差。定理：设下面求解微分方程的数值计算方法定理：设下面求解微分方程的数值计算方法局部截断误差为局部截断误差为p+1阶，且函数阶，且函数 关于关于 y 满足利普希茨条件，满足利普希茨条件，同时初值是准确的，则同时初值是准确的，则整体截断误
12、差为整体截断误差为p阶阶。欧拉公式、后退欧拉公式的整体截断误差为欧拉公式、后退欧拉公式的整体截断误差为 1 阶。阶。欧拉中点公式、欧拉梯形公式的整体截断误差为欧拉中点公式、欧拉梯形公式的整体截断误差为 2 阶。阶。微分方程数值解法的进一步改进。再回到恒等式微分方程数值解法的进一步改进。再回到恒等式如果取如果取 作为节点，将被积函数用插值多作为节点，将被积函数用插值多项式来近似，用插值多项式带到积分中去求出积分，则可以得项式来近似，用插值多项式带到积分中去求出积分，则可以得到所谓的到所谓的亚当斯亚当斯(Adams)显式公式显式公式局部截断误差：局部截断误差：类似地，如果取类似地，如果取 作为作为
13、节点，可得节点，可得亚当斯亚当斯(Adams)隐式公式隐式公式局部截断误差：局部截断误差：进一步，如果我们将恒等式中的积分区间改为进一步，如果我们将恒等式中的积分区间改为 ，并在此区间上用辛普森公式，可得，并在此区间上用辛普森公式，可得1.5 绝对稳定性绝对稳定性一个常微分方程数值解法称为绝对稳定，是指在某一步（如第一个常微分方程数值解法称为绝对稳定，是指在某一步（如第一步）产生的误差（如计算机的存储误差），在计算中会逐一步）产生的误差（如计算机的存储误差），在计算中会逐步减小。步减小。称方程称方程 为为试验方程试验方程，设计算，设计算步长为步长为 h ，设在计算开设在计算开始时产生误差（存储
14、误差），此误差在以后会逐步减弱，我始时产生误差（存储误差），此误差在以后会逐步减弱，我们称该算法相对于们称该算法相对于 是绝对稳定的，是绝对稳定的， 的全体称的全体称为该算法的绝对稳定域。为该算法的绝对稳定域。 欧拉法的绝对稳定区域欧拉法的绝对稳定区域 后退欧拉法的绝对稳定区域后退欧拉法的绝对稳定区域1.6 局部截断误差的实用估计局部截断误差的实用估计（1）用两种阶数相同的算法求解，计算出）用两种阶数相同的算法求解，计算出n+1步的近似值，步的近似值，从而得到局部截断误差估计。从而得到局部截断误差估计。（2）用同样的公式，用不同步长计算出）用同样的公式，用不同步长计算出n+1步的近似值，从步的
15、近似值，从而得到局部截断误差估计。而得到局部截断误差估计。1.7 隐式法隐式法隐式法具有较好的绝对稳定性隐式法具有较好的绝对稳定性！只不过在使用隐式法的时候，需要进行迭代，或者使用预测只不过在使用隐式法的时候，需要进行迭代，或者使用预测校正计算格式。校正计算格式。第二节第二节 泰勒级数法与龙格库塔法泰勒级数法与龙格库塔法对于方程：对于方程：取计算步长取计算步长为为 h ，则则 ，将函数进行泰勒展开，将函数进行泰勒展开如函数如函数 y( x )有有p+1阶导数，容易得到阶导数，容易得到p阶泰勒级数展开法：阶泰勒级数展开法：公式中的导数用下面公式计算：公式中的导数用下面公式计算：例子：例子：步长步
16、长0.1 步长步长0.01龙格库塔法：龙格库塔法：对于常微分方程的数值解法，一个关键在于选择精度高的算法对于常微分方程的数值解法，一个关键在于选择精度高的算法计算下面公式中的积分。计算下面公式中的积分。要高精度的计算积分，常用的方法是适当增加计算节点，不妨要高精度的计算积分，常用的方法是适当增加计算节点，不妨设用设用 m 个节点计算上面积分，节点为个节点计算上面积分，节点为则积分为则积分为将积分改写为：将积分改写为：则得公式：则得公式：这样的公式称为这样的公式称为显式龙格库塔显式龙格库塔公式。公式。确定二确定二阶阶 RK 法法:系数满足：系数满足：此为四个未知数的三个方程，任意取此为四个未知数
17、的三个方程，任意取 ，得，得特别，取特别，取 ，得到，得到通常也称为变形欧拉法，也常写为通常也称为变形欧拉法，也常写为它具有二阶精度，也称为它具有二阶精度，也称为二阶二二阶二级级 RK 方法方法。例子例子三阶三级库塔法三阶三级库塔法局部截断误差为局部截断误差为4阶，整体截断误差为阶，整体截断误差为3阶阶最常用的四级四阶公式，称为最常用的四级四阶公式，称为龙格库塔公式龙格库塔公式：局部截断误差为局部截断误差为5阶，整体截断误差为阶，整体截断误差为4阶。阶。隐式龙格库塔法：隐式龙格库塔法：常用的有二阶常用的有二阶RK法：法：例子：例子：用隐式用隐式二阶二阶RK法：法：用显式用显式二阶二阶RK法：法
18、： DOUBLE PRECISION h,y(0:100),yy(0:100)double precision ak OPEN(20,FILE=OUTPUT5.DAT,STATUS=UNKNOWN)h=1.0/100y(0)=1.0yy(0)=1.0do 10 i=1,100 y(i)=y(i-1)-(h/(1-h/2)*y(i-1) yy(i)=yy(i-1)*(1-h+h*h/2)10 continue do 20 i=0,100 write(20,100) i*h,y(i),yy(i),exp(-i*h)100 format (1x,f4.2,f8.5,f8.5,f8.5)20 cont
19、inueENDh=0.1h=0.01半隐式龙格库塔法：半隐式龙格库塔法：最常用的二级三阶半隐式龙格库塔公式：最常用的二级三阶半隐式龙格库塔公式：微分方程组微分方程组一阶方程组一阶方程组龙格库塔公式龙格库塔公式写成分量形式：写成分量形式：例子：例子： DOUBLE PRECISION h,y1(0:100),y2(0:100)double precision ak1,ak2,ak3,ak4,bk1,bk2,bk3,bk4 OPEN(20,FILE=OUTPUT6.DAT,STATUS=UNKNOWN)h=1.0/10y1(0)=0.0y2(0)=1.0do 10 i=1,10 ak1=h*(-y
20、1(i-1)+y2(i-1) bk1=-h*y2(i-1) ak2=h*(-(y1(i-1)+ak1/2.0)+y2(i-1)+bk1/2.0) bk2=-h*(y2(i-1)+bk1/2.0) ak3=h*(-(y1(i-1)+ak2/2.0)+y2(i-1)+bk2/2.0) bk3=-h*(y2(i-1)+bk2/2.0) ak4=h*(-(y1(i-1)+ak3)+y2(i-1)+bk3) bk4=-h*(y2(i-1)+bk3) y1(i)=y1(i-1)+(ak1+2*ak2+2*ak3+ak4)/6.0 y2(i)=y2(i-1)+(bk1+2*bk2+2*bk3+bk4)/6.
21、010 continue do 20 i=0,10 write(20,100) i*h,y1(i),i*h*exp(-i*h),y2(i),exp(-i*h)100 format (1x,f4.2,f12.5,f12.5,f12.5,f12.5)20 continueEND取取h=0.1节点节点 y1y1精确精确 y1y1数值数值 y2y2精确精确 y2y2数值数值 .00 .00000 .00000 1.00000 1.00000 .10 .09048 .09048 .90484 .90484 .20 .16375 .16375 .81873 .81873 .30 .22224 .22225 .74082 .74082 .40 .26813 .26813 .67032 .67032 .50 .30326 .30327 .60653 .60653 .60 .32929 .32929 .54881 .54881 .70 .34761 .34761 .49659 .49659 .80 .35946 .35946 .44933 .44933 .90 .36591 .36591 .40657 .40657 1.00 .36788 .36788 .36788 .36788高阶微分方程可以转换为一阶方程组进行求解。高阶微分方程可以转换为一阶方程组进行求解。令：令：
