数值求积公式及代数精度

又称为修正的Euler公式
yn+1= yn+ 0.5h[ k1+ k2]
k1=f(xn,yn), k2=f(xn+h, yn+hk1) 局部截断误差概念
设 yn= y(xn), 称Rn+1=y(xn+1) - yn+1为局部截断误差 常表示为: O(hp+1)， p 称为单步法的精度阶数
9/18
二阶Range-Kutta公式一般形式
y1(x0) y0, y2(x0) y0
y (3 )f(x ,y ,y ,y )
y (x 0 )y 0 ,
y (x 0 )y 0 ,1 y (x 0 )y 0 2
令 y1=y, y2=y’ y3=y”
y1’=y2 y2’=y3 y3’=f(x, y1, y2, y2)
12/18
Ex1.推导矩形求积公式
b
bf(n 1 )()
R [f] a [f(x ) L n (x )d ] x a( n 1 )!n 1 (x ) dx
等距结点插值型求积公式称为Newton-Cotes公式, 偶数阶Newton-Cotes公式至少有(n+1)阶代数精度
2/18
1.梯形公式
b f ( x ) d b x a [ f ( a ) f ( b ) ( ] b a ) 3 f ()
y’= f (x) , y(a)=0
结果有: y (b ) h 6 N n 0 1 [f(x n ) 4 f(x n 1 /2 ) f(x n 1 )]
其中, h = (b – a )/N, xn= a + n h ( n = 0,1,2,···, N)
15/18
Ex5.将线性常系数非齐次高阶常微分方程初值问题:
b f( x ) d ( x b a )f( a ) f() ( b a ) 2
a
2
思考:
b f( x ) d ( x b a )f( b ) f( n ) ( b a ) 2
a
2
b f ( x ) d ( b x a ) f ( a b ) f ( n ) ( b a ) 3
yn+1= yn+ h[c1k1+ c2k2]
k1=f (xn, yn) k2=f(xn+λ2h, yn+μ21hk1)
( n=0,1,······ )
四阶龙格-库塔公式(称为经典公式)
yn+1= yn+ h[k1+2k2+2k3+k4]/6
k1=f(xn,yn), k2=f(xn+0.5h, yn+0.5hk1) k3=f(xn+0.5h, yn+0.5hk2), k4=f(xn+h, yn+hk3)
局部截断误差为 O(h5)
( n=0,1,······ )
10/18
初值问题 Euler公式 :
dyf(x,y)
dy (xx0) y0
xx0
y n 1 y n h f ( x ,y n )
y1,n1y1,nh1f(x,y1,n, ,ym,n)
即
ym,n1ym,nhm f(x,y1, n, ,ym,n)
b f( x ) d ( x b a )f( a ) f() ( b a ) 2
a
2
令
u
F(u) f(x)dx
F ( a ) 0 , F ( a ) f ( a ) F ( , ) f ( )
a
F(u)= F(a) + (u-a)F’(a) +0.5(u-a)2F ”()
f(x ) G m ( h ) O ( h 2 (m 1 ))
思考：一阶中心差商和二阶中心差商的外推公式？
7/18
常微分方程初值问题 1. Euler方法
yf(x,y) xx0 y(x0)y0
y y0 n 1 y(y x n 0 )h ,x n ( x f1n ,y x n n ), (n h 0 ,1 ,2 , )
《数值分析》习题课 IV
数值求积公式及代数精度 数值求导方法与截断误差 一阶常微分方程数值法 局部截断误差与精度
插值型求积公式:
b
n
f(x)d x
a
A j f(xj)
n
j0
拉格朗日插值 f(x) lj(x)f(xj)
令 求积余项
j0
A ja b lj(x )d,x (j 0 ,1 ,2 , ,n )
a
2
12
2. Simpson公式
b
b a
a b
af( x ) d x 6[f( a ) 4 f(2) f( b )]
复合梯形求积公式 令h=(b-a)/n
b
h
n 1
f(x )d x[f(a )f(b ) 2f(aj)h ]
a
2
j 1
误差 R [f](b 1 n a 2 2 )3f(), (a,b)
与任意的不超过n次的多项式P(x) 正交,则 wn+1(x)的所 有零点x0, x1 ,······, xn 是Gauss点
正交多项式:
1.勒让德多项式: p0=1, p1=x, p2 =0.5(3 x2 – 1 ) 2.切比雪夫多项式: T0=1, T1=x, T2=2x2 – 1 ,········
5/18
一阶向前差商
f(a )f(a h )f(a ) O (h ) h
一阶向后差商
f(a )f(a ) f(a h ) O (h ) h
一阶中心差商
f( a ) f( a h ) f( a h ) O (h 2 ) 2 h
二阶中心差商 f ( a ) f ( a h ) 2 f h ( 2 a ) f ( a h ) O ( h 2 )
二阶Range-Kutta公式
y n 1 y n 0 . 5 h [ k 1 k 2 ]
k 1f(xn,yn)
k2f(xnh,ynhk1)
11/18
高阶常微分方程初值问题
yf(x,y,y) y(x0)y0, y(x0)y0
令 y1=y, y2=y’
一阶常微分
初值条件:
方程组 yy12yf2(x,y1,y2)
6/18
外推算法 G (h )1[f(x h )f(x h )]
2 h
G ( h ) f ( x ) 1 h 2 2 h 4
G ( h / 2 ) f ( x ) 1 h 2 / 4 2 h 4 / 1 6
4G(h/2)G(h) 3
f(x ) 32 h 4/4
Gm(h)4mGm14 (h 2 m) 1Gm1(h)
2. 梯形公式:
h y n 1 y n 2 [f(x n ,y n ) f(x n 1 ,y n 1 )]
8/18
预测-校正公式
~ y n 1 y n h ( x n , fy n )
y n 1 y n h 2 [f(x n ,y n ) f(x n 1 ,~ y n 1 )]
j0
4/18
Gauss点 如果求积结点x0, x1,···,xn,使插值型求积公
式
1
n
f(x)dx
1
A kf(xk)
k0
代数精度为2n+1,则称该求积公式为Gauss型求积公 式. 称这些求积结点为Gauss点.
定理7.2 如果多项式wn+1(x)=(x – x0) (x – x1)···(x – xn)
a
2 24
13/18
Ex2.利用复合梯形公式计算积分
1s inx
I 0
dx x
使其截断误差不超过 0.5×10-3,应算多少次函数值？
提示：
six n 1
f(x) x 0coxs)d t(t
思考: 给定积分 3ex sinxdx 1
当要求误差小于10-3时用复合梯形公式和Simpson公 式计算时, 需要计算多少次函数值？
3/18
高斯型数值求积公式
1f(x )d xf(1)f(1)
1
33
b f ( x ) d b x a [ f ( b a b a ) f ( b a b a )]
a
2 2 32 2 32
x(t)batba 22
t∈[-1, 1]
思考:三点求积公式
b
n
f(x)d x
a
A j f(xj)
y(n) + a1 y(n-1) + a2 y(n-2) +·······+ an y = f( x, y, ····, y(n-1)) y(x0)=y00, y’(x0)=y01, y”(x0)=y02,···y(n-1)(x0)=y0,n-1 转化为一阶线性常微分方程组问题,并成出矩阵形式
解: 令
y1(x)=y(x), y2(x)=y’(x), y3(x)=y”(x),····, yn(x)=y(n-1)(x),
y1 y2
y2y3
yn(any1 a1yn)f(x,y1, ,yn)
16/18
Fra Baidu bibliotek
Ex6. 初值问题
y axb y(0) 0
有解y(x)=0.5a x2 + b x 。若取 xn = nh，yn为欧
写出梯形公式求解的计算格式，取步长h = 0.2,计 算函数 y(x) 在 x = 1.6 处的近似值 y2.
Ex9.写出用修改的欧拉法求问题
y 2 y 2 yexs ixn 0x3 y(0 )0 ,y(0 ) 1 /2 的计算公式
18/18
14/18
Ex3. 定积分
b
f (x)dx
的计算问题可化为初值问题
a
y’= f (t) , y(a)=0
N1
试证明用Euler公式计算结果为 y(b) f(tn)h
其中, h = (b – a )/N,
tn= a + n h
n0
( n = 0,1,2,···, N)
Ex4. 试证明4阶Range-Kutta公式解[a, b]内初值问题
拉方法得到的数值解，试证明
y(xn) – yn = 0.5 a h xn
Ex7. 初值问题
y axb y(0) 0
若取 xn = nh，yn为用梯形公式计算所得的数值解，
记y(xn)为初值问题的在x=xn处的解析解。试证明:
y(xn) = yn
17/18
Ex8. 初值问题
y83y
y(1)2
(1x2)