第六章 常微分方程初值问题初步1

10
3、中心差商公式
y(x n 1 ) hy(x n 1 )y'(x n 1 ) h 2y''(n)
y n 1y n 1 h(x fn 1 ,y n 1 )
是多步，2阶格式，该格式不稳定
11
4、梯形公式
对微分方程 y'd yf(x,y) , x [a,b] dx
要计算出解函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值 y iy (x i) (i 1 ,..,n .)
通常取节点间距 h i x i 1 x i (i 0 ,.,.n . 1 )为步长，通常采用等 距节点，即取 hi = h (常数)。它适合计算机求解,应用广泛,具 有应用价值。
y ( x i- 1 ) y ( x i- h ) y ( x i) h y ( x i) h 2 2 !y ( x i) ( - n h ! ) n y ( n + 1 ) ( x i)
7
6.1.2 Euler公式
利用等距分割，数值微分来代替导数项，建立差分格式。
y ( x i- 1 ) y ( x i-h ) y ( x i) h y ( x i) h 2 2 !y ( x i) ( - n h ! ) n y ( n ) ( x i)
y ( x i+ 1 ) y ( x i h ) y ( x i) h y ( x i) h 2 2 !y ( x i) h n n ! y ( n + 1 ) ( x i)
y (x n 1 )y (x n) h(x fn 1 ,y (x n 1 ) )h 2 2y ''(n)

【例】 利用四阶泰勒法求解初值问题 1 y ( x y) 2 y(0) 1 其中[a , b] [0,1]。 1 y ( x y) 2 1 y (2 x y ) 4 1 y ( 2 x y ) 8 1 y (4) (2 x y ) 16 解 建立 M 文件 df.m function z=df(x,y) z=[(x-y)/2 (2-x+y)/4 (-2+x-y)/8 (2-x+y)/16]; 在 MATLAB 运行窗口输入 taylor4('df',0,1,1,10)
在Euler算法程序的最后改为 Y1=sqrt(1+2*T); T=[T' Y' Y1']; 可得到解析解和迭代法计算结果比较
6.1.3 微分方程的几何解释 在 xy 平面上，积分曲线上任一点(x,y)的切线斜率等于函数 f(x,y)的值， 如果按函数 f(x,y)在 xy 平面上建立一个向量场，则积分曲线上每一点的切线 方向与向量场在该点的方向相一致。 从初始点 P0 ( x0 , y0 ) 出发，先依向量场在该点的方向推进到 x x1 上一点 P1 ， 再从 P1 依向量场方向推进到 x x2 上一点 P2 ， 这样下去便可作出一条折线 P0 P1 P2 。设已作出折线的顶点为 Pn ，过 Pn ( xn , yn ) 依向量场的方向再推进到 Pn1 ( xn1 , yn1 ) 。 yn 1 yn f ( x n , yn ) x n 1 yn
6.1.7 梯形方法的求解 梯形方法是隐式的，用迭代求解，用欧拉方法提供 的迭代初值的迭代公式为：
(0) yn yn hf ( xn , yn ) 1 ( k 1) h (k ) y y [ f ( x n , yn ) f ( x n 1 , yn n1 n 1 ] 2 ( k 0,1,2, )

第6章 常微分方程初值问题 的数值解法
引言
在实际问题中，常需要求解微分方程(如发电机运动方程)。
只有简单的和典型的微分方程可以求出解析解，而在
实际问题中的微分方程往往无法求出解析解。
在高等数学中我们见过以下常微分方程：
f(x y ,y ) axb y ( a )y 0
f( y x ,y ,y ) a x b y ( a )y ( a ) 0, y
题（1）存在唯一解。 f(x y ,y ) axb 要求它的数值解 对于问题(1) y ( a )y 0 （通常采用等距节点）
a x x x x b 0 1 2 n
就是求未知函数 y (x ) 在区间 [ a ,b ] 上的一系列离散
上函数值 y ( x ) 的近似值 y ( k 1 , 2 , , n ) k k
2、化导数为差商的求解方法思路： 若在点 x k处的导数用差商来近似代替，如向前差商
y x y x k 1 k y x k h 则微分方程初值问题化为
'
y x y x ' k 1 k k y x k 0 , 1 , , n 1 h y a y 0 将近似号改为等号，精确解 yxk 改为近似解 y k ，得
x1
一般形式
y y h f ( x , y ) ( k 0 , ... , n 1 ) k 1 k k k
欧拉公式几何意义 用一条通过初始点的折线 近似表示解曲线 ,亦称为欧拉 折线法 ，或称为矩形法。
局部截断误差和阶数 定义 在假设 yk = y(xk)，即第 k 步计算是精确的前提下，考 虑的截断误差 Rk = y(xk+1) yk+1 称为局部截断误差 。 定义 若某算法的局部截断误差为 O(hp+1) ，则称该算法有 p 阶精度。

解： y Euler公式： n1
0.2 xn yn hf ( xn , yn ) 1.1 yn yn
0.2 0.3 0.4 0.5
xn yn xn yn
0.1
1.1000 1.1918 1.2774 1.3582 1.4351 0.6 0.7 0.8 0.9 1.0
1.5090 1.5803 1.6498 1.7178 1.7848
改进的Euler公式：
0 .2 x n (0) y n 1 1 .1 y n y n y y 0.1 [( y 2 xn ) ( y ( 0 ) 2 xn 1 )] n n n 1 (0) n 1 2 yn y n 1
xn yn xn yn
1 Q ( x , y , h ) Q ( x , y , h) f ( x , y ) f ( x , y ) 2 f ( x h, y hf ( x , y )) f ( x h, y hf ( x , y ))
设L为f关于y的Lipschitz常数,则由上式可得
En1 y( xn1 ) yn1 y( xn1 ) y( xn ) hQ( xn , y( xn ), h) 称为单步法在点 xn1 处的局部截断误差。
定义 若某算法的局部截断误差为 E n1 O( h p 1 ) ，则称该 算法有p 阶精度。 欧拉法的局部截断误差,由Taylor展开:
e n y( x n ) yn
整体截断误差除与 x n步计算有关外,还与 有关
x n1 , , x1的计算
分析计算中的某一步,显式单步法的一般形式可写为:
yn1 yn hQ( xn , yn , h)
其中 Q( xn , yn , h) 称为增量函数。如对于Euler公式其增量函 数 Q( xn , yn , h) f ( xn , yn ) y ' ( xn ) 定义 局部截断误差:设 y( x )是初值问题 (6.1)的解,用单步法计算到 第n步没有误差,即 yn y( xn ) ,则