常微分方程的数值解法 (2)

第七章常微分方程的数值解法,计算机学院陈克建6学时,2/80,本章内容,7.1引言7.2欧拉方法7.3龙格库塔方法7.4阿达姆斯方法小结作业与实验,3/80,本章要求,1.熟悉Euler显格式,梯形法及Euler预校法;2.熟悉局部截断误差及绝对稳定性;3.掌握龙格库塔法。,4/80,7.1引言,本节内容一.问题提出二.两类求定解问题三.数值解法含义返回章节目录,5/80,一.问题提出有一个或多个导数及其函数的方程式称为微分方程，在工程中常遇到求解微分方程的问题。,7.1引言,6/80,很多微分方程的解不能用初等函数来表示，有时即使能够用解析式表示其解，但计算量太大而不实用(表达式过于复杂)。需要用数值方法来求解，一般只要求得到若干个点上的近似值或者解的简单的近似表达式(精度要求满足即可)。,7.1引言,7/80,7.1引言,二.两类求定解问题实际中求解常微分方程的所谓定解问题有两类:初值问题和边值问题定解指已知因变量和/或其导数在某些点上是已知的(约束条件)。1.边解问题约束条件为已知，在自变量的任一非初值上，已知函数值和/或其导数值，如常常可以将边解问题转化为初值问题求解。,8/80,7.1引言,2.初值问题,9/80,初值问题的常见解法,单步法：利用前一个单步的信息(一个点)，在y=f(x)上找下一点yi，有欧拉法，龙格库格法。预测校正法：多步法，利用一个以上的前点信息求f(x)上的下一个yi，常用迭代法，如改进欧拉法，阿当姆斯法。,7.1引言,10/80,7.1引言,11/80,7.1引言,12/80,P230定理1,7.1引言,13/80,三.数值解法含义所谓数值解法,就是设法将常微分方程离散化,建立差分方程,给出解在一些离散点上的近似值。微分方程的数值解:设方程问题的解y(x)的存在区间是a,b,令a=x0x1,7.2欧拉方法,17/80,7.2欧拉方法,18/80,P232,7.2欧拉方法,19/80,7.2欧拉方法,20/80,欧拉公式：/*EulersMethod*/,7.2欧拉方法,21/80,7.2欧拉方法,由于未知数yi+1同时出现在等式的两边，不能直接得到，故称为隐式/*implicit*/欧拉公式，而前者称为显式/*explicit*/欧拉公式。,一般先用显式计算一个初值，再迭代求解。,22/80,7.2欧拉方法,23/80,7.2欧拉方法,当i=0，即x0=0时，y1=0+0.2(1-00)=0.2000当i=1，即x1=0.2时，y2=0.2+0.2(1-0.20.2)=0.3920当i=2，即x2=0.4时，y3=0.392+0.2(1-0.40.392)=0.56064当i=3，即x3=0.6时，y4=0.56064+0.2(1-0.60.56064)=0.6933632当i=4，即x4=0.8时，y4=0.6933632+0.2(10.80.6933632)=0.782425088列表计算如下：,24/80,自己看看是不是,7.2欧拉方法,25/80,7.2欧拉方法,26/80,0)的近似值yk。,x1,x2,x3,积分曲线,45/80,7.3龙格库塔方法,折线法法的几何构造I,x1,x2,x3,Euler法(折线法),斜率近似,y(a),y(x1),y(x2),y(x3),前进Euler法的计算公式,显式公式,46/80,7.3龙格库塔方法,折线法法的几何构造II,x1,x2,x3,斜率近似,y(x1),y(x2),y(x3),y(x4),后退Euler法的公式,隐式公式,47/80,单步法：即利用前一个节点的函数值yi，计算后一个节点的函数值yi+1。目的：建立高精度的单步递推格式。单步递推法的基本思想是从(xi,yi)点出发，以某一斜率沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所能达到的最高精度为2阶。,7.3龙格库塔方法,二.2阶龙格库塔格式,48/80,斜率一定取K1K2的平均值吗？,步长一定是一个h吗？,7.3龙格库塔方法,49/80,7.3龙格库塔方法,首先希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得,Step1:将K2在(xi,yi)点作Taylor展开,50/80,Step2:将K2代入第1式，得到,7.3龙格库塔方法,51/80,Step3:将yi+1与y(xi+1)在xi点的泰勒展开作比较,要求，则必须有：,这里有3个未知数，2个方程。,存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。,注意到，就是改进的欧拉法。,7.3龙格库塔方法,52/80,其中i(i=1,m)，i(i=2,m)和ij(i=2,m;j=1,i1)均为待定系数，确定这些系数的步骤与前面相似。,问题:为获得更高的精度，应该如何进一步推广？,7.3龙格库塔方法,三.高阶龙格库塔格式,53/80,3阶龙格-库塔法,7.3龙格库塔方法,54/80,7.3龙格库塔方法,最常用为四级4阶经典龙格-库塔法/*ClassicalRunge-KuttaMethod*/：,55/80,注：龙格-库塔法的主要运算在于计算Ki的值，即计算f的值。Butcher于1965年给出了计算量与可达到的最高精度阶数的关系：,由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h取小。,7.3龙格库塔方法,56/80,7.3龙格库塔方法,例：隐式龙格-库塔法,而显式14阶方法的绝对稳定区域为,其中2阶方法的绝对稳定区域为,无条件稳定,57/80,7.3龙格库塔方法,例：使用高阶R-K方法计算初值问题解:（1）使用三阶R-K方法(公式),58/80,7.3龙格库塔方法,其余结果如下:,ixik1k2k3yi1.00000.10001.00001.10251.25551.11112.00000.20001.23451.37551.59451.24993.00000.30001.56241.76372.09221.42844.00000.40002.04042.34232.86581.66645.00000.50002.77683.25874.16341.9993,59/80,7.3龙格库塔方法,（2）如果使用四阶R-K方法(公式),60/80,7.3龙格库塔方法,其余结果如下:,ixik1k2k3k4yi1.00000.10001.00001.10251.11331.23511.11112.00000.20001.23461.37561.39211.56331.25003.00000.30001.56251.76391.79082.04231.42864.00000.40002.04082.34282.38922.78051.66675.00000.50002.77773.26003.34764.00572.0000,61/80,62/80,63/80,64/80,65/80,66/80,67/80,三阶Runge-Kutta方法,7.3龙格库塔方法,68/80,四阶(经典)Runge-Kutta方法,7.3龙格库塔方法,69/80,7.4阿达姆斯方法,本节内容一.一般线性多步法二.亚当姆斯显式公式三.亚当姆斯隐式公式四.亚当姆斯预测校正系统返回章节目录,70/80,7.4阿达姆斯方法,一.一般线性多步法由于在计算yn+1时,已经知道yn,yn-1,及(xn,yn),(xn-1,yn-1),利用这些值构造出精度高、计算量小的差分公式就是线性多步法。即用若干节点处的y及y值的线性组合来近似y(xi+1)。,其通式可写为：,当10时，为隐式公式;1=0则为显式公式。,71/80,7.4阿达姆斯方法,基于数值积分的构造法,72/80,7.4阿达姆斯方法,二.亚当姆斯显式公式/*Adamsexplicitformulae*/,Newton插值余项,/*显式计算公式*/,局部截断误差为：,73/80,7.4阿达姆斯方法,例：k=1时有,74/80,7.4阿达姆斯方法,注：一般有，其中Bk与yi+1计算公式中fi,fik各项的系数均可查表得到。,75/80,7.4阿达姆斯方法,76/80,7.4阿达姆斯方法,三.亚当姆斯隐式公式/*Adamsimplicitformulae*/,77/80,7.4阿达姆斯方法,小于Bk,较同阶显式稳定,78/80,7.4阿达姆斯方法,Step1:用Runge-Kutta法计算前k个初值；,Step2:用Adams显式计算预测值；,Step3:用同阶Adams隐式计算校正值。,注意：三步所用公式的精度必须相同。通常用经典Runge-Kutta法配合4阶Adams公式。,4阶Adams隐式公式的截断误差为,外推技术/*extrapolation*/,三.亚当姆斯预测-校正系统/*Ad