计算方法第八章

计算方法第八章第1页，共58页，2023年，2月20日，星期四本章主要研究常微分方程初值问题的数值求解：通常，假设函数f关于第二个变量满足李普希茨条件（L条件），即为存在常数L>0，使得第2页，共58页，2023年，2月20日，星期四第一节一般概念1.1欧拉法及其简单改进方法：选择适当的节点，用差分近似微分，将方程离散化，从而求在这些节点上的函数值。欧拉方法第3页，共58页，2023年，2月20日，星期四例子：下面我们分别取步长为0.1与0.01进行计算，计算结果显示在下面的图中。第4页，共58页，2023年，2月20日，星期四步长为0.1的计算结果。第5页，共58页，2023年，2月20日，星期四步长为0.01的计算结果第6页，共58页，2023年，2月20日，星期四0.01 0.99005 0.990.1 0.90484 0.904380.2 0.81873 0.817910.3 0.74082 0.73970.40.67032 0.668970.41 0.66365 0.662280.59 0.55433 0.552680.6 0.54881 0.547160.90.40657 0.404730.91 0.40252 0.400680.99 0.37158 0.3697310.36788 0.36603第7页，共58页，2023年，2月20日，星期四

DOUBLEPRECISIONh,y(0：100)OPEN(20,FILE='OUTPUT.DAT',STATUS="UNKNOWN") h=1.0/100 y(0)=1.0 do10i=1,100 y(i)=y(i-1)*(1.0-h) write(20,*)i*h,y(i)10continue END第8页，共58页，2023年，2月20日，星期四精度更高的欧拉公式：方法：选择计算导数值的精度更好的差分格式。欧拉中点公式第9页，共58页，2023年，2月20日，星期四利用中点公式求解微分方程时，有一个问题，就是计算时需要两个迭代初值！对于这个问题，我们可以先用欧拉公式，通过给定的初值计算出第一个点的值，然后在利用这两个（第一和第二个点）的值进行计算，直到计算出全部节点上的值。下面，我们用中点公式求解刚才的例子，计算的步长取0.01，可以看到，计算的精度比较高。此时，计算公式为第10页，共58页，2023年，2月20日，星期四计算结果第11页，共58页，2023年，2月20日，星期四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页，2023年，2月20日，星期四1.2欧拉方法的其他改进微分方程数值解的关键在于对导数的处理，可以用差分来近似导数，也可以通过积分，将导数项化掉。对于方程：首先，作出划分设已经求出第n个节点的函数值，在区间上对方程两边积分容易看出，要求第n+1个节点的函数值，关键在于选择适当的积分公式计算积分！第13页，共58页，2023年，2月20日，星期四（1）如选择下矩形公式，则这正是前面的欧拉公式。（2）如选择上矩形公式，则这是所谓的后退欧拉公式。（3）如选择梯形公式，则这是所谓的欧拉梯形公式。第14页，共58页，2023年，2月20日，星期四直接利用已经求得的已知节点上的值计算未知节点上的函数值的算法称为显式法。例如：欧拉公式、欧拉中点公式计算未知节点上的函数值时，用到了未知节点上的函数值，这种算法称为隐式法。例如：后退欧拉法、欧拉梯形公式显然，利用隐式法求微分方程的数值解是，需要从表达式中反解未知节点上的函数值。第15页，共58页，2023年，2月20日，星期四1.3隐式法的具体计算：例如欧拉梯形公式用迭代法计算n+1步的值。（1）简单迭代收敛的条件：第16页，共58页，2023年，2月20日，星期四（2）牛顿迭代必须指出，在真正计算中，常用的是简单迭代，而且只迭代一步，迭代初值称为预测，迭代称为校正，这样组成的一组计算公式称为预测－－校正公式。第17页，共58页，2023年，2月20日，星期四预测－校正公式也称为改进的欧拉法，将上面的组合公式改写为：注意到，将上式进一步改写为：这是我们最终使用的计算格式。第18页，共58页，2023年，2月20日，星期四例子：取步长为0.1计算，结果如图。第19页，共58页，2023年，2月20日，星期四图：第20页，共58页，2023年，2月20日，星期四

DOUBLEPRECISIONh,y(0:10),ak1,ak2OPEN(20,FILE='OUTPUT1.DAT',STATUS="UNKNOWN") h=1.0/10 y(0)=1.0 do10i=1,10

ak1=-h*y(i-1) ak2=-h*(y(i-1)+ak1) y(i)=y(i-1)+(ak1+ak2)/2.010continuedo20i=0,10 write(20,*)i*h,y(i),exp(-i*h)20continue END第21页，共58页，2023年，2月20日，星期四同理，对于后退欧拉公式有预测－校正公式或改写为：第22页，共58页，2023年，2月20日，星期四用此法解前面的例子步长0.1步长0.01第23页，共58页，2023年，2月20日，星期四1.4误差估计定义：利用第n个节点的函数值，通过数值方法计算第n+1个节点的函数近似值，所引起的误差，称为n+1个节点上的局部截断误差。我们记为n+1个节点上函数的精确值，为数值解，则局部截断误差为：如局部截断误差为，称为具有p阶局部截断误差。第24页，共58页，2023年，2月20日，星期四欧拉方法的误差分析：省略余项得公式：第25页，共58页，2023年，2月20日，星期四完全类似的可以得到后退欧拉公式的局部截断误差为：欧拉中点公式的局部截断误差为：欧拉梯形公式的局部截断误差为：第26页，共58页，2023年，2月20日，星期四定义：由初值引起的第n个节点的近似值与精确值之间的误差称为第n个节点整体误差。定理：设下面求解微分方程的数值计算方法局部截断误差为p+1阶，且函数关于y

