07常微分方程初值问题的数值解法

1、第九章常微分方程初值问题的数值解法第一部分内容提要一、数值解的一般概念常微分方程初值问题y(x)二f (x, y)的数值解是指通过一定的近似方法ly(xo) =y得出准确解y = y(x)在一列离散点x0,x1, X2,|( ,Xn,()|上的近似值y, yi,y2,|i(, yn,i。数值解的特征是步进式，即y(x)在1点的近似值yn.i是 由Xn,Xn,|l(等若干点处的近似值yn,yn,川的信息给出的递推公式。若yn -1依 赖于前面k步的值yn, yn二,|1(, yn丄彳，贝U称为k步法；k = 1称为单步法。利用y(x)在Xn,XndH,Xni的精确解y(Xn), y(Xn),|I

2、 ( , y(Xn借助某种 算法计算出yni，则称y(Xnyn 1为该方法的局部截断误差。如果一个算 法的局部截断误差是0(hp1),则称该方法是p阶的；而利用数值解yn,yndyn_i得到的yni与微分方程的精确解之差y(xn-yni称为整体 截断误差，即是该数值方法的误差。对于固定的X Xo,取h =XX0,用某种算法得到yn,如有nim y(x)-yn=0,则称该方法是收敛的。注意，因x是固定的，随着hT 0, 数值解的步数n：:。在实际计算时由于舍入误差不可避免，实际得到数值解是 扎,稳定性即研究yn-yn,是否随着计算步骤n的增加而增加。 通常所提的稳定性是通过模 型方程y、 y(：

3、 ： 0)来讨论的。 若当某一步yn有舍入误差时， 在以后的计算 中误差不会逐步扩大，则称这种稳定性为绝对稳定性。二、简单单步法及其收敛性、稳定性Euler法ynYnhf (Xn,yn)的局部截断误差为O(h2)，整体截断误差为2O(h)，即一阶收敛。对于模型问题y = y(:：0)，当0 ”：h：- 时，Euler扎法是数值稳定的。隐式Euler法yn i=yn hf (Xn i,yn.J的误差与Euler法相同，但是无条件 稳定：即对任意步长h 0，隐式Euler法都是稳定的。梯形法ynynh21 f (Xn, yn) f (x. 1, y. 1) l的误差比Euler法高一阶，也 是无条

4、件稳定的。改进Euler法yn 1二yn三卄风川)f (Xn 1, ynhf (Xn,yn)l是一种预测 校正方法：y訂=yn+hf (Xn, yn)-Euler法预测yc+ =y- f (XnynX f x %卩廿梯形法校正2 -它保持了梯形法的误差阶数，但不是无条件稳定的。三、龙格-库塔方法龙格-库塔类算法采用区间l-Xn,Xni1内若干点的斜率的加权平均来近似整 个区间的平均斜率，一般形式为fsYn厂yn 2KiI匚K厂心讨n)i_1Ki二f (XnCih ,ynCih ajKj小121, sIj=1如经典的4级(s=4)4阶(局部截断误差为0(h5)Runge-Kutta公式为hyn

5、1二ynT(K 0取极限得:2h这证明了数值解的收敛性。求解微分方程初值问题y、-5y，y(0) =1时，试推出其绝对稳定区间。分析：将格式应用于所给方程可导出 误差传播方程，从而求出绝对稳定区间 解：记舍入误差为；n。该隐式单步法应用于方程y=-5y时75hyn50h -15十-=-3 5h 3 5hk1-5 ynk1 =k1-15yn3 5hyn 1 = yn45 +50h -15I 3 + 5h25h2-20h66 10hyn从而误差传播方程为；n 1225h - 20h66 10h解不等式25h2-20h+66 10hyny(xn)、用隐式单步法(该方法也属于隐式Ru nge-Kutt

6、a方法)h(3k1k2)4h依=f (x +-, ynynyn3即1)k2二f (Xnh, ynhki)k -5 ynhb二-5yn故，该方法应用到所给方程的绝对稳定区间为(O,605三、用Taylor展开原理构造形如yn勺=(ynyn _!)h( _0fn 1fn J)的两步法，试确定系数：， +,使方法具有二阶精度，并推导其局部截断误差主 项。分析：本题考察构造多步法的方法；二阶精度即局部截断误差为0(h3)。解：由多步法局部截断误差的定义Tn 1= y(Xn1)一： y(Xn)- ：y( X. _1)- h。y(X.)一“(X. _i)L将y(Xn 1)，y(Xnj)和y(Xn分别在点X

7、n展开得：.2.3h川h叩4y(Xn 1) = y(Xnh)二y(Xn) hy (Xn)y (Xn)y (Xn) O(h )2 6h2hy(x曲(Xn-hry(Xn)-hy(Xn)石八和-厂区)o(h4)ZXnJZXn-hXZXnLhZXnLy区)0(h3)将此三式代入局部截断误差(3)式，我们有Tn1)y(Xn) hyoa1一0- U)h2y (XngG 云y(加一塔)亦)要使方法具有二阶精度，必需1 -2： =0, 1 ：- 一：0- 0, 1_： 2* =0.1解得，=033飞，所以局部截断误差主项为8h3y(Xn)。此时1 * -3四、取 i，试用E“er法求解初值问题；(；：甞(0律0.2分析：先将二阶方程写为一阶方程组，再用Euler法求解 解：弓I进u=y,v = y，并记向量函数y= jU，则原二阶方程变为一阶方程组yj将Euler法yn .1二ynhf(Xn,yn)应用到该方程组得:yi二y。hf(Xo,y0)=+ 0.1广0、广1、即，用Euler法求得原二阶方程的数值解为yi=12=0.