科学与工程计算科学计算14

3 每步计算 f 的值两次， 其截断误差为 O( h )
下面对预估—校正法进行改进，将该公式写成更一般的形式
yn 1 yn h( R1k1 R2 k 2 ) k1 f ( xn , yn ) k f ( x ah, y hbk ) n n 1 2
(3)
选择常数的原则是在假设 yn y( xn ) 的前提下，使 y( xn1 ) yn1 的 阶尽量高。
与二阶龙格—库塔公式的讨论方法类似，只需8个参数满足 R 1 R 2 R 3 1 a 2 b 21 a 3 b 31 b 32 1 a R a 3 R 3 2 2 2 1 2 2 a 2 R 2 a 3 R 3 3 1 a b 2 32 R 3 6 满足条件(6)的一族公式(5)统称为三阶龙格—库塔公式。


(5)
3: 预估一校正法
改进的欧拉公式在实际计算时要进行多次迭代， 计算量较大。 在实用上，对于改进的欧拉公式(5)只迭代一次，即先用欧拉公式 算出 yn1 的预估值 y n1 ， 再用改进的欧拉公式(4)进行一次迭代得到
( 0)
校正值 y n1 ，即 ( 0 ) y y hf ( x ,y ) n 1 n n n h y y f x ,y f x ,y , n 0 , 1 , 2 , n 1 n n n n 1 n 1 2
y n 1 y n k1 k1 hf ( x n , y n )
每步计算 f 的值一次，其截断误差为 O( h 2 )。
对于预估—校正法。
1 1 y n 1 y n 2 k1 2 k 2 k1 hf x n , y n k hf x h, y k n n 1 2
主要内容 常微分数值计算方法简介 偏微分方程数值计算方法 （差分法、有限元法、有限体积法） 统计方法简介 选讲内容 计算流体力学方法简介 计算电磁学方法简介 主要参考书 科学与工程计算方法 北京理工大学出版社 数值计算方法 北京理工大学出版社等 偏微分方程数值解法 清华大学出版社、华中理工大学出版社等 微分方程的数值方法(英文) springer出版社 实用多元统计方法与SAS系统 北京大学出版社 计算电磁学的数值方法 清华大学出版社
0 . 2 1.191818 0 . 3 1.277438 0 . 4 1.358213 0 . 5 1.435133 0 . 6 1.508966 0 . 7 1.580338 0 . 8 1.649783 0 . 9 1.717779 1 1.784771
近似解与准确解比较，欧拉法的结果大致只有两位有效数字， 而预估—校正法的结果则 有3位有效数字。
预估—校正公式也常写成下列形式：
1 1 yn1 yn 2 k1 2 k2 k1 hf (xn , yn ) , n 0, 1 ,2, k hf x h, y k n n 1 2
公式的局部截断误差
定义 若某种微分方程数值解公式的截断误差是 O(h 称这种方法是 k 阶方法。
x n 1
h f ( x , y ( x )) dx f x , y x f x , y x n n n 1 n 1 x n 2
得到改进的欧拉公式
h y y f x , y f x , y n 1 n n n n 1 n 1 2
f f y x x ,y x fx y x y fy fx ff y
2 2 fxx y 2 ff f f ff f fyy xy x y y
所以这种方法实际上不能用来解初值问题。 二 龙格—库塔方法(R-K方法) R-K方法不是通过求导数的方法构造近似公式，而是通过计算不同 点上的函数值，并对这些函数值作线性组合，构造近似公式，再把 近似公式与解的泰勒展开式进行比较，使前面的若干项相同，从而 使近似公式达到一定的阶数。先分析欧拉法与预估—校正法。 对于欧拉法
其中 R1 , R2 , a, b 为待定常数。
如果每步计算三次 f 的值，可将公式写 成下列形式：
y n 1 y n R1 k1 R2 k 2 R3 k 3 k hf x , y 1 n n k 2 hf x n a 2 h, y n b21 k1 k 3 hf x n a3 h, y n b31 k1 , b32 k 2
科学与工程计算方法
常微分方程的初值问题 第一节 几种经典方法
求解常微分方程初值问题
dy f ( x, y ) dx y(x0 ) y0
的数值解，就是寻求准确解
y( x)
在一系列离散节点
x x x x 0 1 2 n
上的近似值 y , y , y , , y , 0 1 2 n
y n1 h2 h k k yxn hy xn y xn 2! k! (2)
则
yxn1 yn1 O h k 1

即公式(2)为k阶方法。
从理论上讲，只要解 y( x) 有任意阶导数，泰勒展开方法就可以构造 任意阶求 yn1 的公式，但由于计算这些导数是非常复杂的。
称为改进的欧拉法，也称为梯形法。
向前欧拉公式是关于 yn 1 的显式，即只要已知 yn , 计算便可得
yn 1 的值，而改进（向后）的欧拉公式是以 yn 1 的隐式方程给出，
不能直接得到 yn 1 。隐式方程通常用迭代法求解，而迭代过程的 实质是逐步显式化。 0 先用欧拉公式
y y h x f , y n 1 n n n
5
的四阶龙格—库塔公式是常用的公式，每步都要
计算四次 f 的值。它的一般形式是
4 h R y n 1 y n iK i i 1 K 1 f (x n, y n) i 1 K x a ,y h b K i f( n ih n ij j) j 1
n 1
k 1
)，则
假设前一步的结果是准确的， 即在 yn y( xn ) 的前提下， 考虑用 yn1 作为 y(x )的近似值而产生的截断误差，称为局部截断误差。
1． 准确解在 xn 附近展开成泰勒公式
y(x n 1 ) y(x n ) hy (x n )
y n hf (x n ,y n )
(5)
一个比较简单的三阶龙格—龙塔公式是
截断误差为 Oh
1 4 1 y n 1 y n 6 R1 6 R 2 6 R3 k hf x , y 1 n n 1 1 k 2 hf x n h, y n k1 2 2 k 3 hf x n h, y n k1 2 k 2
n
0 1 2 3 4 5 6 7 8 9 10
xn
0 0 .1
yn欧拉法Fra bibliotek1 1.1预估-矫正法
1 1.095909 1.184097 1.266201 1.343360 1.416402 1.485956 1.552514 1.616475 1.678166 1.737867
准确解
1 1.095445 1.183216 1.264911 1.341641 1.414214 1.483240 1.549193 1.612452 1.673320 1.732051
(6)
(i 2 ,3 ,4 )
(6)式中13个待定常数需满足下列11
R a a a
1 2 3 4
R b b b
四 显式龙格—库塔法
前面讨论的欧拉法与改进的欧拉法都是一步法，即计算 yn1 时，只 用到前一步值。龙格—库塔(Runge-Kutta)法(简称为 R-K 方法)是 一类高精度的一步法，这类方法与泰勒级数法有着密切的关系。 dy 一 泰勒级数法 f ( x, y ) dx 设有初值问题 y(x0 ) y0 2 h 由泰勒展开式 y x y x h y x h y x y x n 1 n n n n 2 ! k h k k 1 y x O h n k !
给出
yn 1 的迭代初值，然后再用改进的欧拉公式(4)
(0) yn 1 y n hf ( x n , y n ) h ( k 1) (k ) y y f ( x , y ) f ( x , y n 1 n n n n 1 n 1 ) 2 k 0,1,2,
yn 称为问题的数值解，数值解所满足的离散方程统
称为差分格式。 hi xi xi 1 称为步长，实用中常取定 步长。 只有初值问题的解存在且唯一时，使用数值解法才有意义， 这一前提条件由下 面定理保证。
1: 欧拉法(欧拉折线法)
hf (x y n 1 y n n, y n) x0 , n0 , 1 ,2 0 y y
取 h 的线性部分，并用
yn1 yn hf xn , yn
y n 作为 yxn 的近似值，得
dy f x, y 2. 对方程 dx 两边从
xn 到 xn1 积分，得
xn 1 xn
y( xn1 ) yxn
f ( x, y( x))dx
2:
改进的欧拉法(梯形法) 欧拉法形式简单，计算方便，但曲线的曲率较大时，欧拉法 的效果差。于是，对欧拉法进行改进，用梯形公式计算
2x dy y y dx y (0) 1 正法求初值问题
解
该方程为贝努利方程，其精确解为 y 1 2x 。
2 xn y n 1 y n h y n y n 欧拉公式的具体形式为