第6章-常微分方程数值解法课件

数值分析数值分析第第第第6 6章章章章 常微分方程的数值解法常微分方程的数值解法常微分方程的数值解法常微分方程的数值解法 实际中，很多问题的数学模型都是微分方程。我们可实际中，很多问题的数学模型都是微分方程。我们可以研究它们的一些性质。但是，只有极少数特殊的方程有以研究它们的一些性质。但是，只有极少数特殊的方程有解析解。对于绝大部分的微分方程是没有解析解的。解析解。对于绝大部分的微分方程是没有解析解的。因此，需要研究和掌握微分方程的数值解法，即因此，需要研究和掌握微分方程的数值解法，即计算计算解域内离散点上的近似值的方法。解域内离散点上的近似值的方法。考虑考虑一阶一阶常常微分方程的微分方程的初值问题初值问题首先考虑微分方程解的存在性。首先考虑微分方程解的存在性。科大研究生学位课程数值分析数值分析只要只要 f(x,y)在在a,b R1 上连续，且关于上连续，且关于 y 满足满足 Lipschitz 条件条件，即存在与，即存在与 x,y 无关的常数无关的常数 L 使使对任意定义在对任意定义在 a,b 上的上的 y1(x)和和 y2(x)都成立，则上述问都成立，则上述问题题存在唯一解。存在唯一解。称称称称(6.1)(6.1)的解的解的解的解y y(x x)在节点在节点在节点在节点x xi i处的近似值处的近似值处的近似值处的近似值 y yi i y y(x xi i)a a x x1 1 x x2 2 .x xn n=b b.为其为其为其为其数值解数值解数值解数值解，方法称为，方法称为，方法称为，方法称为数值方法数值方法数值方法数值方法。注：注：注：注：考虑等距节点考虑等距节点考虑等距节点考虑等距节点:x xi i=a a+ihih，h h=(=(b b a a)/)/n n.从初始条件从初始条件从初始条件从初始条件y y(a a)=)=y y0 0出发，依次逐个计算出发，依次逐个计算出发，依次逐个计算出发，依次逐个计算y y1 1，y y2 2，y yn n的值，称为的值，称为的值，称为的值，称为步进法步进法步进法步进法。两种：单步法、多步法。两种：单步法、多步法。两种：单步法、多步法。两种：单步法、多步法。科大研究生学位课程数值分析数值分析 微分方程的微分方程的数值方法数值方法。求的是在一系列离散点列上，未知。求的是在一系列离散点列上，未知函数函数y(x)在这些点上的值的近似在这些点上的值的近似。为了考察数值方法提供的数值解，是否有实用价值，需要为了考察数值方法提供的数值解，是否有实用价值，需要知道如下问题：知道如下问题：步长充分小时，所得到的数值解能否逼近问题得真解；步长充分小时，所得到的数值解能否逼近问题得真解；即即收敛性问题收敛性问题 误差估计，主要考虑局部截断误差。误差估计，主要考虑局部截断误差。产生的舍入误差，在以后得各步计算中，是否会无限产生的舍入误差，在以后得各步计算中，是否会无限制扩大；制扩大；稳定性问题稳定性问题科大研究生学位课程数值分析数值分析6.1 欧拉欧拉Euler 法与法与改进欧拉法改进欧拉法1.欧拉法：欧拉法：x0 x1向前差商近似导数向前差商近似导数记为记为在在假假设设 yi=y(xi)，即即第第 i 步步计计算算是是精精确确的的前前提提下下，考考虑虑的截断误差的截断误差 i+1=y(xi+1)yi+1 称为称为局部截断误差局部截断误差若若某某算算法法的的局局部部截截断断误误差差为为O(hp+1)，则则称称该该算算法法有有p 阶阶精度。精度。亦亦称为称为欧拉折线法欧拉折线法 欧拉法的局部截断误差：欧拉法的局部截断误差：欧拉法具有欧拉法具有 1 阶精度。阶精度。科大研究生学位课程数值分析数值分析用向后差商公式代替导数项用向后差商公式代替导数项是隐格式，一般解不得显示表达；要迭代求解是隐格式，一般解不得显示表达；要迭代求解2.隐式隐式 Euler法法从上面知，隐式欧拉法也具有从上面知，隐式欧拉法也具有 1 阶精度。阶精度。（也可利用二元函数的泰勒展开）（也可利用二元函数的泰勒展开）近似有近似有科大研究生学位课程数值分析数值分析例例6.16.1：考虑初值问题：考虑初值问题：取步长取步长h h=0.1=0.1，并与准确解，并与准确解 比较。比较。解：因为解：因为x xi i=0.1=0.1i i，而，而f f(x x,y y)=-)=-y y+x x，故，故EulerEuler格式格式 y yi+i+1 1=y yi i+hf+hf(x xi i,y yi i)=(1-)=(1-h)y yi i+hx xi i (i i=0,1,)=0,1,)隐式隐式EulerEuler格式为格式为 y yi+i+1 1=y=yi i+hf+hf(x xi+i+1 1,y yi+i+1 1)=)=y yi i hy yi i+1+1+hxhxi i+1 +1 (i i=0,1,)=0,1,)整理得整理得 y yi+i+1 1=(=(y yi i+hxhxi i+1+1)/(1+/(1+h h)()(i i=0,1,)=0,1,)用用MatlabMatlab计算：计算：clear;xclear;x=0,y=0;=0,y=0;for i=1:5 for i=1:5 y=(1-h)*y=(1-h)*y+hy+h*x x=*x x=x+hx+h;end;end科大研究生学位课程数值分析数值分析6.1.2 欧拉格式的改进欧拉格式的改进对微分方程对微分方程计算定积分，则：计算定积分，则：用左用左矩形矩形公式公式用右用右矩形矩形公式，得公式，得得得欧拉公式和隐式欧拉公式欧拉公式和隐式欧拉公式科大研究生学位课程数值分析数值分析类似，可以算出梯形格式的误差估计式：类似，可以算出梯形格式的误差估计式：2阶的方法阶的方法所以，有格式为：所以，有格式为：梯形法是二阶、隐式单步的方法，要用迭代法求解。怎么求？梯形法是二阶、隐式单步的方法，要用迭代法求解。怎么求？局部截断误差局部截断误差数值积分使用梯形求积公式数值积分使用梯形求积公式上式称为梯形格式。上式称为梯形格式。科大研究生学位课程数值分析数值分析 改进欧拉格式改进欧拉格式/*modified Eulers Formula*/Step 1:先用显式欧拉格式作预测，算出先用显式欧拉格式作预测，算出),(1 i i i i y x f h y y+=+Step 2:再将再将 代入隐式梯形公式的右边作校正，得到代入隐式梯形公式的右边作校正，得到1+i y),(),(2 1 1 1+=i i i i i i y x f y x f h y y 考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：科大研究生学位课程数值分析数值分析注：注：此法亦称为此法亦称为预测预测-校正法校正法/*predictor-corrector method*/。可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个单单步步递推格式，比隐式公式的迭代求解过程简单。它的递推格式，比隐式公式的迭代求解过程简单。它的稳稳定性高于定性高于显式欧拉法。显式欧拉法。程序见程序见p123科大研究生学位课程数值分析数值分析6.2 龙格龙格-库塔法库塔法建立高精度的单步递推格式。建立高精度的单步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从(xi,yi)点出发，以点出发，以某一某一斜率斜率沿直线达到沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形点。欧拉法及其各种变形所能达到的最高精度为所能达到的最高精度为2阶阶。考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：科大研究生学位课程数值分析数值分析首先希望能确定系数首先希望能确定系数 1、2、p，使得到的算法格式有使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 Step 1:将将 K2 在在(xi,yi)点作点作 Taylor 展开展开将将改进欧拉法推广为：改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+=+=+Step 2:将将 K2 代入第代入第1式，得到式，得到科大研究生学位课程数值分析数值分析Step 3:将将 yi+1 与与 y(xi+1)在在 xi 点点的的泰勒泰勒展开作比较展开作比较要求要求 ，则必须有：，则必须有：这里有这里有 个未知个未知数，数，个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格-库库塔格式。塔格式。注意到，注意到，就是改进的欧拉法。就是改进的欧拉法。科大研究生学位课程数值分析数值分析称之为中点公式称之为中点公式,或可写为或可写为 若取1=0,则得2=1,p=1/2,公式为Q:Q:为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？科大研究生学位课程数值分析数值分析其中其中 i (i=1,m)，i (i=2,m)和和 ij(i=2,m;j=1,i 1)均为待均为待定系数，确定这些系数定系数，确定这些系数的步骤与前面相似。的步骤与前面相似。).,(.),(),(),(.1122112321313312122122111 +=+=+=+=mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 下面列出常用的三阶、四阶下面列出常用的三阶、四阶R-KR-K公式公式.科大研究生学位课程数值分析数值分析 四四阶经典阶经典R-KR-K公式公式 三阶三阶R-KR-K公式公式科大研究生学位课程数值分析数值分析 解解 四四阶阶标准标准R-KR-K公式为公式为例例例例6.4 6.4 用四阶经典用四阶经典用四阶经典用四阶经典R-KR-K方法求初值问题方法求初值问题方法求初值问题方法求初值问题y=y-2x/y ,0 x1y(0)=1的数值解的数值解,取步长取步长h=0.2.计算结果如下计算结果如下:科大研究生学位课程数值分析数值分析n nx xn ny yn ny y(x xn n)n nx xn ny yn ny y(x xn n)0 01 12 20.00.00.20.20.40.41.001.001.18321.18321.34171.34171.001.001.18321.18321.34161.34163 34 45 50.60.60.80.81.01.01.48331.48331.61251.61251.73211.73211.48321.48321.61251.61251.73211.7321 程序见p128科大研究生学位课程数值分析数值分析 求解初值问题求解初值问题的的单步显式方法可一统一写为如下形式单步显式方法可一统一写为如下形式 yn+1=yn+h(xn,yn,h)（*）对于Euler格式,有6.3 6.3 收敛性和稳定性收敛性和稳定性6.3.1 收敛性分析收敛性分析 y=(x,y),axb y(a)=y0 其中(x,y,h)称为增量函数增量函数.(x,y,h)=(x,y)对于改进的Euler格式,有 (x,y,h)=1/2(x,y)+(x+h,y+h(x,y)科大研究生学位课程数值分析数值分析 定义定义6.2 设y(x)是初值问题(6.1)的解,yn是单步法(*)产生的近似解.如果对任意固定的点xn,均有则称单步法(*)是收敛的.可见,若方法(*)是收敛的,则当h0时,整体截断误差整体截断误差en=y(xn)-yn将趋于零.注注 对形式形式简单的方程，可以由差分方程解的表达式的方程，可以由差分方程解的表达式 取极限取极限导出收出收敛性。性。科大研究生学位课程数值分析数值分析用用Euler法得近似解表达式法得近似解表达式例如例如对初初值问题：科大研究生学位课程数值分析数值分析 一般地，有如下定理 定理定理6.1 设单步方法(*)是p1阶方法,增量函数(x,y,h)在区域axb,-yn)的变化均不超过,则称此差分格式是绝对稳定绝对稳定的.讨论数值方法的稳定性,通常仅限于典型的试验方程 y=y (6.22)其中是复数且是复数且Re()0.(书上进一步限定0)在复平面上,当方法稳定时要求变量h的取值范围称为方法的绝对稳定域绝对稳定域,它与实轴的交集称为绝对稳定区间绝对稳定区间.科大研究生学位课程数值分析数值分析 将Euler方法应用于方程y=y,得到 设在计算yn时产生误差n,计算值yn=yn+n,则n将对以后各节点值计算产生影响.记ym=ym+m,mn,由上式可知误差m满足方程 m=(1+h)m-1=(1+h)m-n n,m n 对隐式单步方法也可类似讨论.如将梯形公式用于方程y=y,则有 yn+1=yn+h/2(yn+yn+1)yn+1=(1+h)yn 可见,若要|m|n|,必须且只须|1+h|1,因此Euler法的绝对稳定域为|1+h|1,绝对稳定区间是-2Re()h0.解出yn+1得 科大研究生学位课程数值分析数值分析类似前面分析,可知绝对稳定区域为由于Re()0,所以此不等式对任意步长h恒成立,这是隐式公式的优点.一些常用方法的绝对稳定区间为方方 法法方法的阶数方法的阶数稳稳 定定 区区 间间EulerEuler方法方法梯形方法梯形方法改进改进EulerEuler方法方法二阶二阶R-KR-K方法方法三阶三阶R-KR-K方法方法四阶四阶R-KR-K方法方法1 12 22 22 23 34 4(-2 ,0)(-2 ,0)(-(-,0),0)(-2 ,0)(-2 ,0)(-2 ,0)(-2 ,0)(-2.51 ,0)(-2.51 ,0)(-2.78 ,0)(-2.78 ,0)科大研究生学位课程数值分析数值分析 解解 因y0=1,计算得y10=1024,而y(1)=9.35762310-14.例例4 4 考虑初值问题考虑初值问题考虑初值问题考虑初值问题 y=-30y ,0 x1 y(0)=1取步长h=0.1,利用Euler方法计算y10y(1).y(x)=e-30 x 这是因为h=-3不属于Euler方法的绝对稳定区间.若取h=0.01,计算得y100=3.23447710-16.若取h=0.001,计算得y1000=5.91199810-14.若取h=0.0001,计算得y10000=8.94505710-14.若取h=0.00001,计算得y100000=9.315610-14.科大研究生学位课程数值分析数值分析 单步显式方法的稳定性与步长密切相关单步显式方法的稳定性与步长密切相关,在一种步长下是稳在一种步长下是稳定的差分公式定的差分公式,取大一点步长就可能是不稳定的取大一点步长就可能是不稳定的.收敛性是反映差分公式本身的截断误差对数值解的影响收敛性是反映差分公式本身的截断误差对数值解的影响;稳定稳定性是反映计算过程中舍入误差对数值解的影响性是反映计算过程中舍入误差对数值解的影响.只有即收敛又稳定只有即收敛又稳定的差分格式才有实用价值的差分格式才有实用价值.科大研究生学位课程数值分析数值分析6.4 6.4 线性多步方法（线性多步方法（线性多步方法（线性多步方法（AdamsAdams格式）格式）格式）格式）由于在计算yn+1时,已经知道yn,yn-1,及(xn,yn),(xn-1,yn-1),利用这些值构造出精度高、计算量小的差分公式就是线性多步法.一一.基于基于Taylor展开展开待定参数法构造线性多步方法待定参数法构造线性多步方法 r+1步线性多步方法的一般形式为当-10时,公式为隐式公式,反之为显式公式.参数i,i的选择原则是使方法的局部截断误差为 y(xn+1)-yn+1=O(hr+2)这里,局部截断误差局部截断误差是指,在yn-i=y(xn-i),i=0,1,r的前提下,误差y(xn+1)-yn+1.科大研究生学位课程数值分析数值分析 例 选取参数,0,1,2,使三步方法 yn+1=yn+h(0n+1n-1+2n-2)为三阶方法.解解 设yn=y(xn),yn-1=y(xn-1),yn-2=y(xn-2),则有 n=(xn,y(xn)=y(xn)n-1=(xn-1,y(xn-1)=y(xn-1)=y(xn-h)=y(xn)-hy(xn)+1/2h2y(xn)-1/6h3y(4)(xn)+O(h4)n-2=y(xn)-2hy(xn)+2h2y(xn)-4/3h3y(4)(xn)+O(h4)于是有 yn+1=y(xn)+h(0+1+2)y(xn)-h2(1+22)y(xn)+h3(1/21+22)y(xn)-h4/6(1+82)y(4)(xn)+O(h5)y(xn+1)=y(xn)+hy(xn)+1/2h2y(xn)+1/6h3y(xn)+1/24h4y(4)(xn)+O(h5)若使:y(xn+1)-yn+1=O(h4),只要,0,1,2满足:科大研究生学位课程数值分析数值分析 =1,0+1+2=1,1+22=-1/2,1+42=1/3解之得:于是有三步三阶显式差分格式 yn+1=yn+h/12(23n-16n-1+5n-2)科大研究生学位课程数值分析数值分析设pr(x)是函数(x,y(x)的某个r次插值多项式,则有因为 二二二二.利用数值积分构造线性多步方法利用数值积分构造线性多步方法利用数值积分构造线性多步方法利用数值积分构造线性多步方法 其中 由此,可建立差分格式 选取不同的插值多项式pr(x),就可导出不同的差分公式.下面介绍常用的Adams格式格式.1.Adams显式格式显式格式 科大研究生学位课程数值分析数值分析 设已求得精确解y(x)在步长为h的等距节点xn-r,xn上的近似值yn-r,yn,记k=(xk,yk),利用r+1个数据(xn-r,n-r),(xn,n)构造r次Lagrange插值多项式其中 由此,可建立差分格式 由于 hrj 则有 称之为r+1r+1步步AdamsAdams显式格式显式格式(外推格式外推格式)科大研究生学位课程数值分析数值分析下面列出几个带有局部截断误差主项的Adams显式公式 r=0 yn+1=yn+hn+(1/2)h2y(xn)(一阶(欧拉)格式）r=1 yn+1=yn+(h/2)(3n-n-1)+(5/12)h3y(xn)(二阶格式)r=2 yn+1=yn+(h/12)(23n-16n-1+5n-2)+(3/8)h4y(4)(xn)(三阶格式)r=3 yn+1=yn+(h/24)(55n-59n-1+37n-2-9n-3)+(251/720)h5y(5)(xn)(四阶阶格式)科大研究生学位课程数值分析数值分析 2.Adams隐式格式隐式格式 如果利用r+1个数据(xn-r+1,n-r+1),(xn+1,n+1)构造r次Lagrange插值多项式pr(x)在求积区间xn,xn+1,则可导出数值稳定性好的隐式公式,称为Adams隐式格式隐式格式,其一般形式为其中系数为 下面列出几个带有局部截断误差主项的Adams隐式格式(内插)r=0 yn+1=yn+hn+1-(1/2)h2y(xn)r=1 yn+1=yn+(h/2)(n+n+1)-(1/12)h3y(xn)r=2 yn+1=yn+(h/12)(5n+1+8n-n-1)-(1/24)h4y(4)(xn)r=3 yn+1=yn+(h/24)(9n+1+19n-5n-1+n-2)-(19/720)h5y(5)(xn)(一阶(隐式欧拉）格式)(二阶格式)(三阶格式)(四阶格式)科大研究生学位课程数值分析数值分析 3.Adams预估预估-校正格式校正格式 由显式公式提供一个预估值，再用隐式公式校正一次，求得数值解，称为预估校正方法预估校正方法。校正 yn+1=yn+(h/24)(9f(xn+1,pn+1)+19n-5n-1+n-2)一般预估公式和校正公式都采用同阶公式。例如：预估 pn+1=yn+(h/24)(55n-59n-1+37n-2-9n-3)称为四阶Adams预估校正格式.需要4个初值，实际计算时通常用四阶单步方法(如四阶R-K公式)为它提供起始值y1,y2,y3.例6.12 用四阶Adams预估校正格式求解初值问题 y=y-2x/y,0 x1 y(0)=1取步长h=0.1.科大研究生学位课程数值分析数值分析 解 用四阶R-K公式提供起始值,计算结果如下x xn nR-kR-k法法y yn n预估值预估值y yn n校正值校正值y yn n精确值精确值y y(x xn n)0 00.10.10.20.20.30.30.40.40.50.50.60.60.70.70.80.80.90.91 11 11.0954461.0954461.1832171.1832171.2649121.2649121.3415511.3415511.4140451.4140451.4830171.4830171.5489171.5489171.6121141.6121141.6729141.6729141.7315661.7315661.3416411.3416411.4142131.4142131.4832391.4832391.5491921.5491921.6124501.6124501.6733181.6733181.7320481.7320481 11.0954451.0954451.1832161.1832161.2649911.2649911.3416411.3416411.4142141.4142141.4832401.4832401.5491931.5491931.6124521.6124521.6733201.6733201.7320511.732051科大研究生学位课程数值分析数值分析6.5 一阶微分方程组和高阶微分方程一阶微分方程组和高阶微分方程6.5.1 一阶微分方程组一阶微分方程组一阶微分方程组的初值问题的一般形式为：一阶微分方程组的初值问题的一般形式为：科大研究生学位课程数值分析数值分析以两个方程为例以两个方程为例欧拉格式欧拉格式改进欧拉格式改进欧拉格式科大研究生学位课程数值分析数值分析经典四阶龙格经典四阶龙格-库塔格库塔格式式其中其中科大研究生学位课程数值分析数值分析Adams外插格式外插格式Adams预报预报-校正格式校正格式程序和实例见程序和实例见P140科大研究生学位课程数值分析数值分析6.5.2 高阶常微分方程高阶常微分方程科大研究生学位课程数值分析数值分析以二阶方程为例以二阶方程为例可化为一阶方程组可化为一阶方程组改进欧拉格式改进欧拉格式实例见实例见P142.科大研究生学位课程数值分析数值分析练习题练习题第第143页页 习题习题66-3-6.6,6.8,6.10-11,科大研究生学位课程