表7-6
j bj 0 1 1 -1/2 2 -1/12 3 -1/24 当?k=0?时,(7)式为 ? yn?1?yn?hfn?1 当?k=1?时,(7)式为梯形公式 当 k=3?时,(7)式为? ?即?
yn?1?yn?h24(9fn?1?19fn?5fn?1?fn?2?1)
yn?1?yn?h(fn?1?12?fn?1?112?fn?1?2124?fn?1)3
(8)
类似于亚当斯显式公式,可以导出隐式公式(8.3.8)的局部截断误差为??
y(xn?1)?yn?1??19720hy5(5)(xn)???O(h)5
?所以隐式公式(8)是四阶方法。? 7.3.3亚当斯预估一校正公式
? 把亚当斯显式公式与隐式公式联立使用,构造亚当斯预估一校正公式。以四阶亚当斯方法为例有下列公式:?
yn?1?yn?(0)h24??
(55fn?59fn?1?37fn?2?9fn?3)h24
(9)
yn?1(k?1)?yn?(9f(xn?1,yn?1?19fn?5fn?1?fn?2)(k)迭代过程进行到连续 两次迭代结果之差的绝对值小于给定的精度ε为止。
(k?1)这时取yn?1?yn?1,然后转入下一步计算。?只迭代1次的公式
yn?1?yn?(0)h24h(55fn?59fn?1?37fn?2?9fn?3)
(.10)
yn?1(k?1)?yn?24(9f(xn?1,yn?1?19fn?5fn?1?fn?2)(k)?称为亚当斯预估一校正公式。第1式称为预估公式,第2式称为校正公式。其局部截断误差为??
16
y(xn?1)?yn?1??19720hy5(5)(xn)?O(h)5
下面我们具体讨论误差估计问题。设 由(10)式求得准确解 y(xn+1)的近似解为yn+1,由于??
y(xn?1)?yn?1?(0)25172019720hy5(5)(xn)
y(xn?1)?yn?1??hy5(5)(xn)上面二式相减得??
y(xn?1)?yn?1?(0)270720? 所以??
hy5(5)(xn)
y(xn?1)?yn?1yn?1?yn?1(0)??19270
??于是有?
y(xn?1)?yn?1??19270(yn?1?yn?1)(0)
??若??
?19270(yn?1?yn?1?ε
(0)
则可持续求解yn+2,?,否则,应缩小步长h重新计算。?
与同阶的龙格—库塔方法相比较,亚当斯方法计算量小,公式简单,程序易于实现。但是由于亚当斯方法在计算yn+1时,不仅用到前面一步的信息,而且还要用到更前面几步的信息,因此它不能自动开始。开始的前几个值必须借助于一步法获得。?
例:分别用四阶亚当斯显式公式与预估一校正公式求解初值问题??
dx dy?y?2xy
y(0)?1 0?x?1,h?0.1
解:用第二节例1,按标准四阶龙格—库塔方法求得的结果y1,y2,y3作为开始值,然后用显式公式(6)与预估一校正方法(10) 进行计算,计算结果如表7-7所示。
17
表7-7
xn 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
yn R-K方法 0 1.09544553 1.18321675 1.26491223 显式方法 1.34155176 1.41404642 1.48301819 1.54891888 1.61211634 1.67291704 1.73156976 预估计-校正法 准确解 1.341641136 1.41421383 1.48323983 1.54919338 1.61245154 1.67332000 1.73205072 1 1.09544512 1.18321596 1.26491106 1.34164079 1.41421356 1.48323970 1.54919334 1.61235155 1.67332005 1.73205081 近似解与准确解进行比较,显式方法的结果有4位有效数字,而预估一校正法的结果则有7位有效数字。
第四节 线性多步法
线性多步法的一般公式为
yn?1?0yn?1yn?1=+
+?+?kyn?k+h(??1fn?1??0fn?...??kfn?k)
(1)
当β-1 =0时,公式(1)是显式的;当β-1 ≠0时 ,公式(1)是隐式的。“线性”是指?yn?1是yn,yn?1, ? , yn?k, fn?1,fn,?,fn?k的线性组
合,“多步”是指计算yn+1时 ,不仅用到前一步的结果yn,而且用到更前面几步的结果。?所有形如(8.4.1)式的公式都 可以利用泰勒展开的方法构造出来。作为例子,我们推导形如??
yn?1?0yn?1yn?1?2yn?2=
+
+
+h(
?0fn??1fn?1??2fn?2??3fn?3)
(2)
的数值计算公式。?
18
假设 yi=y(xi), i=n, n-1, n-2 fi=f(xi,yi)?i=n,n-1,n-2,n-3 将上述各项在?xn?处泰勒展开 Yn-i=y(xn-ih)=y(xn)-ih
fn?iy(xn)'+
(ih)2!2y(xn)-
''(ih)3!3y(xn)+?, i=1,2
'''=
y(xn?i)y(xn)''=
y(xn?ih)'
2 =-
ihy?(xn)+
(ih)2!y???(xn)-
(ih)3!3y(4)(xn)+?, i=1,2,3
?代入(2)式 得?
yn?1=(?0??1??2)y(xn)+(??0?2?2??0??1??2??3)hy?(xn)+
(?1?4?2?2?1?4?2?6?3)(
h22!y??(xn)+ h3??1?8?2?3?1?12?2?27?3)
3!y???(xn)+
4(
?1?16?2?4?1?32?2?108?3)
h4!y(4)(xn)+... (3)
将?y(xn?1)在xn处泰勒展开??
y(xn?1)?y(xn)?hy?(xn)?h22!y??(xn)?h33!y???(xn)?h33!y???(xn)?... ?(4)
比较(3)式与(4)式,让h的同次幂的系数相等,得?? ? ?0??1??2?1
??1?2?2??0??1??2??3?1
?1?4?2?2?1?4?2?6?3?1 (5) ??1?8?2?3?1?12?2?27?3?1 ?1?16?2?4?1?32?2?108?3?1 ??
如果(5)的前7个方程求出未知数?0,?1,?2,?0,?1,?2,?3,则其局部截断误差为?
O(h)?.当然,也可以从前k
7个方程解出未知数(k≤7),其局部截断误差为
19
O(h).比如从前5 个方程解出
k7个未知数,因此,有2个自由未知数,若令
?1??2?0,则从(5)的前5 个方程解得?? ?
?0?1,?0?5524,?1??5924,?2?3724,?3??924?
?对应的公式是??
yn?1?yn?h24(55fn?59fn?1?37fn?2?9fn?3)?
它正好是四阶阿达姆斯显式公式。?
类似地, 若设计算公式为?
yn?1= ?0yn+?1yn?1+?2yn?2+h(??1fn?1??0fn??1fn?1??2fn?2)?
为使这类公式为四阶方法,系数应满足??? ?0??1??2?1
??1?2?2???1??0??1??2?1
?1?4?2?2??1?2?1?4?2?1 (7) ??1?8?2?3??1?8?1?12?2?1 ?1?16?2?4??1?4?1?32?2?1 特别地,取?1??2?0,则从(7)式解得?? ?
?0?1,??1?924,?0??h241924,?1??524,?2?124?
对应的公式是
yn?1?yn?(9fn?1?19fn?5fn?1?fn?2)???
它正好是四 阶亚当斯隐式公式。
第五节 方程组与高阶方程的数值解法
?前面几节介绍了一阶微分方程初值问题的各种数值解法,这些解法同样适用 于一阶微分方程组与高阶方程。为了避免书写的复杂,下面仅就两个未知函数的方程组和二阶方程为例说明各种方法的计算公式。 7.5.1一阶方程组
?设有一阶微分方程组初值问题
20