微分方程数值解法

第五章第五章：常微分方程数值解（NumericalMethodsforOrdinaryDifferentialEquations）引欧拉方§5.1§5.1引•本节内问题提定解问数值解法•在许多实际•在许多实际问题中，往往不能找出所需要的数关系，但是根据问题所提供的情况，有时例（1，2），且在该曲线上任意一曲线通过点点Mx,y处的切线斜率为横坐标的两倍，求这曲线的方程。什么是微分方•表示未知什么是微分方•表示未知函数、未知函数的导数及自变量之的关系的方程。（未知函数的导数必须出现––判断下列方程是否为微分方程xyxy3y微分方程的•微分方程中所出现的未知函数的最高阶导数的阶微分方程的•微分方程中所出现的未知函数的最高阶导数的阶数。dyx2yxy4y2y12y5ysin微分方程的一般形微分方程的一般形1、一阶微分方yfx,yFx,y,y或2、二阶微分方Fx,y,y,yyfx,y,或为什为什微分方程有着深刻而生动的实际背在自然界中，很多时候为了刻划客观对象的运动规律或变化规律，经常需要描述变量之但针对实际问题，我们通常很难直接找到这种函数关系，却容易建立起变量所满足的微如果方程可求解，则可以得到描述客观对象•初值问求一阶微分方程F初值问求一阶微分方程Fx,y,y0满足初值的特解的问题，称为一阶xx0分方程的初值问题x,y,y记y0xy=xy=x+yxy(0)其解析解为yx1本章重点讨论定解问题(初值问题f(x,y(本章重点讨论定解问题(初值问题f(x,y(x)定解条件（初始条件00f(x,是否能够找到定解问题的解取决如等x2y'sin(xy21xsiny',只要函数fy)适当光滑连续，且关于y只要函数fy)适当光滑连续，且关于y满足李普希兹(Lipschitz)条件即存在L，使得f(x,y)f(x,y)Ly由常微分方程理论知：初值问题的解必存在数数值解法所谓数值解法就是设法将常微分方程离散化建立差分方程,给出解在一些离散点上的近似值。微分方程的数值解:设方程问题的解y(x)的存在区间是[a,b令ax0x1xn=b其中hk=xk+1-xk如是等距节点h=(b-a)/n,h称为步长。似值用yk表示即yk≈y(xk这样y0,y1yn称为微本章：求(式1)的数值解，即解函y(x)在一些离散点x0x1x2本章：求(式1)的数值解，即解函y(x)在一些离散点x0x1x2y(x1)y(x2)y(xn)y(xn1)上的的近似解y()nnyy n112nyyf(xxxxxxh0单步法单步法利用前一个单步的信息(一个点)，y=f(x)上找下一点有欧拉法，龙格－库格法初解预测校正法多步法，利用一个以上的前点信息f(x)上的下一个斯法微分方程的数微分方程的数值解法需要解决的主要问基本方法有：有限差分法（数值微分）如何估计迭代公式的局部截断误差与整体误差如何保证迭代公式的稳定性与收敛性如果找不到解函数学界还数值解的思（1）将连续变量x[a,b]离散ax0如果找不到解函数学界还数值解的思（1）将连续变量x[a,b]离散ax0x1xkxn存在yy(x)x（2） y(xkk,,,工程师yk的光滑的振yk解的周期yx稳定数学界关解的混沌§5.2欧拉方§5.2欧拉方•本节内欧拉格梯形格––Euler预估—误差估计、收敛性和稳定欧拉方18世纪最杰出的数学家之一，13欧拉方18世纪最杰出的数学家之一，13时入读巴塞尔大学，1516岁获1727年-1741年（20岁-34岁）24岁晋1735年（28岁）右眼失明1741-1766（34岁-59岁）长，任职25年。在1741-1766（34岁-59岁）长，任职25年。在1766年应沙皇礼聘重回彼得堡，在1771年（64岁）明Euler是数学史上最多产的数学家，平均以每年800页的度写出创造性论文。他去世后，人们用35果74•一欧拉（Euler）式中hbn设节点为xaih(i,•一欧拉（Euler）式中hbn设节点为xaih(i,i方法一：alor展开y(i)()y(x)y(x)(x)xy(iiiiiiiy(i)y(xi)hf(xi,y(xi))2yihf(xi,yi(i,nyi式式2叫Euler显格式，可循环求解。方法二：数值微分法——实质Taylor展开法（方法三：数值积分法将方程fx,y)x]积分iixxi1i1yf(方法三：数值积分法将方程fx,y)x]积分iixxi1i1yf(x,y(xy(xi1)y(x)i1f(x,y(式i对右端积分采用左矩y(xi1)y(xi)hf(xi,yiyihf(xi,yiyi欧拉方f(x)dxba[f(a)欧拉方f(x)dxba[f(a)bf2abf(x)dx(ba)foxabbf(x)dx(ba)f右矩形公a中矩f(x)dx(ba)f(ab2ba左矩形a梯形公方法四：几何方法Euler法的几何意义过x0方法四：几何方法Euler法的几何意义过x0yx0x0y0作切线y斜率kyx0fx0yx0kf(x0,y0yy0fx0y0xx0)与xy0k(xx0x1求交点，纵坐标记为y1，(x0,y0过x1,y1)以fx1,y1)为斜率作直yy1fx1y1xx1)与xx2求交点，纵坐标记为y2y2yyy(yyy(20y欧拉公式：/*Euler’sMethody(x1)y(x0欧拉公式：/*Euler’sMethody(x1)y(x0y(x)0h记y(x)y(x)hy(x)yhf(x,y1000 1yihf(xi,yi(i0,...,nyiy(x1)y(y(x1)y(x0y(x)1hy(x1)y0hf(x1,y(x1hf(xi1,yi1 (i0,...,nyi由于未知数同时出现在等式的两边，不能直接得到，故称为/*implicit*/欧拉公式，而前者称为显/*explicit*/欧拉公式。例：用欧拉法解初值dy1例：用欧拉法解初值dy1(0x取步h0.2。计算过程4位小数解：因为fx,y)1xy，h欧拉公式为：y(xi1)yihf(xi,yiyi0.2(1xiyiyi(i当i=0，即x0=0时y1=0+0.2(1-0=当i=0，即x0=0时y1=0+0.2(1-0=0当i=1，即x1=0.2y2=0.2+0.2(1-当i=2，即x2=0.4y3=0.392+0.2(1-当i=3，即x3=0.6y4=0.56064+0.2(1-当i=4，即x4=0.8i01234501例：利用Euler方法求初值问题y12y2，0x1例：利用Euler方法求初值问题y12y2，0x1x的数值解，此问题的yx)解：此时的Euler公式为x/(1x21y2y2iii1xi,i 分别取步长h0.2,0.1,0.05，计算结果如下表：h h 7.00 0.---------------梯形格 )y(x)h[f（x梯形格 )y(x)h[f（x,y(x)),iii2yh[f（x,y)y, i 2(i,n(式显、隐式两种算法的平yihf(xi,yi(i0,...,nEuleryiyihf(xi,yi(i0,...,nEuleryiyihf(xi1,后退Euler yh[f(x,y)f((i0,...,n, i 2Euler格式的y(xn1p2后退Euler格式xn中点欧拉公式/*midpointformulay(中点欧拉公式/*midpointformulay(x)y(xy(x)20中心差商近似导1y(x2)y(x0)2hf(x1,y(x1yi12hf(xi,yii1,...,n三Euler校正法/*predictor-correctormethod•将Euler格式和梯形格式综合。三Euler校正法/*predictor-correctormethod•将Euler格式和梯形格式综合。改进的欧拉格精度低精度提高hyif(xi,yiStep1先用显式hyif(xi,yiStep1先用显式欧拉公式作预测Step2yh[f(x,y)f(,ii ii 2yhf(x,y)f,yhf(x,y(i0,...,nyii iiiii2hf(xi,yihhf(xi,yih [f(x,y)f,(式 i 2(i,n当yi已知时，有第一式预估出初值yi1，再代入第二改进的Euler方法也可以写成yh(K)改进的Euler方法也可以写成yh(K)ii122f(xi,yif(xih,yihK1KK,i,n27y，0xy例：求初值问题的数值解，取步长h0.1。（yx)(127y，0xy例：求初值问题的数值解，取步长h0.1。（yx)(12解：1.1yi0.2xi/yi(1)利用Euler方法,i yi0.05(K1K22x/y1iii2(xi(2)利用改进Euler方法y2i1yi1,illiii精确解012345678901111一些一些基本概yyn12hf(xn,yn§5.3本节§5.3本节内—库塔方•2高阶龙—库塔格—库塔格数值解法数值解法[a,b分成n得1ax0xn计算y(xk)(k>0)的近似值ykabEuler法(折线法—库塔方abEuler法(折线法—库塔方abyihf(xi,yiyi(i,n§7.3—库塔方折ab§7.3—库塔方折abyiyihf(xi1,yi1(i,n•2—库•2—库塔格单步递推法的基本思想xiyi点出发，以xi+1yi+1到的最高精度为2考察改进的欧拉法，ykhf(xk1,yk1h[f()y,f(考察改进的欧拉法，ykhf(xk1,yk1h[f()y,f(,kkkkkk2可将其改写为：斜率一定取K112yih2yi2K1f(xi,yif(xh,hK2ii1步长一定是一个h库塔方K yi1库塔方K yi1h[1K1Kf(xi,yif(xiph,yiphK1yiy(xiRiy(xi1)yi1O(h3Step1K2xiyiTaylorK2f(xiph,yiphK1f(xi,yi)phfx(xi,yi)phK1fy(xiy(xi)phy(xi)O(hyi)O(hy(x)d f(x,y)f(x,y)f(x,y)dyf(x,y)f(x,y)f(x,y) Step2K2代入第1yStep2K2代入第1yhy(x)[y(x)phy(x)O(h)]2yii1i2iiyi(12)hy(xi)2ph2y(xi)O(h3Step3yi+1yxi+1xiy()hy(x)phStep3yi+1yxi+1xiy()hy(x)phy(x)O(h23yii12i2i2y( )y(x)hy(x) y(x)O(h3iiii2O(h3y()Riii1 p1222-库塔格式注意到，p1,1就是改进的欧拉法212•高阶龙—库塔格•高阶龙—库塔格:yi1yih[1K12K2...mKm]K1f(xi,yi)K2f(xi2h,yi21K3f(xi3h,yi31hK132Kmf(ximh,ym1hK1m2hK2...m其中ii1,m)，ii2m(i=2,m;j=1,…,i1)y (K4KK3y (K4KK3hi126f(xi,yif(x,hhK)ii2221f(xih,yihK12hK2/*ClassicalRunge-KuttaMethod*/yi/*ClassicalRunge-KuttaMethod*/yi (K2K2KK4h3126f(xi,yif(xif(xih,h)221hK2h,22hK3f(xih,采用低阶算法而将步长h取小。每步须算Ki234567nO(h4O(hn2例：使用§7-K题塔方y0例：使用§7-K题塔方y0xy(0)取h解（1）使用三阶R-KnK1y200.1 (2201230120.1(yK)101236—库塔—库塔方其余结果§7: 库塔方（2）阶-KnK1y20库塔方（2）阶-KnK1y200.1 (20120.1 (30223 (2400.1(yK)1012346::—库塔方 —库塔方yh(K)i—库塔方yh(K)ii1236 f(xi,yif(xh,hKK2ii122Kf(xih,yih(2K2K1y(x 0yh(K)ii12346yh(K)ii12346f(xi,yi1h,2h,h2)f(iKh2)f(f(i32K4h,yihK3y(x0n级显式Runge-Kuttanyi1bjn级显式Runge-Kuttanyi1bjkjjf(x,ykiijhajmkm),j2,3,...,f(xicjh,jajmcj二阶Runge-Kutta方法的局部截断误二阶Runge-Kutta方法的局部截断误只能达三阶Runge-Kutta方法的局部截断误O(h4只能达四阶Runge-Kutta方法的局部截断误只能达五阶Runge-Kutta方法的局部截断误O(h6只能达几种方法的数值几种方法的数值计例Euler改进的Euler例四阶经典Lunge-Kutta例四阶经典Lunge-Kutta几种方法的结果与误x精确Euler0几种方法的结果与误x精确Euler01Runge-Kuttay(xn1y(xn,Runge-Kuttay(xn1y(xn,h存在01,y(xn1)y(xn)=ynh由所给的方程y'f(xy)y(xn1)=y(xn)+hf(xnh,y(xnfh,y(xnh)),故可以将K取K**区间[xn,xn1]对平均斜率给出不同的算法，就能有（*对平均斜率给出不同的算法，就能有（*）导相应的计算公式譬如取简单地取xn处的斜率值K1f(xnyn)作平均斜率值K*就得到Euler公式,精度当然很低.作为K*就得到改进Euler公式.这就给我们以启示,可以通过多预测几个点的斜量值,然后取加值作为K*,得到更高精度的公式所求的方法称为2阶RungeKutta方法ynh所求的方法称为2阶RungeKutta方法ynh(1K12K2f(x,)1nnf(xnph,ynK21112,222阶RungeKutta方法需要计算2个函数值令p12,则1=0，2令p12,则1=0，2ynf(x,)1nnh1,2 hK1)f(xn1=1/令p23,3/3Kh(11=1/令p23,3/3Kh(1y1n24)4f(x,1nn22K2f(xn+3h,1K3这个算法称为：Heun方单步方法的收敛性与稳and收敛/*单步方法的收敛性与稳and收敛/*Convergency定xxix0ih，h0(同收敛的i)时yiy(xi)，则称该算法y例：就初值问收敛性考察欧拉显式格式y(0)0y(x)解：该问题的精确解yy(x)解：该问题的精确解yihyi(1h)yi欧拉公式lim(1h)1/xxiih，y0h)1/i/]0yexiy(x0i稳定/*Stabilityy(x)30y(例：稳定/*Stabilityy(x)30y(例：在区间[0,y(0)上的解.分别用欧拉显、隐式格式和改进的欧拉式计算数y欧拉欧拉隐改进欧拉精确xi定若某算法在计算过程定若某算法在计算过程中任一步产生的误差在以后的计算中都逐步衰减，则称该算法是/*absolutelystable*/一般分析时为简单起，方equationy当步长取初值产生误差时，将某当步长取初值产生误差时，将某算法应用于上式，并假设只y0，则若此误差以后逐步衰减的全体构就称该算法相对于绝对稳定区域绝对稳定我们称算法比算法稳定，就是的绝对稳B的大hyi(1h)i1y0例：考察显式欧拉hyi(1h)i1y0例：考察显式欧拉0y0(1h)i1(1h)i1--0i1yi1由此可见，要保证初始误差以后逐步衰减h|1h|hyiyi例：考察隐式欧拉hyiyi例：考察隐式欧拉1yii1hi1i1h可见绝对稳定区域为|1h|注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。 例：隐式龙格-库塔yih[1K1...mKm例：隐式龙格-库塔yih[1K1...mKmf(xh,yhK...)jiji1jm(j1,...,yi的绝对稳其中2阶方Kf(xh,yhK1ii122无条件稳0而4阶方法的绝对稳定区域而4阶方法的绝对稳定区域------§5.4亚§5.4亚当姆斯方•本节内一般线性多亚当姆斯显式亚当姆斯隐式亚当姆斯预测—•一般线性多由于在计算yn+1时已经知道yn•一般线性多由于在计算yn+1时已经知道ynyn-1及(xn-1,yn-1利用这些值构造出精度高、计算即用若干节点处yy’值的线性组合来近fjf(xj,yj )0yi1yi1fi1fifi1...fik10时，为隐式公式;1=0则为显式公式yfxyxyfxyxixi1]y(xi1)y(xi)f(x,y( fx,yx))dxkxiyi过yi同的计算公式•/*Adamsexplicitformulaefi,fi1,...,fi•/*Adamsexplicitformulaefi,fi1,...,fiNk(xith)t[0,1]11f(x,y(x))dx(xth)R(xth)N kix00iNewton1yi N(xth)dtyiki01局部截断误差为：Riy(xi1)yi1 R(xth)ki0例：k1N1(xith)1fit fit(fifi1例