计算方法第八章

计算方法第八章第一页，共五十八页，2022年，8月28日本章主要研究常微分方程初值问题的数值求解：通常，假设函数f关于第二个变量满足李普希茨条件（L条件），即为存在常数L>0，使得第二页，共五十八页，2022年，8月28日第一节一般概念1.1欧拉法及其简单改进方法：选择适当的节点，用差分近似微分，将方程离散化，从而求在这些节点上的函数值。欧拉方法第三页，共五十八页，2022年，8月28日例子：下面我们分别取步长为0.1与0.01进行计算，计算结果显示在下面的图中。第四页，共五十八页，2022年，8月28日步长为0.1的计算结果。第五页，共五十八页，2022年，8月28日步长为0.01的计算结果第六页，共五十八页，2022年，8月28日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第七页，共五十八页，2022年，8月28日

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

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第二十一页，共五十八页，2022年，8月28日同理，对于后退欧拉公式有预测－校正公式或改写为：第二十二页，共五十八页，2022年，8月28日用此法解前面的例子步长0.1步长0.01第二十三页，共五十八页，2022年，8月28日1.4误差估计定义：利用第n个节点的函数值，通过数值方法计算第n+1个节点的函数近似值，所引起的误差，称为n+1个节点上的局部截断误差。我们记为n+1个节点上函数的精确值，为数值解，则局部截断误差为：如局部截断误差为，称为具有p阶局部截断误差。第二十四页，共五十八页，2022年，8月28日欧拉方法的误差分析：省略余项得公式：第二十五页，共五十八页，2022年，8月28日完全类似的可以得到后退欧拉公式的局部截断误差为：欧拉中点公式的局部截断误差为：欧拉梯形公式的局部截断误差为：第二十六页，共五十八页，2022年，8月28日定义：由初值引起的第n个节点的近似值与精确值之间的误差称为第n个节点整体误差。定理：设下面求解微分方程的数值计算方法局部截断误差为p+1阶，且函数关于y

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

个节点计算上面积分，节点为则积分为第三十七页，共五十八页，2022年，8月28日将积分改写为：则得公式：这样的公式称为显式龙格－库塔公式。第三十八页，共五十八页，2022年，8月28日确定二阶R－K法:系数满足：此为四个未知数的三个方程，任意取，得第三十九页，共五十八页，2022年，8月28日特别，取，得到通常也称为变形欧拉法，也常写为它具有二阶精度，也称为二阶二级R－K方法。第四十页，共五十八页，2022年，8月28日例子第四十一页，共五十八页，2022年，8月28日三阶三级库塔法局部截断误差为4阶，整体截断误差为3阶第四十二页，共五十八页，2022年，8月28日最常用的四级四阶公式，称为龙格－库塔公式：局部截断误差为5阶，整体截断误差为4阶。第四十三页，共五十八页，2022年，8月28日隐式龙格－库塔法：常用的有二阶R－K法：第四十四页，共五十八页，2022年，8月28日例子：用隐式二阶R－K法：用显式二阶R－K法：第四十五页，共五十八页，2022年，8月28日

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第四十六页，共五十八页，2022年，8月28日h=0.1第四十七页，共五十八页，2022年，8月28日h=0.01第四十八页，共五十八页，2022年，8月28日半隐式龙格－库塔法：第四十九页，共五十八页，2022年，8月28日最常用的二级三阶半隐式龙格－库塔公式：第五十页，共五十八页，2022年，8月28日微分方程组一阶方程组第五十一页，共五十八页，2022年，8月28日龙格－库塔公式第五十二页，共五十八页，2022年，8月28日写成分量形式：第五十三页，共五十八页，2022年，8月28日例子：第五十四页，共五十八页，2022年，8月28日

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第五十五页，共五十八页，2022年，8月28日