计算方法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.1 y2=y1+hf(x1,y1)=1.1+0.1(1.1-20.1/1.1)=1.191818 y3=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 yc=1+0.1X1.12=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/2 k1=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+c3k3 k1=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 及N Step 2: (b-a)/N=h,a=x,y0=y Step 3: 输出 (x,y) Step 4: For I=1 T0 N hf(x,y)=k1 hf(x+h/2,y+ k1/2)= k2 hf(x+h/2,y+k2/2)=k3 hf(x+h,yk3)=k4 y+(k1+2k2+2k3+k4)/6=y x+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法，标准四阶龙格库法也收敛。,4.5.2 稳定性,在讨论收敛性时，一般认可：数值方法本身计算过程是准确的。实际上，并非如此： 初始值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时恒成立 梯形公式恒稳定。,y