数值分析常微分方程数值解法市公开课一等奖百校联赛特等奖课件.pptx

数值分析数值分析Numerical Analysis第八章常微分方程数值解法郑州大学硕士课程郑州大学硕士课程 （-年第一学期）年第一学期）ISCM,Beijing China1第1页第八章 常微分方程数值解法 8.1 8.1 引言引言8.2 8.2 欧拉欧拉(Euler)(Euler)法法8.3 8.3 改进欧拉改进欧拉(Euler)(Euler)方法方法8.4 8.4 单步法稳定性单步法稳定性ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第2页8.1 引言问题提出问题提出 倒葫芦形状容器壁上刻度问题倒葫芦形状容器壁上刻度问题.对于圆柱形状容器对于圆柱形状容器壁上容积刻度壁上容积刻度,能够利用圆柱体体积公式能够利用圆柱体体积公式其中直径D为常数.因为体积V与相对于容器底部任意高度H函数关系明确,所以在容器上能够方便地标出容器刻度,而对于几何形状不是规则容器,比如倒葫芦形状容器壁上怎样标出刻度呢?ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第3页8.1 引言下表是经过测量得到部分容器高度与直径关系下表是经过测量得到部分容器高度与直径关系.H 0 0.2 0.4 0.6 0.8 1.0D 0 0.11 0.26 0.56 1.04 1.17依据上表数据依据上表数据,能够拟合出倒葫能够拟合出倒葫芦形状容器图芦形状容器图,建立如图所表示建立如图所表示坐标轴后坐标轴后,问题即为怎样依据任问题即为怎样依据任意高度意高度x x标出容器体积标出容器体积V V刻度刻度,由由微元思想分析微元思想分析可知可知ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第4页8.1 引言其中其中x x表示高度表示高度,直径直径D D是高度是高度x x函数函数,记为记为D(x),D(x),所所以得到以下微分方程初值问题以得到以下微分方程初值问题只要求解上述方程只要求解上述方程,就可求出体积就可求出体积V V与高度与高度x x之间函之间函数关系数关系,从而可标出容器壁上容积刻度从而可标出容器壁上容积刻度,但问题是函但问题是函数数D(x)D(x)无解析表示式无解析表示式,我们无法求出其解析解我们无法求出其解析解.ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第5页8.1 引言 包含自变量、未知函数及未知函数导数或微分包含自变量、未知函数及未知函数导数或微分方程称为方程称为微分方程微分方程。在微分方程中。在微分方程中,自变量个数只自变量个数只有一个有一个,称为称为常微分方程常微分方程。自变量个数为两个或两。自变量个数为两个或两个以上微分方程叫个以上微分方程叫偏微分方程偏微分方程。微分方程中出现未。微分方程中出现未知函数最高阶导数阶数称为微分方程知函数最高阶导数阶数称为微分方程阶阶数。假如未数。假如未知函数知函数y y及其各阶导数及其各阶导数都是一次都是一次,则称它是则称它是线性线性,不然称为不然称为非线性非线性。ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第6页常微分方程常微分方程(ODEs 未知函数是一元函数未知函数是一元函数)偏微分方程偏微分方程(PDEs 未知函数是多元函数未知函数是多元函数)ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第7页同一个微分方程同一个微分方程,含有不一样初含有不一样初始条件始条件ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第8页当当x=0时时,y=1,可得可得c=1特特解解当当x=0时时,y=1,可得可得c=-1特解特解两边积分两边积分通解通解ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第9页8.1 引言 在高等数学中，对于常微分方程求解，给出了在高等数学中，对于常微分方程求解，给出了一些经典方程求解析解基本方法，如一些经典方程求解析解基本方法，如可分离变量法、可分离变量法、常系数齐次线性方程解法、常系数非齐次线性方程常系数齐次线性方程解法、常系数非齐次线性方程解法解法等。但能求解常微分方程依然是有限，大多数等。但能求解常微分方程依然是有限，大多数常微分方程是不可能给出解析解。常微分方程是不可能给出解析解。ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第10页8.1 引言 待求解问题待求解问题：一阶一阶常微分方程常微分方程初值问题初值问题/*Initial-Value Problem*/:解存在唯一性解存在唯一性（“常微分方程常微分方程”理论）：只要理论）：只要 f(x,y)在在a,b R1 上连续，且关于上连续，且关于 y 满足满足 Lipschitz 条件条件，即存在，即存在与与 x,y 无关常数无关常数 L 使使对任意定义在对任意定义在 a,b 上上 y1(x)和和 y2(x)都成立，则上述都成立，则上述IVP存存在唯一解在唯一解。ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第11页解析解法解析解法：（常微分方程理论）：（常微分方程理论）只能求解极少一类常微分方程；实际中给定问题不一定只能求解极少一类常微分方程；实际中给定问题不一定是解析表示式，而是函数表，无法用解析解法。是解析表示式，而是函数表，无法用解析解法。怎样求解怎样求解ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第12页ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第13页ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第14页8.2 欧拉(Euler)法 欧拉（欧拉（Euler）方法是解初值问题最简单）方法是解初值问题最简单数值方法。初值问题数值方法。初值问题解解y=y(x)y=y(x)代表经过点代表经过点 一条称之为微分一条称之为微分方程方程积分曲线积分曲线。积分曲线上每一点。积分曲线上每一点 切线斜率切线斜率 等于函数等于函数 在这点在这点值。值。ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第15页Euler法求解过程是法求解过程是:从初始从初始点点P0(即点即点(x(x0 0,y,y0 0)出发出发,作积分曲线作积分曲线y=y(x)y=y(x)在在P0点上点上切线切线 (其斜率为其斜率为 ),),与与x=xx=x1 1直线直线相交于相交于P1点点(即点即点(x(x1 1,y,y1 1),),得到得到y y1 1作为作为y(xy(x1 1)近似值近似值,如如上图所表示。过点上图所表示。过点(x(x0 0,y,y0 0),),以以f(xf(x0 0,y,y0 0)为为斜率切线方斜率切线方程为程为 当当 时时,得得 这么就取得了这么就取得了P P1 1点坐标。点坐标。ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第16页一样一样,过过点点P1(x x1 1,y,y1 1),),作积分曲线作积分曲线y=y(x)y=y(x)切线切线交直线交直线x=xx=x2 2于于P2点点,切线切线 斜率斜率 直线方程为直线方程为当当 时时,得得 ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第17页当当 时时,得得由此取得了由此取得了P P2 2坐标。重复以上过程坐标。重复以上过程,就可取得一系列就可取得一系列点点:P P1 1,P P1 1,P Pn n。对已求得点对已求得点以以 为斜率作直线为斜率作直线 取取ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第18页 从图形上看从图形上看,就取得了一条近似于曲线就取得了一条近似于曲线y=y(x)y=y(x)折线折线 。这么这么,从从x x0 0逐一算出逐一算出对应数值解对应数值解 ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第19页8.2 欧拉(Euler)法通常取通常取 (常数常数),),则则Euler法计算格式法计算格式 i=0,1,n (8.2)ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第20页8.2 欧拉(Euler)法还可用以下方法推导还可用以下方法推导EulerEuler格式：格式：数值微分数值微分 数值积分法数值积分法对微分方程离散，能够对微分方程离散，能够有各种思绪，但最基本有各种思绪，但最基本想法是想法是“以直代曲以直代曲”ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第21页8.2 欧拉(Euler)法(1)用差商近似用差商近似导导数数差分方程初差分方程初值问题值问题向前向前Euler方法方法ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第22页8.2 欧拉(Euler)法若若用向后差商近似用向后差商近似导导数数，即，即向后向后Euler方法方法ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第23页8.2 欧拉(Euler)法（2）用数）用数值积值积分方法分方法ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第24页8.2 欧拉(Euler)法若若对积对积分用梯形公式，分用梯形公式，则则得得梯形欧拉公式梯形欧拉公式ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第25页例例8.2.1 用欧拉法解初值问题用欧拉法解初值问题 取步长取步长h=0.2,h=0.2,计算过程保留计算过程保留4 4位小数位小数 解解:h=0.2,:h=0.2,欧拉迭代格式欧拉迭代格式 当当 k=0,x1=0.2时，已知时，已知x0=0，y0=1，有，有 y(0.2)y1=0.21(401)0.8当当 k=1,x2=0.4时，已知时，已知x1=0.2，y1=0.8，有，有 y(0.4)y2=0.20.8(40.20.8)0.6144当当 k=2,x3=0.6时，已知时，已知x2=0.4，y2=0.6144，有，有 y(0.6)y3=0.20.6144(4-0.40.6144)=0.4613 ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第26页8.2 欧拉(Euler)法解作解作为为微分方程初微分方程初值问题值问题数数值值解，即解，即以差分方程初以差分方程初值问题值问题ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第27页8.2 欧拉(Euler)法ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第28页8.2 欧拉(Euler)法x0 x1x2x3y y0h h h h h h 欧拉折线法欧拉折线法ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第29页解：解：Euler公式公式为为当当h=0.5时时ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第30页当当h=0.25时时ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第31页00.50.751.010.25h=0.5h=0.25ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第32页8.2 欧拉(Euler)法欧拉方法收敛性ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第33页8.2 欧拉(Euler)法局部截断误差称为局部截断误差ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第34页8.2 欧拉(Euler)法欧拉方法收敛性定义 若给定方法局部截断误差满足则称该方法是 P 阶，或称为含有 P 阶精度。ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第35页8.2 欧拉(Euler)法整体截断误差ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第36页8.2 欧拉(Euler)法欧拉方法收敛性ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第37页由此知，当 ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第38页8.2 欧拉(Euler)法 注 整体截断误差与局部截断误差关系：ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第39页8.2 欧拉(Euler)法 向后欧拉公式隐式欧拉法或向后欧拉法隐式欧拉法或向后欧拉法/*implicit Euler method or backward Euler method*/xn+1点向后差商近似导数点向后差商近似导数隐式或后退欧拉公式隐式或后退欧拉公式ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第40页8.2 欧拉(Euler)法 向后欧拉公式因为未知数因为未知数 yn+1 同时出现在等式两边，故称为同时出现在等式两边，故称为隐式隐式/*implicit*/欧拉公式，而前者称为欧拉公式，而前者称为显式显式/*explicit*/欧拉公欧拉公式。隐式公式不能直接求解，普通需要用式。隐式公式不能直接求解，普通需要用Euler显式公式显式公式得到初值，然后用得到初值，然后用Euler隐式公式迭代求解隐式公式迭代求解。所以隐式公。所以隐式公式较显式公式计算复杂，但稳定性好（后面分析）。式较显式公式计算复杂，但稳定性好（后面分析）。隐式欧拉公式中未知数隐式欧拉公式中未知数 yn+1 可经过以下迭代法求解：可经过以下迭代法求解：ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第41页8.2 欧拉(Euler)法 向后欧拉公式迭代法求隐式欧拉格式中迭代法求隐式欧拉格式中yn+1收敛性收敛性ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第42页 见上图，见上图，显然，这种近似也有一定误差，显然，这种近似也有一定误差，怎样预计这种误差怎样预计这种误差y(xn+1)yn+1？方法同上，基于方法同上，基于Taylor展开展开预计预计局部截断误差局部截断误差。不过注意，隐式公式中右边含有不过注意，隐式公式中右边含有f(xn+1,yn+1),因为因为yn+1不准确，所以不能直接用不准确，所以不能直接用y(xn+1)代替代替f(xn+1,yn+1)设已知曲线上一点设已知曲线上一点 Pn(xn,yn),过过该点作弦线，斜率为该点作弦线，斜率为(xn+1,yn+1)点方向场点方向场f(x,y)。若步长。若步长h充分小，充分小，可用弦线和垂线可用弦线和垂线x=xn+1交点近似曲交点近似曲线与垂线交点。线与垂线交点。几何意义xnxn+1PnPn+1xyy(x)ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第43页 隐式隐式欧拉法局部截断误差：欧拉法局部截断误差：ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第44页ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第45页 隐式隐式欧拉法局部截断误差：欧拉法局部截断误差：即隐式欧拉公式含有即隐式欧拉公式含有 1 阶精度。阶精度。隐式隐式欧拉法局部截断误差：欧拉法局部截断误差：ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第46页8.2 欧拉(Euler)法 向后欧拉公式比较欧拉显式公式和隐式公式及其局部截断误差显式公式隐式公式ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第47页 若将这两种方法进行算术平均，即可消除误差若将这两种方法进行算术平均，即可消除误差主要部分主要部分/*leading term*/而取得更高精度而取得更高精度,称为梯形法称为梯形法 梯形公式梯形公式/*trapezoid formula*/显、隐式两种算法显、隐式两种算法平均平均注：确有局部截断误差 ，即梯形公式具有2 阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第48页例例8.2.3 对初值问题对初值问题 证实用梯形公式求得近似解为证实用梯形公式求得近似解为 并证实当步长并证实当步长h h0 0时时,y,yn n收敛于准确解收敛于准确解证实证实:解初值问题梯形公式为解初值问题梯形公式为 整理成显式整理成显式 重复迭代重复迭代,得到得到 ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第49页公式公式局部截断误差局部截断误差精精度度显显隐隐稳稳定定性性步数步数欧拉显欧拉显式公式式公式1 1阶阶显显差差单步单步欧拉隐欧拉隐式公式式公式1 1阶阶隐隐好好单步单步梯形梯形公式公式2 2阶阶隐隐好好单步单步欧拉法小结欧拉法小结ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第50页8.3 改进欧拉(Euler)方法ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第51页8.3 改进欧拉(Euler)方法ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第52页8.3 改进欧拉(Euler)方法 显显式式欧欧拉拉公公式式计计算算工工作作量量小小,但但精精度度低低。梯梯形形公公式式虽虽提提升升了了精精度度,但但为为隐隐式式公公式式,需需用用迭迭代代法法求求解解,计计算算工工作作量量大大。综综合合欧欧拉拉公公式式和和梯梯形形公公式式便便可可得得到到改改进进欧拉公式。欧拉公式。结合已经有格式优点，以结合已经有格式优点，以得到计算方便、计算量降得到计算方便、计算量降低且精度保持数值格式低且精度保持数值格式ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第53页8.3 改进欧拉(Euler)方法 先用欧拉公式先用欧拉公式(8.2)求出一个初步近似值求出一个初步近似值,称称为为预预测测值值,它它精精度度不不高高,再再用用梯梯形形公公式式对对它它校校正正一一次次,即即迭迭代代一一次次,求求得得yn+1,称称为为校校正正值值,这这种种预预测测-校校正方法称为正方法称为改进欧拉公式改进欧拉公式：称称为为Euler公式与梯形公式公式与梯形公式预测预测校正系校正系统统。ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第54页8.3 改进欧拉(Euler)方法实际计实际计算算时时，常改写成以下形式，常改写成以下形式几何几何解释解释xnxn+1ABPn+1=(A+B)/2欧拉法欧拉法改进欧拉法改进欧拉法梯形法梯形法ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第55页predictorcorrectorISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第56页8.3 改进欧拉(Euler)方法 能能够够证证实实,改改进进欧欧拉拉公公式式精精度度为为二二阶阶。这这是是一个一步一个一步显式格式显式格式,它能够表示为嵌套形式它能够表示为嵌套形式。ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第57页例例8.3.18.3.1ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第58页ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第59页8.3 改进欧拉(Euler)方法ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第60页改进欧拉法算法ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第61页8.4 单步法稳定性 稳稳定定性性在在微微分分方方程程数数值值解解法法中中是是一一个个非非常常主主要要问问题题。因因为为微微分分方方程程初初值值问问题题数数值值方方法法是是用用差差分分格格式式进进行行计计算算，而而在在差差分分方方程程求求解解过过程程中中，存存在在着着各各种种计计算算误误差差，这这些些计计算算误误差差如如舍舍入入误误差差等等引引发发扰扰动动，在在传传输输过过程程中中，可可能能会会大大量量积积累累，对对计计算算结结果果准确性将产生影响准确性将产生影响。这就包括到算法稳定性问题。这就包括到算法稳定性问题。ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第62页例：例：考查初值问题考查初值问题 在区间在区间0,0.5上解。上解。分别用欧拉显、隐式格式和改进欧拉格式计算数值解。分别用欧拉显、隐式格式和改进欧拉格式计算数值解。0.00.10.20.30.40.5准确解准确解改进欧拉法改进欧拉法 欧拉隐式欧拉隐式欧拉显式欧拉显式 节点节点 xi 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.101 1.00002.5000 10 1 6.2500 10 21.5625 10 23.9063 10 39.7656 10 41.00002.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 78.4 单步法稳定性ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第63页8.4 单步法稳定性ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第64页若若某某算算法法在在计计算算过过程程中中任任一一步步产产生生误误差差在在以以后后计计算算中中都都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定绝对稳定/*absolutely stable*/。普通分析某算法稳定性时，为简单起见，只考虑普通分析某算法稳定性时，为简单起见，只考虑模型方程模型方程或或试验方程试验方程/*test equation*/8.4 单步法稳定性ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第65页记数值误差为：引进试验方程：ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第66页向前欧拉公式稳定性试验方程欧拉公式若每步计算有舍入误差，则 ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第67页当为复数时当为实数时ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第68页本章小结本章小结 本本章章介介绍绍了了常常微微分分方方程程初初值值问问题题基基本本数数值值解解法法。包包含含单单步步法法和和多多步步法法。单单步步法法主主要要有有欧欧拉拉法法、改改进进欧欧拉拉法法和和龙龙格格库库塔塔方方法法。多多步步法法是是亚亚当当姆姆斯斯法法。它它们们都都是是基基于于把把一一个个连连续续定定解解问问题题离离散散化化为为一一个个差差分分方方程程来来求求解解，是是一一个个步步进进式式方方法法。用多步法求常微分方程数值解可取得较高精度。用多步法求常微分方程数值解可取得较高精度。实际应用时，选择适当算法有一定难度，既要考虑算法实际应用时，选择适当算法有一定难度，既要考虑算法简易性和计算量，又要考虑截断误差和收敛性、稳定性简易性和计算量，又要考虑截断误差和收敛性、稳定性。ISCM,Beijing China/69 郑州大学硕士-年课程 数值分析 Numerical Analysis第69页