计算方法4常微分方程数值解法.ppt

常微分方程的数值解法Numerical Solutions to Ordinary Differential Equations,对象,一阶常微分方程初值问题:,一阶常微分方程组初值问题:,高阶常微分方程初值问题:,(4.1),一阶常微分方程初值问题:,实际工程技术、生产、科研上会出现大量的微分方程问题很难得到其解析解，有的甚至无法用解析表达式来表示，因此只能依赖于数值方法去获得微分方程的数值解。,用数值方法，求得y(x)在每个节点xk 的值y(xk)的近似值，用yk 表示，即yk y(xk)，这样y0,y1,.,yn称为微分方程的数值解求y(x)求y0,y1,.,yn,？,微分方程的数值解法：不求y=y(x)的精确表达式,而求离散点x0,x1,xn处的函数值设(4.1)的解y(x)的存在区间是a,b，初始点x0=a，取a,b内的一系列节点x0,x1,.,xn。a=x0 x1 xn=b，一般采用等距步长,思路,计算过程：,方法：采用步进式和递推法,将a,bn等分,a=x0 x1 xn=b,步长h=(b-a)/n,xk=a+kh,怎样建立递推公式？Taylor公式数值积分法,4.1 Euler 公式,思想:用向前差商近似代替微商.,(4.2),欧拉公式(Euler Scheme）,几何意义,y(x)过点P0(x0,y0)且在任意点(x,y)的切线斜率为f(x,y)y(x)在点P0(x0,y0)的切线方程为：y=y0+f(x0,y0)(x-x0)在切线上取点P1(x1,y1)y1=y0+f(x0,y0)(x1-x0)y1正是Euler 公式所求,4.类似2，过P1以f(x1,y1)为斜率作y(x)的切线，在其上取点 P2(x2,y2)，依此类推,5.折线P0 P1 P2 Pn作为曲线y(x)的近似,欧拉折线法,思想:用向后差商近似代替微商.,欧拉法（续）,用隐式欧拉法，每一步都需解方程（或先解出yn+1的显式表达式），但其稳定性好。,隐式欧拉公式,(4.3),整体误差ek=y(xk)-yk，下面对其加以分析,y1=y0+hf(x0,y0)=1+0.1(1-0/1)=1.1y2=y1+hf(x1,y1)=1.1+0.1(1.1-20.1/1.1)=1.191818y3=y2+hf(x2,y2)=1.277438其精确解为,欧拉法（续）,思想:用中心差商近似代替微商.,注：计算时，先用欧拉法求出y1，以后再用二步欧拉法计算。,二步欧拉法,(4.4),欧拉法（续）,公式,单步否,显式否,单步,显式,单步,隐式,二步,显式,截断误差y(xn+1)-yn+1,截断误差,Def4.1 设y(xn)是(4.1)式的精确解，yn是(4.2)式欧拉法得到的近似解，称y(xn)-yn为欧拉法的整体截断误差.,Def4.3 若某算法的局部截断误差为O(hp+1),称该算法有p阶精度.,Def4.2 假设yn=y(xn),即第n步计算是精确的前提下,称Rn+1=y(xn+1)-yn+1为欧拉法的局部截断误差.,分析：证明其局部截断误差为O(h2),可通过Taylor展开式分析。,证明：Euler 公式为,令yn=y(xn)，下证：y(xn+1)-yn+1=O(h2),由 y(x)=f(x,y(x),定理4.4 欧拉法的精度是一阶。,二步欧拉法的局部截断误差,Def4.5 假设yn=y(xn),yn1=y(xn1),称Rn+1=y(xn+1)-yn+1为二步欧拉法的局部截断误差.,定理4.6 隐式欧拉法的精度是一阶,二步欧拉法的精度是二阶。,证明：对二步欧拉法进行证明，考虑其局部截断误差，,令yn=y(xn),yn1=y(xn1),将上两式左右两端同时相减：,二步欧拉法的局部截断误差为O(h3)，其精度是二阶。,数值积分法,对右端的定积分用数值积分公式求近似值：,（1）用左矩形数值积分公式：,（2）用梯形公式：,梯形公式,梯形公式:将显示欧拉公式,隐式欧拉公式平均可得,梯形公式是隐式、单步公式,其精度为二阶,证：令yn=y(xn),由Talor公式有,分析：梯形公式是隐式公式，证明其局部截断误差为O(h3)要用到 二元函数的Taylor公式。,f(xn+1,yn+1)=f(xn+1,y(xn+1)+(yn+1-y(xn+1)=f(xn+1,y(xn+1)+fy(xn+1,)(yn+1-y(xn+1)，(xn，xn+1)=y(xn+1)+fy(xn+1,)(yn+1-y(xn+1)=y(xn)+hy”(xn)+O(h2)+fy(xn+1,)(yn+1-y(xn+1)=f(xn,yn)+hy”(xn)+fy(xn+1,)(yn+1-y(xn+1)+O(h2)又y(xn+1)=y(xn+h)=y(xn)+hy(xn)+h2y”(xn)/2+O(h3)=yn+hf(xn,yn)+h2y”(xn)/2+O(h3)=yn+hf(xn,yn)/2+hf(xn,yn)+hy”(xn)/2+O(h3),定理4.7 梯形公式的精度是二阶。,从而y(xn+1)=yn+1+h fy(xn+1,)y(xn+1)-yn+1/2+O(h3)y(xn+1)-yn+1=h fy(xn+1,)y(xn+1)-yn+1/2+O(h3)y(xn+1)-yn+1=O(h3)/1-hfy(xn+1,)/2=O(h3)梯形公式的截断误差为O(h3)，其精度是2阶。,f(xn+1,yn+1)=f(xn,yn)+hy”(xn)+fy(xn+1,)(yn+1-y(xn+1)+O(h2)y(xn+1)=yn+hf(xn,yn)/2+hf(xn,yn)+hy”(xn)/2+O(h3)=yn+hf(xn,yn)/2+h f(xn+1,yn+1)-fy(xn+1,)(yn+1-y(xn+1)+O(h2)/2+O(h3)=yn+hf(xn,yn)+f(xn+1,yn+1)/2+h fy(xn+1,)(y(xn+1)-yn+1)/2+O(h3),解:取h=0.01 x0=0 y0=y(0)=1 则 y(0.01)y1 f(x,y)=y，由梯形公式,基于稳定性考虑,解析解,欧拉公式的比较,4.2 改进的Euler法,Euler公式 计算量小，精度低梯形公式 计算量大，精度高,综合两个公式，提出预报校正公式：,预报,校正,改进的Euler法,写成嵌套公式：,平均化形式：,例4.4 用改进的Euler法解初值问题在区间0,0.4上，步长h=0.1的解，并比较与精确解(y=1/(1-x)的差异。,解：Euler法的具体形式为：yn+1=yn+hyn2，改进的Euler法的具体形式为：,x0=0,h=0.1,则 x1=0.1,x2=0.2,x3=0.3,x4=0.4 计算y1:yp=y0+0.1y02=1+0.112=1.1 yc2=1.121 y1=(1.1+1.121)/21.1118,同样可求y2,、y3、,y4,注：（1）令y(xn)=yn，可推导改进的Euler法的局部截断误差 为O(h3),具有二阶精度。,（2）改进的Euler法也可写成如下平均化形式,4.3龙格库塔方法,如何构造高阶的方法？,(4.5),为了构造函数使得(4.5)式成为高阶方法，Taylor,对于一般的显式单步法,若讨论其精度（阶数），即考虑,令,则有,考虑到,上述方法求高阶微分，实际上不实用,改进的Euler公式:,Euler公式：yn+1=yn+hf(xn,yn),写成,精度：一阶,精度：二阶,由Lagrange中值定理，,称为y(x)在xn,xn+1上的平均斜率,取,Euler公式,取,改进Euler公式,Euler公式用一个点的值k1作为k*的近似值，而改进的Euler公式用二个点的值k1和k2的平均值作为k*近似值，改进的Euler法比Euler法精度高。,(4.6),Runge-Kutta法的思想：在xn,xn+1内多预报几个点的ki值并用其加权平均值作为k*近似值。而构造出具有更高精度的计算公式。,多个斜率加权平均可提高精度,Runge-Kutta法一般形式：,以k1与k2的加权平均来近似k*，设,其中c1,c2,2,b21为待定参数。适当选取参数，使(*)式的精度为二阶，即使其局部截断误差为O(h3),令y(xn)=yn，由泰勒公式：,二阶龙格库塔方法,(*),由多元函数的泰勒公式和(*)式,则有,比较(a)与(b)要使 y(xn+1)-yn+1=O(h3),(a),(b),上述方程组有四个未知量，只有三个方程，有无穷多组解。取任意一组解便得一种二阶龙库公式。,当c1=c2=1/2,a2=b21=1时二阶Runge-Kutta公式为,yn+1=yn+k1/2+k2/2k1=hf(xn,yn)k2=hf(xn+h,yn+k1),此即改进Euler法,取c2=0，c2=1，a2=1/2，b21=1/2,yn+1=yn+k2 k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2),此为中点法或变形的 Euler公式,三阶龙格库塔法是用三个值k1,k2,k3的加权平均来近似k*，即有,yn+1=yn+c1k1+c2k2+c3k3k1=hf(xn,yn)k2=hf(xn+a2h,yn+b21k1)k3=hf(xn+a3h,yn+b31k1+b32k2),要使其具有三阶精度，必须使局部截断误差为O(h4),类似二阶龙格库塔法的推导，c1,c2,c3,a2,a3,b21,b31,b32应满足,由该方程组任意解可得三阶龙格-库塔公式,例：Kutta公式,类似可推出四阶龙格-库塔公式，常用的有 例：经典Runge-Kutta法,局部截断误差 O(h5),还有:Gill公式及m(m4)阶龙格库塔法。m4时：计算量太大，精确度不一定提高，有时会降低。,Gill公式,节省存储单元控制舍入误差,对于经典的四阶Runge-Kutta法给出如下算法：,算法4.2求解：dy/dx=f(x,y)axb y(a)=y0,Step 1:输入a,b,y0 及NStep 2:(b-a)/N=h,a=x,y0=yStep 3:输出(x,y)Step 4:For I=1 T0 Nhf(x,y)=k1hf(x+h/2,y+k1/2)=k2hf(x+h/2,y+k2/2)=k3hf(x+h,yk3)=k4y+(k1+2k2+2k3+k4)/6=yx+h=x输出(x,y)Step 5:停止,例4.3用四阶经典RungeKutta方法解初值问题：,(1)求，,(2)求，,自适应龙格库塔法,用户提出问题I：,问题：如何判断|y(xn)-yn|精度值y(xn)未知。：如何取h=?解：如用p阶龙格库塔法计算，局部截断误差为O(hp+1),如 xn-xn+1,令 yn=y(xn)yn+1(h)则 y(xn+1)-yn+1(h)chp+1,步长折半xnxn+h/2xn+1分两步计算y(xn+1)的近似值yn+1(h/2)。则 y(xn+1)-yn+1(h/2)2c(h/2)p+1,定理：对于问题I若用P阶龙格库塔法计算y(xn+1)在步长折半前后的近似值分别为yn+1(h),yn+1(h/2)则有误差公式,注：10 误差的事后估计法 20 停机准则：(可保证|y(xn+1)-yn+1(h/2)|),解：h取大，局部截断误差chp+1大，不精确 h取小，运算量大（步多），舍入误差积累大解决策略：变步长龙格库塔法 If（）将步长折半反复计算，直至为止，最后一次运算的前一次计算结果即为所需。,4.5 收敛与稳定性,对于常微分方程初值问题(4.1),例：Euler法，改进的Euler法，龙格库塔法都是单步法。,数值解法：,单步法：计算yn+1时只用到前一步的结果yn。,显式单步法：yn+1=yn+h(xn,yn,h)(x,y,h)为增量函数，它依赖于f,仅是xn,yn,h的函数,Def:若某数值方法对于任意固定的xn=x0+nh,当 h-0(n-)时，yn-y(xn)，则称该法是收敛的。,即(xn=x0+nh为固定值),改进Euler法的收敛性,判定改进Euler法的收敛性：yn+1=yn+hf(xn,yn)+f(xn+h,yn+hf(xn,yn)/2 则(x,y,h)=f(x,y)+f(x+h,y+hf(x,y)/2,若：yo=y(x0),f(x,y)关于y满足李-条件:,则：,限定hh0，则,即(x,y,h)满足李普希兹条件,由Th改进的Euler法收敛,同样可验证，若f(x,y)关于y满足李普希兹条件,且 y0=y(x0)时，Euler法，标准四阶龙格库法也收敛。,稳定性,在讨论收敛性时，一般认可：数值方法本身计算过程是准确的。实际上，并非如此：初始值y0有误差0y0-y(x0).计算的每一步有舍入误差。这些初始误差在计算过程的传播中，是逐步衰减的，还是恶性增长，这就是稳定性问题。,Def4.1若一种数值方法在节点xn处的数值解yn的扰动n0，而在以后的各节点值ym(mn)上有扰动m.当|m|n|，(m=n+1,n+2,)，则称该数值算法是稳定的。,分析某算法的数值稳定性，通常考察模型方程 y=y，(0),Def4.1：设在节点xn处用数值法得到的理想数值解为yn，而实际计算得到的近似值为，称 为第n步数值解的扰动,Euler法的稳定性,Euler法：yn+1=yn+hf(xn,yn)考察模型方程 y=y，(0)即yn+1=(1+h)yn,假设在节点值yn上有扰动n，在yn+1上有扰动n+1，且n+1仅由n引起(计算过程不再引进新的误差),Euler法稳定,Euler法的稳定的条件为 0h-2/,隐式Euler法稳定性,隐式Euler法：yn+1=yn+hf(xn+1,yn+1)对于模型方程 y=y(yn+1=yn/(1-h),yn+1的扰动,0 上式均成立，隐式Euler法无条件稳定,梯形公式稳定性,梯形公式 yn+1=yn+hf(xn,yn)+f(xn+1,yn+1)/2，设模型方程为y=y(0)则 yn+1=yn+hyn+yn+1/2,0时恒成立 梯形公式恒稳定。,yn+1的扰动,Taylor公式,数值积分法,RungeKutta方法,局部截断误差,整体截断误差,收敛与稳定性,4.6 多步法,某一步解和前若干步解的值有关,m步线性多步法的一般形式,其中a0,am-1,b0,bm为和h无关的常数初值为y0,y1,ym-1若bm0，方法是隐式的(implicit)否则是显式的(explicit),(4.7),AdamsBashforth四阶方法,AdamsMoulton四阶方法,多步法的截断误差,4.7 常微分方程组和高阶微分方程的数值解法,形如,欧拉方法的计算公式,隐式欧拉方法的计算公式,梯形公式,Adams公式,R-K方法,对于高阶微分方程，总可以化成方程组的形式，例如,可以化成,例,令y(x)=z，则有,用4阶RK方法，有,4.8 常微分方程边值问题的数值解法,形如,第一类边界条件,(4.8),第二类边界条件,第三类边界条件,定理 假定问题（4.8）式中的函数f在集合,上连续，且偏导数fy和fy也在D上连续。若有,则边值问题有唯一解。,例 边值问题,有,因此边值问题有唯一解。,定理 若线性边值问题,满足,则边值问题有唯一解。,打靶法（Shooting Method）,其中，y1(x)是初值问题,的解,的解，y2(x)是初值问题,有限差分法(Finite-Difference Methods),以差商代替导数，将微分方程离散成为差分方程,用Taylor展开方法,两式相加得到,由中值定理得,线性边值问题可写成,然后可得到有限差分公式,上式可写成,写成矩阵形式,其中,