《计算方法_常微分方程数值解法课件》由会员分享，可在线阅读，更多相关《计算方法_常微分方程数值解法课件（82页珍藏版）》请在金锄头文库上搜索。
1、常微分方程的数值解法 Numerical Solutions to Ordinary Differential Equations,对象,一阶常微分方程初值问题:,一阶常微分方程组初值问题:,高阶常微分方程初值问题:,(4.1),对象,一阶常微分方程初值问题:,f(x,y)在a,b上连续,且满足Lipschitz条件:,实际工程技术、生产、科研上会出现大量的微分方程问题很难得到其解析解，有的甚至无法用解析表达式来表示，因此只能依赖于数值方法去获得微分方程的数值解。,求函数y=y(x)满足： y(x) =f(x, y(x),则初值问题Problem I有唯一解,微分方程的数值解法： 不求y=y(
2、x)的精确表达式,而求离散点x0,x1,xn处的函数值 设Problem I的解y(x)的存在区间是a,b，初始点x0=a，取a,b内的一系列节点x0, x1,.,xn。a= x0 x1 xn =b，一般采用等距步长。,用数值方法，求得y(x)在每个节点xk 的值y(xk ) 的近似值，用yk 表示，即yk y(xk)，这样y0, y1,.,yn称为微分方程的数值解。,求y(x)求y0, y1,.,yn,？,计算过程：,方法：采用步进式和递推法,将a,bn等分， a= x0 x1 xn =b，步长h= ，xk=a+kh,怎样建立递推公式？ Taylor公式 数值积分法,思想:用向前差商近似代替
3、微商.,欧拉公式（Euler Scheme）,4.1 欧拉公式,几何意义,y(x)过点P0(x0,y0)且在任意点(x,y)的切线斜率为f(x,y),2. y(x)在点P0(x0,y0)的切线方程为：,3. 类似2，过P1以f (x1,y1)为斜率作y(x)的切线，在其上取点 P2(x2,y2)，依此类推,4.折线P0 P1 P2 Pn作为曲线y(x)的近似,欧拉折线法,在切线上取点P1 (x1,y1),y1正是Euler 公式所求。,用向后差商近似代替微商：,隐式欧拉公式,欧拉法（续）,here,注：用隐式欧拉法，每一步都需解方程（或先解出yn+1的显式表达式），但其稳定性好。,用中心差商近
4、似代替微商：, 二步欧拉法,注：计算时，先用欧拉法求出y1 ，以后再用二步欧拉法计算。,欧拉法（续）,公式,单步否,显式否,单步,显式,单步,隐式,二步,显式,截断误差y(xn+1)-yn+1,局部截断误差,定义1 假设yn=y(xn) ,即第n步计算是精确的前提下,称Rn+1=y(xn+1)-yn+1为欧拉法的局部截断误差.,定理 欧拉法的精度是一阶。,定义2 若某算法的局部截断误差为O(hp+1),称该算法有p阶精度.,注：无yn=y(xn) 前提下，称Rn+1为整体截断误差。,？,分析：证明其局部截断误差为O(h2),可通过Taylor展开式分析。,证明： Euler 公式为,令yn=y
5、(xn)，下证： y(xn+1)-yn+1 = O(h2),由 y(x) =f(x, y(x),二步欧拉法的局部截断误差,定义3 假设yn=y(xn) , yn1=y(xn1),称Rn+1=y(xn+1)-yn+1为 二步欧拉法的局部截断误差.,定理 隐式欧拉法的精度是一阶，二步欧拉法的精度是二阶。,证明： 对二步欧拉法进行证明，考虑其局部截断误差，,令yn=y(xn) , yn1=y(xn1),将上两式左右两端同时相减：,二步欧拉法的局部截断误差为O(h3)，其精度是二阶。,例: 求 的近似值。,解:这儿 ，,由欧拉公式 得:,又其精确解为,整体误差 ，下面对其加以分析,从表中看出误差在逐步
6、增加、积累,局部截断误差,而误差是,数值积分法,对右端的定积分用数值积分公式求近似值：,（1）用左矩形数值积分公式：,（2）用梯形公式：,梯形公式, 梯形公式:将显示欧拉公式,隐式欧拉公式平均可得, 梯形公式是隐式、单步公式,其精度为二阶,梯形公式的精度,证：令yn=y(xn),由Talor公式有,分析：梯形公式是隐式公式，证明其局部截断误差为O(h3) 要用到 二元函数的Taylor公式。,定理：梯形公式 的精度是2阶的.,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) ， (
7、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),从而y(xn+1)=yn+1+h fy(xn+1,)y(xn+1)-yn+
8、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)/
9、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),梯形公式的应用,例4.1 用梯形公式求初值问题的解在x=0.01上的值y(0.01).,解: 取h=0.01 x0=0 y0=y(0)=1 则 y(0.01)y1 f(x,y)=y，由梯形公式,基于稳定性考虑,解析解,欧拉公式的比较,4.2 改进的Euler法,Euler公式 yn+1=yn+hf(xn,yn) 显式 一阶 梯形公式 yn+1=yn+hf(xn,yn)+f(xn+1,yn+1)/2 隐式 二阶,Euler公式 计算量
10、小，精度低。 梯形公式 计算量大，进度高。,综合两个公式，提出预报校正公式：,预报,校正,改进的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)/2
11、1.1118,同样可求y2,y3,y4 ，见P95表: y(xn)-yn随n增大而增大，表明误差积累。,注： （1）令y(xn)=yn，可推导改进的Euler法的局部截断误差 为O(h3),具有二阶精度。,（2） 改进的Euler法也可写成如下平均化形式,HW: p.117 #3，4,龙格库塔方法,改进的Euler公式 :,Euler公式：yn+1=yn+hf(xn,yn),写成,精度：一阶,精度：二阶,由Lagrange中值定理，,称为y(x)在xn,xn+1上的平均斜率,取,Euler公式,取,改进Euler公式,Euler公式用一个点的值k1作为k*的近似值，而改进的Euler公式用二个
12、点的值k1和k2的平均值作为k*近似值，改进的Euler法比Euler法精度高。,龙格库塔法,Runge-Kutta法的思想：在xn,xn+1内多预报几个点的ki值并用其加权平均值作为k*近似值。而构造出具有更高精度的计算公式。,4.3.2 二阶龙格库塔方法,以k1与k2的加权平均来近似k*，设,其中w1,w2,为待定参数。适当选取参数，使(*)式的精度为二阶，即使其局部截断误差为O(h3),令y(xn)=yn，由泰勒公式：,二阶龙格库塔法,由多元函数的泰勒公式,则有,注：上述方程组有四个未知量，只有三个方程，有无穷多组解。,比较与要使 y(xn+1)-yn+1=O(h3),二阶龙格库塔法,取
13、任意一组解便得一种二阶龙库公式。,当w1=w2=1/2, =1时二阶Runge-Kutta公式为,yn+1=yn+k1/2+k2/2 k1=hf(xn,yn) k2=hf(xn+h,yn+k1),此即改进得Euler法,取w2=0 ，w2=1，=1/2，=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
14、=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时：计算量太大，精确度不一定提高，有时会降低。,龙格库塔法,对于经典的四阶Runge-Kutta法给出如下算法：,Step 1: 输入a,b,y0 及N Step 2: (b-a
15、)/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) END,龙格库塔法,例：用四阶经典RungeKutta方法解初值问题：,(1)求 ，,龙格库塔法,(2)求 ，,变步长龙格库塔法,用户提出问题I ：,问题：如何判断|y(xn)-yn| 精度值y(xn)未知。 ：如何取h=? 解：如用p阶龙格库塔法计算，局部截断误差为O(hp+1),如
16、 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（） 将步长折半反复计算，直至为止，最后一次运算的前一次