满足利普希茨条件，同时初值是准确的，则整体截断误差为p阶。欧拉公式、后退欧拉公式的整体截断误差为1阶。欧拉中点公式、欧拉梯形公式的整体截断误差为2阶。第27页，共58页，2023年，2月20日，星期四微分方程数值解法的进一步改进。再回到恒等式如果取作为节点，将被积函数用插值多项式来近似，用插值多项式带到积分中去求出积分，则可以得到所谓的亚当斯(Adams)显式公式局部截断误差：第28页，共58页，2023年，2月20日，星期四类似地，如果取作为节点，可得亚当斯(Adams)隐式公式局部截断误差：第29页，共58页，2023年，2月20日，星期四进一步，如果我们将恒等式中的积分区间改为，并在此区间上用辛普森公式，可得第30页，共58页，2023年，2月20日，星期四1.5绝对稳定性一个常微分方程数值解法称为绝对稳定，是指在某一步（如第一步）产生的误差（如计算机的存储误差），在计算中会逐步减小。称方程为试验方程，设计算步长为h，设在计算开始时产生误差（存储误差），此误差在以后会逐步减弱，我们称该算法相对于是绝对稳定的，的全体称为该算法的绝对稳定域。第31页，共58页，2023年，2月20日，星期四欧拉法的绝对稳定区域后退欧拉法的绝对稳定区域第32页，共58页，2023年，2月20日，星期四1.6局部截断误差的实用估计（1）用两种阶数相同的算法求解，计算出n+1步的近似值，从而得到局部截断误差估计。（2）用同样的公式，用不同步长计算出n+1步的近似值，从而得到局部截断误差估计。1.7隐式法隐式法具有较好的绝对稳定性！只不过在使用隐式法的时候，需要进行迭代，或者使用预测－校正计算格式。第33页，共58页，2023年，2月20日，星期四第二节泰勒级数法与龙格－库塔法对于方程：取计算步长为h，则，将函数进行泰勒展开如函数y(x)有p+1阶导数，容易得到p阶泰勒级数展开法：第34页，共58页，2023年，2月20日，星期四公式中的导数用下面公式计算：例子：第35页，共58页，2023年，2月20日，星期四步长0.1步长0.01第36页，共58页，2023年，2月20日，星期四龙格－库塔法：对于常微分方程的数值解法，一个关键在于选择精度高的算法计算下面公式中的积分。要高精度的计算积分，常用的方法是适当增加计算节点，不妨设用m

个节点计算上面积分，节点为则积分为第37页，共58页，2023年，2月20日，星期四将积分改写为：则得公式：这样的公式称为显式龙格－库塔公式。第38页，共58页，2023年，2月20日，星期四确定二阶R－K法:系数满足：此为四个未知数的三个方程，任意取，得第39页，共58页，2023年，2月20日，星期四特别，取，得到通常也称为变形欧拉法，也常写为它具有二阶精度，也称为二阶二级R－K方法。第40页，共58页，2023年，2月20日，星期四例子第41页，共58页，2023年，2月20日，星期四三阶三级库塔法局部截断误差为4阶，整体截断误差为3阶第42页，共58页，2023年，2月20日，星期四最常用的四级四阶公式，称为龙格－库塔公式：局部截断误差为5阶，整体截断误差为4阶。第43页，共58页，2023年，2月20日，星期四隐式龙格－库塔法：常用的有二阶R－K法：第44页，共58页，2023年，2月20日，星期四例子：用隐式二阶R－K法：用显式二阶R－K法：第45页，共58页，2023年，2月20日，星期四

DOUBLEPRECISIONh,y(0:100),yy(0:100) doubleprecisionakOPEN(20,FILE='OUTPUT5.DAT',STATUS="UNKNOWN") h=1.0/100 y(0)=1.0 yy(0)=1.0 do10i=1,100

y(i)=y(i-1)-(h/(1-h/2))*y(i-1) yy(i)=yy(i-1)*(1-h+h*h/2)10continuedo20i=0,100 write(20,100)i*h,y(i),yy(i),exp(-i*h)100format(1x,f4.2,f8.5,f8.5,f8.5)20continue END第46页，共58页，2023年，2月20日，星期四h=0.1第47页，共58页，2023年，2月20日，星期四h=0.01第48页，共58页，2023年，2月20日，星期四半隐式龙格－库塔法：第49页，共58页，2023年，2月20日，星期四最常用的二级三阶半隐式龙格－库塔公式：第50页，共58页，2023年，2月20日，星期四微分方程组一阶方程组第51页，共58页，2023年，2月20日，星期四龙格－库塔公式第52页，共58页，2023年，2月20日，星期四写成分量形式：第53页，共58页，2023年，2月20日，星期四例子：第54页，共58页，2023年，2月20日，星期四

DOUBLEPRECISIONh,y1(0:100),y2(0:100) doubleprecisionak1,ak2,ak3,ak4,bk1,bk2,bk3,bk4OPEN(20,FILE='OUTPUT6.DAT',STATUS="UNKNOWN") h=1.0/10 y1(0)=0.0 y2(0)=1.0 do10i=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.010continuedo20i=0,10 write(20,100)i*h,y1(i),i*h*exp(-i*h),y2(i),exp(-i*h)100format(1x,f4.2,f12.5,f12.5,f12.5,f12.5)20continue END第55页，共58页，2023年，2