7--微分方程求解

第七章常微分方程的差分法,本章内容,7.1引言7.2欧拉法7.3龙格库塔法7.4亚当斯法7.5收敛性与稳定性小结,本章要求,1.熟悉各种Euler格式2.掌握龙格库塔法3.了解亚当斯法3.了解差分格式的收敛性和稳定性,7.1引言,本节内容一.问题提出二.两类求定解问题三.数值解法含义四.差商公式返回章节目录,一.问题提出有一个或多个导数及其函数的方程式称为微分方程，在工程中常遇到求解微分方程的问题。,7.1引言,很多微分方程的解不能用初等函数来表示，有时即使能够用解析式表示其解，但计算量太大而不实用(表达式过于复杂)。需要用数值方法来求解，一般只要求得到若干个点上的近似值或者解的简单的近似表达式(精度要求满足即可)。,7.1引言,二.两类求定解问题实际中求解常微分方程的所谓定解问题有两类:初值问题和边值问题定解指已知因变量和/或其导数在某些点上是已知的(约束条件)。1.边值问题约束条件为已知，在自变量的任一非初值上，已知函数值和/或其导数值，如常常可以将边解问题转化为初值问题求解。,7.1引言,2.初值问题,7.1引言,初值问题的常见解法,单步法：利用前一个单步的信息(一个点)，在y=f(x)上找下一点yi，有欧拉法，龙格库格法等。多步法：利用前面两步或两步以上的信息，在y=f(x)上找下一点yi，有改进欧拉法，亚当斯法等。,7.1引言,7.1引言,7.1引言,7.1引言,李普希兹条件:一条曲线上任意两点连线的斜率的绝对值都小于某一个数L,三.数值解法含义所谓数值解法,就是设法将常微分方程离散化,建立差分方程,给出解在一些离散点上的近似值。微分方程的数值解:设方程问题的解y(x)的存在区间是a,b,令a=x0x1,7.2欧拉方法,7.2欧拉方法,19/47,欧拉格式：/*EulersMethod*/,7.2欧拉方法,7.2欧拉方法,由于未知数yi+1同时出现在等式的两边，不能直接得到，故称为隐式/*implicit*/欧拉格式，而前者称为显式/*explicit*/欧拉格式。,一般先用显式计算一个初值，再迭代求解。,3.2欧拉方法,23/47,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列表计算如下：,7.2欧拉方法,25/47,7.2欧拉方法,显、隐式两种算法的平均,7.2欧拉方法,二.梯形格式,欧拉两步格式/*midpointformula*/,中心差商近似导数,7.2欧拉方法,7.2欧拉方法,/*predictor-correctormethod*/,改进的欧拉格式/*modifiedEulersmethod*/,Step1:先用显式欧拉格式作预测，算出,7.2欧拉方法,7.2欧拉方法,7.2欧拉方法,32/47,7.2欧拉方法,33/47,7.2欧拉方法,计算结果如下:,四.欧拉方法的精度,7.2欧拉方法,7.3龙格库塔方法,本节内容一.设计思想二.二阶龙格库塔格式三.高阶龙格库塔格式返回章节目录,7.3龙格库塔方法,一.设计思想,K*-平均斜率,7.3龙格库塔方法,7.3龙格库塔方法,39/47,7.3龙格库塔方法,折线法法的几何构造I,x1,x2,x3,Euler法(折线法),斜率近似,y(a),y(x1),y(x2),y(x3),前进Euler法(欧拉显式)的计算公式,显式公式,40/47,7.3龙格库塔方法,折线法法的几何构造II,x1,x2,x3,斜率近似,y(x1),y(x2),y(x3),y(x4),后退Euler法(欧拉隐式)的公式,隐式公式,取K1K2的算术平均值作为平均斜率,7.3龙格库塔方法,7.3龙格库塔方法,龙格库塔法的设计思想,7.3龙格库塔方法,二.二阶龙格库塔格式,斜率一定取K1K2的平均值吗？,步长一定是一个h吗？,7.3龙格库塔方法,7.3龙格库塔方法,7.3龙格库塔方法,这里有2个未知数,存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。,问题:为获得更高的精度，应该如何进一步推广？,7.3龙格库塔方法,三.高阶龙格库塔格式,3阶龙格-库塔法,7.3龙格库塔方法,7.3龙格库塔方法,最常用为四级4阶经典龙格-库塔法/*ClassicalRunge-KuttaMethod*/：,注：龙格-库塔法的主要运算在于计算Ki的值，即计算f的值。Butcher于1965年给出了计算量与可达到的最高精度阶数的关系：,由于龙格-库塔法的导出基于多项式展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h取小。,7.3龙格库塔方法,6,4,3,5,4,2,51/47,7.3龙格库塔方法,例：使用高阶R-K方法计算初值问题解:（1）使用三阶R-K方法(公式),52/47,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,53/47,7.3龙格库塔方法,（2）如果使用四阶R-K方法(公式),54/47,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,三阶Runge-Kutta方法,7.3龙格库塔方法,四阶(经典)Runge-Kutta方法,7.3龙格库塔方法,7.4亚当斯方法,本节内容一.设计思想二.二阶亚当斯格式三.亚当斯预报校正系统四.改进的亚当斯预报校正系统五.四阶亚当斯格式(略)返回章节目录,7.4亚当斯方法,一.设计思想由于龙格-库塔法在计算yn+1时,需要预报xn,xn+1区间上的若干斜率值。比如四阶龙格-库塔法，每一步都需要预报4个K值，即求4次函数f的值，计算量大。如果考虑利用yn+1之前的一系列“老”节点xn,xn-1，上的斜率值，则可以大大减少计算量。亚当斯法属于线性多步法,龙格-库塔法(4阶),亚当斯法(4阶),7.4亚当斯方法,欧拉格式和隐式欧拉格式属于一阶亚当斯格式,7.4亚当斯方法,二.二阶亚当斯格式,平均斜率K*,用xn-1,xn两点的斜率值加权平均作为平均斜率K*,7.4亚当斯方法,二阶显式亚当斯格式,7.4亚当斯方法,类似地,改用xn,xn+1两点的斜率值加权平均作为平均斜率K*，得到二阶隐式亚当斯格式,二阶隐式亚当斯格式(梯形格式),Step1:用Runge-Kutta法计算前k个初值；,Step2:用Adams显式计算预报值；,Step3:用同阶Adams隐式计算校正值。,注意：1-3步所用公式的精度必须相同。通常用Runge-Kutta法配合Adams公式。,三.亚当斯预报-校正系统/*Adamspredictor-correctorsystem*/,7.4亚当斯方法,四.改进的亚当斯预报-校正系统/*ImprovedAdamspredictor-correctorsystem*/,7.4亚当斯方法,注意预报和校正格式：,都具有二阶精度，能否将二者进行一次松弛计算，加工成具有三阶精度的格式？,考察差分格式对应的近似关系式,7.4亚当斯方法,7.4亚当斯方法,预报时Cn+1尚未知，用上一步的Pn-Cn代替此项,7.4亚当斯方法,改进的二阶亚当斯预报-校正系统:,7.5收敛性与稳定性,例：就初值问题考察欧拉显式格式的收敛性。,解：该问题的精确解为,欧拉公式为,对任意固定的x=xi=ih，有,7.5收敛性与稳定性,71/47,7.2欧拉方法,例：考察初值问题在区间0,0.5上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。,1.00002.00004.00008.00001.60001013.2000101,1.00002.50001016.25001021.56251023.90631039.7656104