计算方法第八章优秀PPT.ppt

计算方法第八章现在学习的是第1页，共58页本章主要研究常微分方程初值问题的数值求解：本章主要研究常微分方程初值问题的数值求解：通常，假设函数通常，假设函数 f 关于第二个变量满足李普希茨条件（关于第二个变量满足李普希茨条件（L条条件），即为存在常数件），即为存在常数 L 0，使得使得现在学习的是第2页，共58页第一节第一节 一般概念一般概念1.1 欧拉法及其简单改进欧拉法及其简单改进方法：选择适当的节点，用差分近似微分，将方程离散化，方法：选择适当的节点，用差分近似微分，将方程离散化，从而求在这些节点上的函数值。从而求在这些节点上的函数值。欧拉欧拉方法方法现在学习的是第3页，共58页例子：例子：下面我们分别取步长为下面我们分别取步长为0.1与与0.01进行计算，进行计算，计算结果显示在下面的图中。计算结果显示在下面的图中。现在学习的是第4页，共58页步长为步长为0.1的计算结果。的计算结果。现在学习的是第5页，共58页步长为步长为0.01的计算结果的计算结果现在学习的是第6页，共58页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.37158 0.36973 1 0.36788 0.36603 现在学习的是第7页，共58页 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现在学习的是第8页，共58页精度更高的欧拉公式：精度更高的欧拉公式：方法：选择计算导数值方法：选择计算导数值的精度更好的差分格式。的精度更好的差分格式。欧拉中点公式欧拉中点公式现在学习的是第9页，共58页利用中点公式求解微分方程时，有一个问题，就是计算时需要两利用中点公式求解微分方程时，有一个问题，就是计算时需要两个迭代初值！个迭代初值！对于这个问题，我们可以先用欧拉公式，通过给定的初值计算出第一对于这个问题，我们可以先用欧拉公式，通过给定的初值计算出第一个点的值，然后在利用这两个（第一和第二个点）的值进行计算，个点的值，然后在利用这两个（第一和第二个点）的值进行计算，直到计算出全部节点上的值。直到计算出全部节点上的值。下面，我们用中点公式求解刚才的例子，计算的步长取下面，我们用中点公式求解刚才的例子，计算的步长取0.01，可以看，可以看到，计算的精度比较高。到，计算的精度比较高。此时，计算公式为此时，计算公式为现在学习的是第10页，共58页计计算算结结果果现在学习的是第11页，共58页0.02 0.9802 0.98020.1 0.90484 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.40663现在学习的是第12页，共58页1.2 欧拉方法的其他改进欧拉方法的其他改进微分方程数值解的关键在于对导数的处理，可以用差分来近似微分方程数值解的关键在于对导数的处理，可以用差分来近似导数，也可以通过积分，将导数项化掉。导数，也可以通过积分，将导数项化掉。对于方程：对于方程：首先，作出划分首先，作出划分设已经求出第设已经求出第 n 个个节点的函数值节点的函数值 ，在区间，在区间 上对方程两边积上对方程两边积分分容易看出，要求第容易看出，要求第 n+1 个个节点的函数值，关键在于节点的函数值，关键在于选择适当的积分公式选择适当的积分公式计算积分计算积分！现在学习的是第13页，共58页（1）如选择下矩形公式，则）如选择下矩形公式，则这正是前面的这正是前面的欧拉公式欧拉公式。（2）如选择上矩形公式，则）如选择上矩形公式，则这是所谓的这是所谓的后退欧拉公式后退欧拉公式。（3）如选择梯形公式，则）如选择梯形公式，则这是所谓的这是所谓的欧拉梯形公式欧拉梯形公式。现在学习的是第14页，共58页直接利用已经求得的已知节点上的值计算未知节点上的函数值的算法称为直接利用已经求得的已知节点上的值计算未知节点上的函数值的算法称为显式法显式法。例如：欧拉公式、欧拉中点公式例如：欧拉公式、欧拉中点公式计算未知节点上的函数值时，用到了未知节点上的函数值，这种算法称为计算未知节点上的函数值时，用到了未知节点上的函数值，这种算法称为隐式法。隐式法。例如：后退欧拉法、欧拉梯形公式例如：后退欧拉法、欧拉梯形公式显然，利用隐式法求微分方程的数值解是，显然，利用隐式法求微分方程的数值解是，需要从表达式中反解未知需要从表达式中反解未知节点上的函数值节点上的函数值。现在学习的是第15页，共58页1.3 隐式法的具体计算：隐式法的具体计算：例如欧拉梯形公式例如欧拉梯形公式用迭代法计算用迭代法计算n+1步的值。步的值。（1）简单迭代）简单迭代收敛的条件：收敛的条件：现在学习的是第16页，共58页（2）牛顿迭代）牛顿迭代必须指出，在真正计算中，常用的是简单迭代，而且只迭代一步，必须指出，在真正计算中，常用的是简单迭代，而且只迭代一步，迭代初值迭代初值 称为预测，迭代称为校正，这样组成的一组计称为预测，迭代称为校正，这样组成的一组计算公式称为算公式称为预测校正公式预测校正公式。现在学习的是第17页，共58页预测校正公式也称为改进的欧拉法，将上面的组合公式改写为：预测校正公式也称为改进的欧拉法，将上面的组合公式改写为：注意到注意到 ，将上式进一步改写为：，将上式进一步改写为：这是我们最终使用的计算格式。这是我们最终使用的计算格式。现在学习的是第18页，共58页例子：例子：取步长为取步长为0.1计算，结果如图。计算，结果如图。现在学习的是第19页，共58页图：图：现在学习的是第20页，共58页 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现在学习的是第21页，共58页同理，对于后退欧拉公式同理，对于后退欧拉公式有预测校正公式有预测校正公式或改写为：或改写为：现在学习的是第22页，共58页用此法解前面的例子用此法解前面的例子步长步长0.1步长步长0.01现在学习的是第23页，共58页1.4 误差估计误差估计定义：定义：利用第利用第n个节点的函数值，通过数值方法计算第个节点的函数值，通过数值方法计算第n+1个节点的个节点的函数近似值，所引起的误差，称为函数近似值，所引起的误差，称为n+1个节点上的个节点上的局部截断误差局部截断误差。我们记我们记 为为n+1个节点上函数的精确值，个节点上函数的精确值，为数值解，为数值解，则局部截断误差为：则局部截断误差为：如局部截断误差为如局部截断误差为 ，称为具有，称为具有p阶局部截断误差。阶局部截断误差。现在学习的是第24页，共58页欧拉方法的误差分析：欧拉方法的误差分析：省略余项得公式：省略余项得公式：现在学习的是第25页，共58页完全类似的可以得到完全类似的可以得到后退欧拉公式的局部截断误差为：后退欧拉公式的局部截断误差为：欧拉中点公式的局部截断误差为：欧拉中点公式的局部截断误差为：欧拉梯形公式的局部截断误差为：欧拉梯形公式的局部截断误差为：现在学习的是第26页，共58页定义：由初值引起的定义：由初值引起的第第 n 个节点的近似值与精确值之间的误差个节点的近似值与精确值之间的误差称为第称为第 n 个节点整体误差。个节点整体误差。定理：设下面求解微分方程的数值计算方法定理：设下面求解微分方程的数值计算方法局部截断误差为局部截断误差为p+1阶，且函数阶，且函数 关于关于 y 满足利普希茨条件，满足利普希茨条件，同时初值是准确的，则同时初值是准确的，则整体截断误差为整体截断误差为p阶阶。欧拉公式、后退欧拉公式的整体截断误差为欧拉公式、后退欧拉公式的整体截断误差为 1 阶。阶。欧拉中点公式、欧拉梯形公式的整体截断误差为欧拉中点公式、欧拉梯形公式的整体截断误差为 2 阶。阶。现在学习的是第27页，共58页微分方程数值解法的进一步改进。再回到恒等式微分方程数值解法的进一步改进。再回到恒等式如果取如果取 作为节点，将被积函数用插值多作为节点，将被积函数用插值多项式来近似，用插值多项式带到积分中去求出积分，则可以得项式来近似，用插值多项式带到积分中去求出积分，则可以得到所谓的到所谓的亚当斯亚当斯(Adams)显式公式显式公式局部截断误差：局部截断误差：现在学习的是第28页，共58页类似地，如果取类似地，如果取 作为作为节点，可得节点，可得亚当斯亚当斯(Adams)隐式公式隐式公式局部截断误差：局部截断误差：现在学习的是第29页，共58页进一步，如果我们将恒等式中的积分区间改为进一步，如果我们将恒等式中的积分区间改为 ，并在此区间上用辛普森公式，可得，并在此区间上用辛普森公式，可得现在学习的是第30页，共58页1.5 绝对稳定性绝对稳定性一个常微分方程数值解法称为绝对稳定，是指在某一步（如第一步）产一个常微分方程数值解法称为绝对稳定，是指在某一步（如第一步）产生的误差（如计算机的存储误差），在计算中会逐步减小。生的误差（如计算机的存储误差），在计算中会逐步减小。称方程称方程 为为试验方程试验方程，设计算，设计算步长为步长为 h，设在计算开始时产设在计算开始时产生误差（存储误差），此误差在以后会逐步减弱，我们称该算法生误差（存储误差），此误差在以后会逐步减弱，我们称该算法相对于相对于 是绝对稳定的，是绝对稳定的，的全体称为该算法的绝对的全体称为该算法的绝对稳定域。稳定域。现在学习的是第31页，共58页欧拉法的绝对稳定区域欧拉法的绝对稳定区域 后退欧拉法的绝对稳定区域后退欧拉法的绝对稳定区域现在学习的是第32页，共58页1.6 局部截断误差的实用估计局部截断误差的实用估计（1）用两种阶数相同的算法求解，计算出）用两种阶数相同的算法求解，计算出n+1步的近似值，从而得步的近似值，从而得到局部截断误差估计。到局部截断误差估计。（2）用同样的公式，用不同步长计算出）用同样的公式，用不同步长计算出n+1步的近似值，从而得到局步的近似值，从而得到局部截断误差估计。部截断误差估计。1.7 隐式法隐式法隐式法具有较好的绝对稳定性隐式法具有较好的绝对稳定性！只不过在使用隐式法的时候，需要进行迭代，或者使用预测校正计算只不过在使用隐式法的时候，需要进行迭代，或者使用预测校正计算格式。格式。现在学习的是第33页，共58页第二节第二节 泰勒级数法与龙格库塔法泰勒级数法与龙格库塔法对于方程：对于方程：取计算步长取计算步长为为 h，则则 ，将函数进行泰勒展开，将函数进行泰勒展开如函数如函数 y(x)有有p+1阶导数，容易得到阶导数，容易得到p阶泰勒级数展开法：阶泰勒级数展开法：现在学习的是第34页，共58页公式中的导数用下面公式计算：公式中的导数用下面公式计算：例子：例子：现在学习的是第35页，共58页步长步长0.1 步长步长0.01现在学习的是第36页，共58页龙格库塔法：龙格库塔法：对于常微分方程的数值解法，一个关键在于选择精度高的算法计对于常微分方程的数值解法，一个关键在于选择精度高的算法计算下面公式中的积分。算下面公式中的积分。要高精度的计算积分，常用的方法是适当增加计算节点，不妨要高精度的计算积分，常用的方法是适当增加计算节点，不妨设用设用 m 个节点计算上面积分，节点为个节点计算上面积分，节点为则积分为则积分为现在学习的是第37页，共58页将积分改写为：将积分改写为：则得公式：则得公式：这样的公式称为这样的公式称为显式龙格库塔显式龙格库塔公式。公式。现在学习的是第38页，共58页确定二确定二阶阶 RK 法法:系数满足：系数满足：此为四个未知数的三个方程，任意取此为四个未知数的三个方程，任意取 ，得，得现在学习的是第39页，共58页特别，取特别，取 ，得到，得到通常也称为变形欧拉法，也常写为通常也称为变形欧拉法，也常写为它具有二阶精度，也称为它具有二阶精度，也称为二阶二二阶二级级 RK 方法方法。现在学习的是第40页，共58页例子例子现在学习的是第41页，共58页三阶三级库塔法三阶三级库塔法局部截断误差为局部截断误差为4阶，整体截断误差为阶，整体截断误差为3阶阶现在学习的是第42页，共58页最常用的四级四阶公式，称为最常用的四级四阶公式，称为龙格库塔公式龙格库塔公式：局部截断误差为局部截断误差为5阶，整体截断误差为阶，整体截断误差为4阶。阶。现在学习的是第43页，共58页隐式龙格库塔法：隐式龙格库塔法：常用的有二阶常用的有二阶RK法：法：现在学习的是第44页，共58页例子：例子：用隐式用隐式二阶二阶RK法：法：用显式用显式二阶二阶RK法：法：现在学习的是第45页，共58页 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 continueEND现在学习的是第46页，共58页h=0.1现在学习的是第47页，共58页h=0.01现在学习的是第48页，共58页半隐式龙格库塔法：半隐式龙格库塔法：现在学习的是第49页，共58页最常用的二级三阶半隐式龙格库塔公式：最常用的二级三阶半隐式龙格库塔公式：现在学习的是第50页，共58页微分方程组微分方程组一阶方程组一阶方程组现在学习的是第51页，共58页龙格库塔公式龙格库塔公式现在学习的是第52页，共58页写成分量形式：写成分量形式：现在学习的是第53页，共58页例子：例子：现在学习的是第54页，共58页 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*(-y1(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.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现在学习的是第55页，共58页取取h=0.1现在学习的是第56页，共58页节点节点 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现在学习的是第57页，共58页高阶微分方程可以转换为一阶方程组进行求解。高阶微分方程可以转换为一阶方程组进行求解。令：令：现在学习的是第58页，共58页