数值分析方法第六章
《数值分析》第六章答案
习题61．求解初值问题y x y +=' )10(≤≤x 1)0(=y取步长2.0=h ，分别用Euler 公式与改进Euler 公式计算，并与准确解xe x y 21+-=相比较。
解： 1) 应用Euler 具体形式为 )(1i i i i y x h x y ++=+，其中i x i 2.0= 10=y 计算结果列于下表i i x i y )(i x y i i y x y -)( 1 0.2 1.200000 1.242806 0.042806 2 0.4 1.480000 1.583649 0.103649 3 0.6 1.856000 2.044238 0.188238 4 0.8 2.347200 2.651082 0.303882 5 1.0 2.976640 3.436564 0.4599242) 用改进的Euler 公式进行计算，具体形式如下： 10=y)()(1i i i D i y x h y y ++=+ )()(11)(1D i i i C i y x h y y +++++= )(21)(1)(11c i D i i y y y ++++= 4,3,2,1,0=i计算结果列表如下i i x i y )(1D i y + )(1c i y + i i y x y -)( 0 0.0 1.000000 1.200000 1.280000 0.000000 1 0.2 1.240000 1.528000 1.625600 0.002860 2 0.4 1.576800 1.972160 2.091232 0.006849 3 0.6 2.031696 2.558635 2.703303 0.012542 4 0.8 2.630669 3.316803 3.494030 0.020413 5 1.0 3.405417 0.0311473. 对初值问题1)0(=-='y y y)0(>x ，证明用梯形公式所求得的近似值为ii hh y ih y )22()(+-=≈ ),2,1,0( =i并证明当0→h 时，它收敛于准确解ix e y -=，其中ih x i =为固定点。
数值分析课后参考答案06
第六章习题解答1、设函数01(),(),,()n x x x φφφ 在[,]a b 上带权()x ρ正交，试证明{}()nj j x φ=是线性无关组。
证明：设0()nj jj l x φ==∑，两端与01()(,,,)kx k n φ= 作内积，由()jx φ的正交性可知，200(),()((),())((),())()()n n b k j j j k j k k k k k a j j x l x l x x l x x l x x dx φφφφφφρφ==⎛⎫==== ⎪⎝⎭∑∑⎰， 于是有001(,,,)k l k n == ，即{}()nj j x φ=是线性无关组。
2、试确定系数,a b 的值使22(()cos )ax b x dx π+-⎰达到最小。
解：定义02,[,]f g C π∈上的内积为20fgdx π⎰，取011(),()x x x ϕϕ==，()s x ax b =+，()cos f x x =，则法方程为0001010111(,)(,)(,)(,)(,)(,)f a f b ϕϕϕϕϕϕϕϕϕϕ⎛⎫⎛⎫⎛⎫= ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭ 其中()2000112,dx ππϕϕ=⨯=⎰，()2201018,xdx ππϕϕ=⨯=⎰，()3211024,x xdx ππϕϕ=⨯=⎰，()2001,cos f xdx πϕ==⎰，()21012,cos f x xdx ππϕ==-⎰，于是方程组为22312812824a b πππππ⎛⎫⎛⎫ ⎪⎛⎫ ⎪ ⎪= ⎪ ⎪- ⎪ ⎪⎝⎭⎝⎭ ⎪⎝⎭，解之得1158506644.,.a b ==-。
3、已知函数11()(,)f x x =∈-，试用二类Chebyshev 多项式()n U x 构造此函数的二次最佳平方逼近元。
解：法一、取20121(),(),(),x x x x x ϕϕϕ===()()()00112222235,,,,,ϕϕϕϕϕϕ===，()()()011202203,,,,ϕϕϕϕϕϕ===，同时由二类Chebyshev 多项式的性质知 ()()()11101211028,,,,,f f f x ππϕϕϕ---======⎰⎰⎰于是可得法方程为0122203220003220835c c c ππ⎛⎫⎛⎫⎪ ⎪⎛⎫ ⎪ ⎪ ⎪⎪= ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭ ⎪ ⎪⎝⎭ ⎪⎝⎭，解之得0121.0308,0,0.7363c c c ===-， 于是()f x 的二次最佳逼近元是2001122() 1.03080.7363x c c c x ϕϕϕϕ=++=-法一、二类Chebyshev 多项式2012()1,()2,()41U x U x x U x x ===-，取内积权函数()()x f x ρ==，于是11200114(,)(1)3f U fU dx x dx ρ--==-=⎰⎰，1121111(,)2(1)0f U fU dx x x dx ρ--==-=⎰⎰，112222114(,)(41)(1)15f U fU dx x x dx ρ--==--=-⎰⎰ 由()n U x 正交性及(,)2n n U U π=可得0000(,)8(,)3f U c U U π==，1111(,)0(,)f U c U U ==，2222(,)8(,)15f U c U U π==-， 于是()f x 的二次最佳逼近元为001122()x c U c U c U ϕ=++=21632515x ππ- 4、设012{(),(),()}L x L x L x 是定义于[0,)+∞上关于权函数()xx eρ-=的首项系数为1的正交多项式组，若已知01()1,()1L x L x x ==-，试求出二次多项式2()L x 。
数值分析第六章课件
a(1) 1n
x1
b(1) 1
a(1) 21
a(1) 22
a(1) 2n
x2
b(1) 2
.
a(1) m1
a(1) m2
a(1) mn
xn
b(1) m
将(2.1)记为A(1)x=b(1)，其中
a(1) 11
a(1) 12
a(1) 1n
a11
a12
a1n
A(1)
a(1) 21
5.2 高斯消去法
本节介绍高斯消去法(逐次消去法)及消去法和 矩阵三角分解之间的关系. 虽然高斯消去法是一种 古老的求解线性方程组的方法(早在公元前250年 我国就掌握了解方程组的消去法)，但由它改进、 变形得到的选主元素消去法、三角分解法仍然是目 前计算机上常用的有效方法.我们在中学学过消去 法,高斯消去法就是它的标准化的、适合在计算机 上自动计算的一种方法.
有的问题的数学模型中虽不直接表现为含线性方 程组,但它的数值解法中将问题“离散化”或“线性 化”为线性方程组.因此线性方程组的求解是数值分 析课程中最基本的内容之一.
关于线性方程组的解法一般有两大类：
1. 直接法 经过有限次的算术运算,可以求得方程组的精确解( 假定计算过程没有舍入误差).如线性代数课程中提到 的克莱姆算法就是一种直接法.但该法对高阶方程组 计算量太大,不是一种实用的算法.
下面讨论求解一般线性方程组的高斯消去法.由
a11 x1 a12 x2 a1n xn b1 a21 x1 a22 x2 a2n xn b2
am1 x1 am2 x2 amn xn bm
a(1) 11
a(1) 21
x1 x1
a(1) 12
x2
数值分析第六章_ppt课件
12 j i n n
% | D l = l ( e ) a | = e , j = 1 , 2 , L , n j| j
1 n
E 设 10,n 20 , 则 101 .若 在 A 的上三角位置，则特征值并无扰动.
湘潭大学数学与计算科学学院 上一页 下一页 8
定理3 设 A 为 n阶实对称矩阵，其特征值
为 1 2 n，则
n x 0
1 x 0
m i n R ( x ) m i n ( A x , x )
m a x R ( x ) m a x ( A x , , x )
通过求它的根来求矩阵的特征值,实际计算中并不采用。
数值方法 基本思想: 直接从矩阵A 或者对A 做一系列的 相似变换后得到的具有更简单形式的矩阵入手,
设计迭代过程;最后求得A 的近似特征值和相应
的特征向量.
湘潭大学数学与计算科学学院
上一页
下一页
4
特征值的估计及扰动问题
1、特征值的估计
D ( A ) { z c :| z a | | a | }, i 1 , 2 , ,n i ii ij i
为矩阵 A关于向量 x 的Rayleigh（雷利）商.
A 为 n 阶实对称矩阵，则其特征值皆为实数，
记做
1 2
， n
并且存在规范正交特征向量系满足：
Au u ,i 1 , 2 , , n , ( u , u ) ,i , j 1 , 2 , , n i i i i j ij
3 4
湘潭大学数学与计算科学学院
上一页
下一页
7
实对称矩阵的极大-极小定理：
数值分析Newton法
迭代法。
© 2009, Henan Polytechnic University §3 Newton法
2 2
第六章 方程求根
6.3.1 牛顿迭代法 基本思想：将非线性方程f(x)=0 线性化
f ( ) 2 f ( x ) f ( x0 ) f ( x0 )( x x0 ) ( x x ) 0 Newton 2 迭代公式 f ( x ) 0 0 f ( x ) f ( x0 ) f ( x0 )( x x0 ) x x0 f ( x0 ) f ( x0 ) x1 x0 作为第一次近似值 f ( x0 ) f ( xk ) 重复上述过程 xk 1 xk f ( xk )
x* x
x 2 x 1 x0
牛顿法也称为切线法
© 2009, Henan Polytechnic University §3 Newton法
f ( x1 ) x2 x1 f ( x1 )
4 4
第六章 方程求根
6.3.3 牛顿法的收敛性与收敛速度 （局部收敛性定理）设 f (x)C2[a, b]，若 x* 为 f (x) 在[a, b]上的根,且 f (x*) 0，则存在 x* 的 邻域 U ( x ) 使得任取初始值 x0 U ( x )， Newton 法产生的序列{ xk } 收敛到 x*，且满足
2
则由Newton迭代公式
xn c 1 c x n 1 x n ( xn ) 2 xn 2 xn
选取初始值 x0 c即可
© 2009, Henan Polytechnic University §3 Newton法
x0 c可以吗？
1212
第六章 方程求根
《数值分析》第六章
有局部收敛性.
证 明 . 由 连 续 函 数 的 性 质 ， 存 在 x* 的 邻 域
Δ : x − x* ≤ δ
，使 ∀x ∈ Δ 成立 ϕ '( x) ≤ L < 1 ，此外，
对于任意 x ∈ Δ ，总有 ϕ ( x) ∈ R ，这是因为
15 16
迭代法不一定收敛. 对同一个问题，不同的迭代法， 可能有的收敛，有的不收敛. 如下例.
Th 1 假定函数 ϕ (x) 满足： 1 对任意 x ∈ [a, b] 有， ϕ ( x) ∈ [a, b] （即，映像入内）
∀x ∈ [a, b] ， ϕ '( x ) ≤ L < 1 2 存在非负数 L < 1 使得， （压 缩映射）
k → ∞ 时成立下列渐近关系式
= xk − x * 当
求根 x * 的邻近连续，并且满足：
ϕ '( x* ) = ϕ ''( x* ) = L = ϕ ( p −1) ( x * ) = 0 ， ϕ ( p ) ( x * ) ≠ 0
ek +1 → C ( C ≠ 0) e kp
则称该迭代过程是 p 阶收敛的. 特别地， p = 1 时称为线性收敛，
*
* *
* * 假设 x , y ∈ [a, b] 是任意的两个根，因为
xk = x * . 故 lim k →∞
x* − y* = ϕ ( x* ) − ϕ ( y* ) = ϕ '(ξ )( x* − y* ) ≤ L x* − y*
* * 故 x = y , 即， x = ϕ ( x ) 在[a,b]上有唯一的根.
数值分析(李庆扬)第六章资料
（n1） B （n） g
若收敛
x x { （k）} * ，则
x* Bx* g
n 0，1，2，
即
(I B)x* g D1Ax* D1b
Ax* b
故如果序列收敛, 则收敛到解.B 称迭代矩阵.
例：用Jacobi迭代法求解 1x01x1 10xx2222xx337823 x1 x2 5x3 42
k
k
即x是方程组Ax b的解。
引入误差向量
(k 1) x(k 1) x
所以 lim x(k) x 等价于 lim (k) 0
k
k
由
x(k 1) Mx (k ) g
x Mx g
则可得
(k 1) M (k )
(k ) M (k 1) M k (0)
问题是在什么条件下
满足
x(k1) Bx(k) g (k 0,1, 2, )
此过程所给出的迭代法称为Jacobi迭代法，又称简单
迭代法。
Jacobid迭代的矩阵形式
0
B
b21
b12
0
b1n 1
b2n
0
0 1
0 1
0
b21
b 12 1
b1n b 2n
b b
n1
n2
0
0
0
1
, n).
0 b12 b13 若记 B b21 0 b23
bn1 bn2 bn3
b1n1 b1n
g1
b2n1
b2
n
g
g
2
bnn1 0
gn
则方程组可简记为 x Bx g
选初值向量x(0)代入 x(1) , x（1） Bx（0） g，代入x(1)
数值分析第六章_数值插值方法
M n1 (n 1)!
n1 ( x)
说明：
n=1时，
R1 ( x)
1 2
f
( )2 (x)
1 2
f
( )(x
x0 )(x
x1)
n=2时，
( [x0 , x1])
R2 (x)
1 6
f
( )(x
x0 )(x
x1)(x
x2 )
( [x0 , x2 ])
,
x1,
Hale Waihona Puke xn)1
x1
x12
x1n
n
( xi
ni j1
xj)
1 xn xn2 xnn
因 xi x j (i j) 故上式不为0。
据Cramer法则，方程组解存在且唯一。 故Pn (x)存在且唯一。虽然直接求解上述方程组 可求得插值多项式，但繁琐复杂，一般不用。
得关于a0,a1,…,an的n+1阶线性方程组
a0 a1x0 a0 a1x1
an x0n an x1n
y0 y1
a0 a1xn an xnn yn
其系数行列式是Vandermonde行列式
1 x0 x02 x0n
V
( x0
jk jk
(j,k=0,1)
称l0 (x)及l1 (x)为线性插值基函数。
2. 抛物插值：n=2情形
假定插值节点为x0, x1, x2 ，求二次插值多项式 L2 (x)，使 L2(xj)=yj (j=0,1,2) y= L2 (x)的几何意义就是过 (x0, y0)，(x1, y1) ， (x2, y2)三点的抛物线。 采用基函数方法，设
数值分析第六章插值法
P( x) an x an1 x
n
n1
a1 x a0
是函数 y f ( x) 在区间[a, b]上的n+1个互异的节点 x i (i=0,1,2,…,n )上的插值多项式,则求插值多项式P(x) 的问题就归结为求它的系数 a i(i=0,1,2,…,n )。
,现要求用线性函数 p( x) ax b 近似地代替f(x)。选 择参数a和b, 使 p( xi ) f ( xi )(i 0,1)。称这样的线性函数 P(x)为f(x)的线性插值函数 。
线性插值的几何意义:用
通过点
A( x0 , f ( x0 ))
和 B( x1 , f ( x1 ))
由插值条件: p( xi ) f ( xi ) (i=0,1,2,…,n),可得
a n x0 n a n 1 x0 n 1 a1 x0 a 0 f ( x0 ) n n 1 a n x1 a n 1 x1 a1 x1 a 0 f ( x1 ) a x n a x n 1 a x a f ( x ) n 1 n 1 n 0 n n n
拉格朗日插值多项式
两个插值点可求出一次插值多项式,而三 个插值点可求出二次插值多项式。插值点增加到n+1 个时,也就是通过n+1个不同的已知点 ( xi , yi )(i 0,1,, n) ,来构造一个次数为n的代数多项式P(x)。与推导线性插 值的基函数类似,先构造一个特殊n次多项式 li ( x) 的插 值问题,使其在各节点 x i 上满足
n
例6.2 已知y=f(x)的函数表 X 1 y 1
3 2
求线性插值多项式, 并计算x=1.5 的值 解: 由线性插值多项式公式得 x x0 x x1 p( x) y0 y1 x 0 x1 x1 x 0
东南大学_数值分析_第六章_常微分方程数值解法
第六章 常微分方程数值解法——RK 4法、AB 4法******（学号） *****（姓名）上机题目要求见教材P307,23题。
一、算法原理题目要求采用RK 4法和AB 4法求解最简单的常微分方程初值问题(,),()y f x y a x by a η'=≤≤⎧⎨=⎩ (1)为求解式(1)，采用离散化方法，就是寻求解)(x y 在区间],[b a 上的一系列点<<<<<n x x x x 321上的近似值 ,,,,21n y y y 。
记1(1,2,)i i i h x x i -=-=表示相邻两个节点的间距，称为步长。
求微分方程数值解的主要问题：(1) 如何将微分方程(,)y f x y '=离散化,并建立求其数值解的递推公式; (2) 递推公式的局部截断误差、数值数n y 与精确解)(n x y 的误差估计； (3) 递推公式的稳定性与收敛性. a) Runge-Kutta 方法基本思想：通过在1[,]i i x x +多预报几个点求斜率，并将其加权平均作为k *的近似值，以此构造更高精度的计算公式。
如果每步计算四次函数 的值，完全类似的，可以导出局部截断误差为)(5h O 的四阶Runge-Kutta 公式(RK 4)：1123412132431(22),6(,),(,),221(,),22(,).n n n n n n n n n n y y k k k k k f x y h h k f x y k h k f x h y k k f x h y hk +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ (2)b) Adams 显式公式Runge-Kutta 方法是单步法，计算1+n y 时，只用到n y ， 而已知信息1-n y 、2-n y 等没有被直接利用。
可以设想如果充分利用已知信息1-n y ，2-n y ，…来计算1+n y ，那么不但有可能提高精度，而且大大减少了计算量，这就是构造所谓线性多步法的基本思想。
数值分析第六章课后习题答案
第六章课后习题解答(1)()()123(1)()213(1)()()312(01.21125551154213351010(1,1,1),17( 4.0000186,2.99999k k k k k k k k k Tx x x x x x x x x x x+++ìïï=---ïïïïïï=-+íïïïïï=-++ïïïî==-(17)解：(a )因系数矩阵按行严格对角占优，故雅可比法与高斯-塞德尔均收敛。
（b )雅可比法的迭代格式为取迭代到次达到精度要求(1)()()123(1)(1)()213(1)(1)(1)312(0)(8)15,2.0000012)21125551154213351010(1,1,1),8( 4.0000186,2.9999915,2.0000012)Tk k k k k k k k k TTx x x x x x x x x x++++++-ìïï=---ïïïïïï=-+íïïïïï=-++ïïïî==-高斯塞德尔法的迭代格式为x 取迭代到次达到精度要求1212:00.40.4.0.400.80.40.80||(0.8)(0.80.32)()1.09282031,00.40.4()00.160.6400.0320.672DL U I BD L U l l l l--骣--÷ç÷ç÷ç÷ç÷=+=--ç÷ç÷÷ç÷ç÷--÷ç桫-=-+-=>-æ--çççç=-=-ççççèlJJJS解(a )雅可比法的迭代矩阵B（）BB故雅可比迭代法不收敛高斯塞德尔法迭代矩阵131()||||0.81022101220||022023002SJBDL U I BD L Ul l¥--ö÷÷÷÷÷÷÷÷÷÷ç÷ø?<骣-÷ç÷ç÷ç÷ç÷=+=--ç÷ç÷÷ç÷ç÷--ç÷桫-=骣-÷ç÷ç÷ç÷ç÷=-=-ç÷ç÷÷ç÷ç÷ç桫llSJJ SB故高斯-塞德尔迭代法收敛。
数值分析6
第六章 常微分方程的数值解法§6.1 引言在高等数学的有关常微分方程的内容中，我们主要讨论的是一些典型方程的解析解的基本求法.然而在生产实际和科学研究中遇到的微分方程往往比较复杂，在大多数情况下都不能得到问题的解析解表达式；有时即使是一些已有了求解的基本方法的典型方程，在实际使用时也有很大的困难.如求解线性常系数微分方程，当方程的阶数较高时，就涉及高次代数方程的求根问题.因此在实际问题中，对于求微分方程，一般只要求得到解在若干点上的近似值或者解的便于计算的近似表达式.本章主要考虑一阶初值问题的单步法和多步法，及其误差分析和单步法的稳定性分析；介绍二点边值问题的差分法，打靶法和有限元法.一阶初值问题的基本形式是：⎩⎨⎧∈=≤≤='],[)(),(000b a x y x y b x a y x f y 二点边值问题的一般形式为：⎩⎨⎧==≤≤'=''βα)(,)(),,(b y a y b x a y y x f y 在讨论中我们总假定初值问题和边值问题的解存在唯一，初值问题是好条件的.§6.2 初值问题的数值解法一、Euler 方法及其截断误差1. Euler 公式考虑微分方程初值问题：⎩⎨⎧∈=≤≤='],[)(),(000b a x y x y b x a y x f y （6.1）设b x x x a n =<<<= 10是区间],[b a 的一个分划.所谓微分方程的数值解，即寻找方程在节点n x x x ,,1,0 上的近似值n y y y ,,,10 .相邻2个节点间的距离i i i x x h -=+1称为步长.一般取i h 为常 数，即i h =h ，于是kh x x k +=0，n k ,,2,1,0 =将（6.1）的两端在区间],[1+i i x x 上积分，得到⎰+'1i i x x dx y =⎰+1),(i i x x dx y x f 即)(1+i x y =)(i x y +⎰+1),(i i x x dx y x f 对积分⎰+1),(i i x x dx y x f 应用左矩形公式⎰+1),(i i x x dx y x f ≈h x y x f i i ))(,(，)(1i i x x h -=+ 得到离散的近似等式)(1+i x y =)(i x y +h x y x f i i ))(,(或者1+i y =i y +h y x f i i ),( （6.2）若利用右矩形近似公式⎰+1),(i i x x dx y x f ≈h x y x f i i ))(,(11++得到相应的离散化近似等式)(1+i x y =)(i x y +h x y x f i i ))(,(11++或者1+i y =i y +h y x f i i ),(11++ （6.3）从0x 的初值0y 开始，按（6.2）式逐点计算以后各点的上的值.因此（6.2）式称为显式Euler 公式.（6.3）式中因右端式子中含有待求函数值1+i y ，不能逐步显式计算故称为隐式Euler 公式.如果将两式坐算术平均，就得到梯形公式：1+i y = i y +)],(),([211+++i i i i y x f y x f h （6.4） 公式（6.2），（6.3），（6.4）统称为差分公式，由于它们都是用i y 去计算1+i y ，故称它们为单步法. 显式Euler 公式的具有明显的几何意义：在区间],[10x x 上用过),(00y x P 以),(00y x f 为斜率的直线 y =0y +))(,(000x x y x f -， P 近似代替)(x y ，用该直线与直线1x x =的交点),(111y x P 的纵坐标01231y =0y +h y x f ),(00近似代替)(1x y .然后在区间],[21x x 上用过),(111y x P ，以),(11y x f 为斜率的直线y =1y +))(,(111x x y x f -近似代替)(x y ，用该直线与直线2x x =的交点),(222y x P 的纵坐标2y =1y +h y x f ),(11作为)(2x y 的近似值.一般地设折线已推进到),(i i i y x P ，则在区间],[1+i i x x 上用过点),(i i i y x P ，斜率为),(i i y x f 的直线y =i y +))(,(i i i x x y x f -近似替代)(x y ….依次类推，得到一条折线.因此显式Euler 公式有时又称折线法.2.隐式公式的计算对于隐式方法，如果),(y x f 是y 的线性函数则可显化计算.如5+='xy y ，其隐式Euler 公式为 1+i y =i y +h y x i i )5(11+++，简单化简得：1+i y =)1()5(1+-+i i hx h y ，依此公式可计算各点处的函数值.但当),(y x f 是y 的非线性函数时，则相应的隐式Euler 公式不能显化，而需要通过其它的方法（比如迭代法）求解.对于隐式公式的求解，每前进一步都要进行若干次迭代，而且需要提供好的迭代初始值，使迭代尽快收敛.通常有两种方法：一是先由显式的Euler 公式提供一个迭代初值，再用隐式公式进行简单迭代.以梯形公式为例，其简单迭代为： ⎪⎩⎪⎨⎧++=+=+++++)],(),([2),()(11)1(1)0(1k i i i i i k i i i i i y x f y x f h y y y x hf y y ， ,2,1,0=k （6.5） 反复迭代直到满足精度ε<-+++)(1)1(1k i k i y y但是要使这种迭代收敛，其步长h 应满足一定条件.这可以由),(y x f 关于y 满足Lipschtiz 条件得到.称函数),(y x f 关于y 满足Lipschtiz 条件，是指存在常数0>L ，使得 2121),(),(y y L y x f y x f -≤-对任意的1y ，2y 成立.于是由（6.4）与（6.5）得到)1(11+++-k i i y y =)],(),([2)(1111k i i i i y x f y x f h +++++ 从而)1(11+++-k i i y y =)],(),([2)(1111k i i i i y x f y x f h ++++- ≤)(112k i i y y L h ++-≤…≤)0(11)2(++-i i k y y L h 这说明，只有当12<hL ，即L h 2<迭代序列才收敛. 二是可以采用所谓的预测-校正技术.其思想是在每步迭代前用显式方法预测一个值p i y 1+，然后用隐式公式校正一次（迭代一次）.如显式Euler 公式预测，梯形公式校正，其格式如下：⎪⎩⎪⎨⎧++=+=++++)],(),([2),(1111p i i i i i c i i i i p i y x f y x f h y y y x hf y y （6.6） 再用c i y 1+近似替代1+i y .这等价于1+i y = i y +))],(,(),([21i i i i i i y x hf y x f y x f h +++ （6.7） 公式（6.7）显然是一个显式的差分方程，称为改进的Euler 公式.注意这里按公式（6.7）计算得到的1+i y 实际上是真实的1+i y 的一个近似c i y 1+，在实际应用中为了便于编程计算，可将公式（6.7）变形为 ⎪⎪⎪⎩⎪⎪⎪⎨⎧+=+=+=++)(21),(),(11c p i p i i c i i i p y y y y x hf y y y x hf y y （6.8）1Exp 求解初值问题（取步长为0.1）⎩⎨⎧=≤≤-='1)0(102y x xy y解：不难求得该初值问题的准确解为2x e y -=，下面用显式Euler 公式及改进的Euler 公式分别计算并比较精度.用显式方法计算：显式Euler 公式的具体形式为 1+i y =i y +)2(i i y x h -=i i y x )2.01(-现将计算结果列于下表：如果采用改进的公式，其具体形式为：⎪⎪⎪⎩⎪⎪⎪⎨⎧+=+-=-+=-=-+=++)(21)1.0(2.0)2()2.01()2(11c p i p i i p i i c i i i i i p y y y y x y y x h y y x y y x h y y将计算结果列于下表：比较两表的计算结果，可以看出改进的公式明显提高了计算精度.可见不同的方法具有不同的精度，需要讨论这种方法的误差.3. 局部截断误差和方法的阶一般地，单步公式的形式可表示为⎩⎨⎧=+=+++00111)(),,,,(y x y h y y x x h y y i i i i i i ϕ，1,,2,1,0-=n i （6.9） 称),,,(1h y y x i i i +ϕ为增量函数.从0x 开始计算，每一步都会产生误差，如果迭代到n x ，则此时有误差n ε=n n y x y -)(.n ε称为方法在n x 点的整体截断误差.显然要求出整体截断误差n ε是复杂的，因为这里牵涉到多次多步迭代，所以我们转而考虑从n x 到1+n x 的局部情况.假定在n x 处n y 没有误差，即n y =)(n x y ，用1+n T 表示从n x 到1+n x 的局部误差，则由（6.9）式，1+n T 可表示成1+n T =)(1+n x y -1+n y=)(1+n x y -)(n x y -)),(),(,,(11h x y x y x x n n n n ++ϕ于是得到局部截断误差的概念.定义1：设)(x y 是微分方程的精确解，则1+n T =)(1+n x y -)(n x y -)),(),(,,(11h x y x y x x n n n n ++ϕ称为单步法的局部截断误差.局部截断误差与整体截断误差是紧密联系的，在一定条件下，如果局部截断误差是)(1+p h O ，则整体截断误差是)(p h O .为此我们给出定义2.定义2：如果给定方法的局部截断误差满足1+n T =)(1+p h O其中1≥p 为整数，则称该方法是p 阶的，或具有p 阶精度.一般来说，一个差分公式的阶p 越大，方法的精度越高.定义3：若一个p 阶单步法的局部截断误差为；1+n T =1))(,(+p n n h x y x g +)(2+p h O则称其第一非零项1))(,(+p n n h x y x g 为该方法的局部截断误差主项.下面我们推导Euler 公式的局部截断误差.对显式Euler 公式，其局部截断误差为：1+n T =)(1+n x y -)(n x y -h x y x f n n ))(,(=)(1+n x y -)(n x y -)(n x y h '=)(n x y +)(n x y h '+)(22n x y h ''+)(3h O -)(n x y -)(n x y h ' =)(22n x y h ''+)(3h O 这里应用了对精确解)(n x y 满足原微分方程)(n x y '=))(,(n n x y x f .因此，显式Euler 公式是具有一阶精度的方法，其局部截断误差的主项为)(22n x y h ''.采用相似的方法可以推出隐式Euler 公式也是一阶方法其截断误差主项为-)(22n x y h ''，与显式相比只差一个符号.但梯形公式是二阶方法，其截断误差主项为-)(123n x y h '''.二、龙格-库塔)(Kutta Runge -方法1.Kutta Runge -方法的基本思想显然，方法的阶数越高，其精度也就越高.显式、隐式的Euler 公式都是一阶方法，梯形公式是二阶方法，能否进一步提高方法的阶呢？我们从研究差商hx y x y i i )()(1-+开始.由微分中值定理hx y x y i i )()(1-+=)(h x y i θ+'，)10(<<θ 利用微分方程),(y x f y ='，得到hx y x y i i )()(1-+=))(),(h x y h x f i i θθ++ 或等价于)(1+i x y -)(i x y =))(),(h x y h x f h i i θθ++⇒)(1+i x y =)(i x y +))(),(h x y h x f h i i θθ++ （6.10）（6.10）式中的))(),(h x y h x f i i θθ++称为区间),(1+i i x x 上的平均斜率，记为*k ，即 *k =))(),(h x y h x f i i θθ++因此只要对平均斜率*k 提供一种算法，则由（6.10）便可以得到微分方程的一个数值计算公式.用这个观点来考虑Euler 公式及改进的Euler 公式，我们发现由于Euler 公式中仅取i x 一个点的函数值),(i i y x f 作为平均斜率*k 的近似值，因而精度较低.而改进的Euler 公式却用了两个点i x 与1+i x 处的函数值),(i i y x f 与),(1+i i y x f 的平均数作为平均斜率*k 的近似值，因而精度达到了二阶.这启发了我们，如果在区间],[1+i i x x 上多取几个函数值，取它们的加权平均作为平均斜率*k 的近似值，则有可能构造出具有更高精度的数值微分计算公式.这就是Kutta Runge -方法的基本思想.2.二阶Kutta Runge -公式首先推广改进的Euler 公式，我们考虑在区间],[1+i i x x 内任取一点l i x +=i x +h l ，)10(≤<l我们希望取i x 与l i x +两点处的斜率值1k ，2k 的加权平均作为平均斜率*k 的近似值，即 *k ≈11k λ+22k λ从而离散后的微分方程为1+i y =i y +h （11k λ+22k λ）这里1λ，2λ是待定参数.同Euler 公式一样，我们任取1k =),(i i y x f ，问题是如何预测l i x +处的斜率值2k ？仿照改进的Euler 公式做法，我们可以先用l i x +提供)(l i x y +的预测值l i y +=i y +h l ),(i i y x f然后，通过预测值l i y +计算f 来产生斜率值2k =),(l i l i y x f ++.其具体计算公式如下：⎪⎩⎪⎨⎧+===++=++++),(),(),()(12122111lhk y x f y x f k y x f k k k h y y i l i l i l i i i i i λλ （6.11） 公式（6.11）中含有3个待定参数1λ，2λ，l .我们希望适当选取这些参数值，产生具有二阶精度的数值微分计算公式.为此需要考虑其局部截断误差1+i T =)(1+i x y -)(i x y -))],(,())(,([21i i i l i i i y x lhf y x f x y x f h +++λλ利用)),(,(i i i l i y x lhf y x f ++的二元泰勒展开公式及))(,()(i i i x y x f x y ='可以得到：1+i T =)()1(21i x y h '--λλ+)(2122i x y l h ''⎪⎭⎫⎝⎛-λ+)(3h O 要使公式具有二阶精度，即1+i T =)(3h O ，只需 ⎪⎩⎪⎨⎧==+211221λλλl （6.12） 注意这里共有3个待定参数，但只有两个方程，因此有一个自由度，即满足公式（6.12）的参数不是一组，而是很多，我们统称为二阶Kutta Runge -公式.公式（6.12）的解为 l 2111-=λ，l212=λ 特别当l =1时，即l i x +=1+i x ，1221λλ==，二阶Kutta Runge -公式成为改进的Euler .如果取21=l ，则12=λ，1λ=0，这时的二阶Kutta Runge -公式称为变形的Euler 公式（也称中点公式）.其形式是：1+i y =i y +))],(2,2([i i i i y x f h y h x f h ++ 如果为了使用方便，也可写成方程组形式： ⎪⎪⎪⎩⎪⎪⎪⎨⎧++===+=+++)2,21(),(),(121212121k h y x f y x f k y x f k hk y y i i i i i i i i （6.13）从表面上看，中点公式仅含一个斜率值2k ，但2k 是通过1k 计算得到的，因此每完成一步仍然需要两次计算函数值，工作量和改进的Euler 公式相同.。
数值分析第六章1
若Ln ( x )是一个次数不超过n的多项式 Ln ( x ) a0 a1 x an x n (其中ai 为实数) 则称Ln ( x )为f ( x )的n次Lagrange插值多项式.
定理1 (插值多项式存在唯一性) 满足插值条件 Ln ( xk ) yk , k 0,1,, n 的次数 不超过n的多项式
( n1) | f ( x) | 在区间[a, b]的一个上界. 其中M是
拉格朗日插值公式的缺点： 当对于给定的若干个节点求出其插值函 数后，如果想再增加一个节点并求出新的插 值函数时，所有的插值基函数都必须重新计 算. 因此只适用于插值节点已给定的情形, 若 插值过程中需要增加新节点, 就需要考虑其 它插值公式.
M3 | R2 ( x ) | |(x x0 )(x x1 )(x x2 )|, 6 M 3 cos( x0 ) 0.828, 1 | R2 (0.3367) | 0.828 0.0167 0.0033 0.0233 6 0.178 10 6.
定理3 | x i x i 1 | 在定理2条件成立的情况下, 令 h max 1 i n M n 1 | Rn ( x ) | h 则 4( n 1)
M2 | R1 ( x ) | | ( x x0 )( x x1 ) |, 2
M 2 max | f ( x ) | sin x1 0.3335,
x0 x x 2
1 | R1 (0.3367) | 0.3335 0.0167 0.0033 0.92 10 5. 2
数 值 分 析
第六章 插值方法
一、 Lagrange插值
二、 Newton插值
三、 Hermite插值 四、 分段低次插值 五、 三次样条插值
数值分析课件第6章
a22 x2
a2n xn
b2
an1 x1 an2 x2 ann xn bn
就是 Ax b 进行矩阵分裂 A M - N
且M为可选择的非奇异矩阵，且使得Mx d容易求解于是
Ax b x M 1Nx M 1b
可得一阶定常迭代法
取 初 始 向 量x (0)， x(k1) Bx(k) f , k 0,1,, 其中B M 1N M 1 (M A) I M 1A, f M 1b.
取 1.3, 第11次迭代 结果为
x(11) (0.99999646, 1.00000310, 0.99999953,
0.99999912)T
此时， (11) 0.46 105 2
注：迭代法求解线性方程组在计算机实现时可用
x(k 1) x(k )
max
1in
x(k 1) i
x(k) i
机动 上页 下页 首页 结束
本章将介绍迭代法的一般理论及雅可比迭代法、高斯—塞
德尔迭代法、超松弛迭代法，研究它们的收敛性。
例6-1 求解线性方程组
8x1 3x2 2x3 4x1 11x2 x3
20, 33,
6x1 3x2 12x3 36.
解：方程记为 Ax b，即
8 A 4
6
3 2
(1
x(k 1) 1
4 x2(k )
x(k ) 3
x4(k ) ) /
4
x(k 1) 3
x(k ) 3
(1
x(k 1) 1
x(k 1) 2
4x3(k )
x4(k ) ) / 4
x4(k 1)
x(k ) 4
(1
x(k 1) 1
x(k 1) 2
数值分析第六章小结【计算方法】
第六章 数值积分学习小结一、本章知识梳理 求积公式及其代数精度：数值求积公式的一般形式：()0()()nbn k k ak f x dx f x λ=≈∑⎰截断误差：)()(0k bank k n x f dx x f R ⎰∑=-=λ数值求积公式是一种近似方法，因此，要求它对尽可能多的被积函数f 能准确计算积分的值，这就有了代数精度的概念。
定义：对于上面所列的求积公式，当()f x 为任何次数不高于m 的多项式时都成为等式，而当()f x 为某个m+1次多项式时不能成为等式，则称它具有m 次代数精度。
插值型求积公式：()0()()nbn k k ak f x dx f x λ=≈∑⎰其中()()(0,1,...,)bn kk al x dx k n λ==⎰截断误差：(1)0()[()](1)!n nbn j aj f R x x dx n ξ+==-+∏⎰定理：n+1个节点的插值型求积公式至少具有n 次代数精度。
Newton —Cotes 求积公式：求积公式：()0()()nbn kak b a f x dx f a k nλ=-≈+∑⎰其中()()()(0,1,...,)n n k k b a c k n λ=-=()00(1)[()]!()!n kn n n kj j kct j dt k n k n -=≠-=--∏⎰ 2(1)00()[()](1)!n n n n n j h R f t j dt n ξ++==-+∏⎰ 梯形公式(n=1)：()[()()]2bab af x dx f a f b -≈+⎰31()''(),(,)12b a R f a b ηη-=-∈Simpson 公式(n=2)：()[()4()+()]62bab a a bf x dx f a f f b -+≈+⎰5(4)2()(),(,)2880b a R f a b ηη-=-∈Simpson3/8公式(n=3)：22()[()3()3()()]833bab a a b a bf x dx f a f f f b -++≈+++⎰5(4)3()(),(,)6480b a R f a b ηη-=-∈Cotes 公式(n=4)：33()[7()32()12()32()7()]90424bab a a b a b a bf x dx f a f f f f b -+++≈++++⎰7(6)4()(),(,)1935360b a R f a b ηη-=-∈复化求积法：复化梯形公式：(0,1,...,),k b ax a kh k n h n-=+==11()[()()2()]2n bak hf x dx f a f b f a kh -=≈+++∑⎰2()''(),[,]12T b a R h f a b ηη-=-∈复化Simpson 公式：(0,1,...,2),2k b ax a kh k m h m-=+==121211()[()()4()2()]3m m bi i ai i hf x dx f a f b f x f x --==≈+++∑∑⎰4(4)()(),[,]180S b a R h f a b ηη-=-∈ Gauss 型求积公式：一般理论：1()()()nbi i ai x f x dx A f x ρ=≈∑⎰()()(1,2,...,)()'()bn i ai n i x x A dx i n x x x ρωω==-⎰()()()()!n bn af R x x dx n ξρω=⎰几种Gauss 型求积公式：Gauss-Legendre 求积公式：111()()ni i i f x dx A f x -=≈∑⎰222(1,2,...,)(1)['()]i i n i A i n x L x ==-(2)242()2(!)2(1,1)(2)![(2)!]21n n f n R n n n ηη=••∈-+， Gauss-Laguerre 求积公式：1()()nxi i i e f x dx A f x ∞-=≈∑⎰22(!)(1,2,...,)['()]i i n i n A i n x U x == 二、本章测验题已知Legendre 正交多项式)(x L n 有三项递推关系式：⎪⎪⎩⎪⎪⎨⎧=+-++===-+ ,2,1)(1)(112)()(,1)(1110n x L n n x xL n n x L x x L x L n nn试推导两点Gauss-Legendre 求积公式)()()(221111x f A x f A dx x f +≈⎰-的求积系数和节点，并用此公式计算下列积分的近似值。
数值分析第六章线性方程组迭代解法
数值分析第六章线性方程组迭代解法线性方程组是数值分析中的重要内容之一，其求解方法有很多种。
其中一种常用的方法是迭代解法，即通过不断迭代逼近方程组的解。
本文将介绍线性方程组迭代解法的基本思想和常用方法。
线性方程组可以用矩阵形式表示为Ax=b，其中A是系数矩阵，b是常数向量，x是未知向量。
线性方程组的解可以是唯一解，也可以是无穷多个解。
迭代解法的基本思想是通过不断迭代，并利用迭代序列的极限，逼近线性方程组的解。
迭代解法适用于大型的线性方程组，而直接求解法则适用于小型的线性方程组。
常用的迭代解法有雅可比迭代法、高斯-赛德尔迭代法和逐次超松弛迭代法。
雅可比迭代法是最简单的线性方程组迭代解法之一、它的基本思想是将线性方程组的每个方程都单独表示为未知数x的显式函数，然后通过不断迭代求解。
雅可比迭代法的迭代公式为：x(k+1)=D^(-1)(b-(L+U)x(k))其中，D是A的对角元素构成的对角矩阵，L是A的下三角矩阵，U 是A的上三角矩阵，x(k)是第k次迭代的解。
高斯-赛德尔迭代法是雅可比迭代法的改进版。
它的基本思想是将每个方程的解带入到下一个方程中，而不是等到所有方程都迭代完毕后再计算下一组解。
高斯-赛德尔迭代法的迭代公式为：x(k+1)=(D-L)^(-1)(b-Ux(k))其中，D是A的对角矩阵，L是A的下三角矩阵（除去对角线），U是A的上三角矩阵（除去对角线），x(k)是第k次迭代的解。
逐次超松弛迭代法是对高斯-赛德尔迭代法的改进。
它引入了松弛因子w，通过调节松弛因子可以加快收敛速度。
逐次超松弛迭代法的迭代公式为：x(k+1)=(D-wL)^(-1)[(1-w)D+wU]x(k)+w(D-wL)^(-1)b其中，D是A的对角矩阵，L是A的下三角矩阵（除去对角线），U是A的上三角矩阵（除去对角线），w是松弛因子，x(k)是第k次迭代的解。
线性方程组迭代解法需要设置迭代停止准则，通常可以设置迭代次数上限或者设置一个精度要求。
数值分析方法第六章
向前差商
f ( x0 h) f ( x0 ) f ' ( x0 ) h
由Taylor展开
h2 f ( x0 h) f ( x0 ) hf ' ( x0 ) f ' ' ( ), x0 x0 h 2!
误差
f ( x0 h) f ( x0 ) h R( x) f ' ( x0 ) f ' ' ( ) O(h) h 2!
7
向后差商
f ( x0 ) f ( x0 h) f ' ( x0 ) h
由Taylor展开
h2 f ( x0 h) f ( x0 ) hf ' ( x0 ) f ' ' ( ), x0 x0 h 2!
误差
f ( x0 ) f ( x0 h) h R( x) f ' ( x0 ) f ' ' ( ) O(h) h 2!
xi 1 0.02k (k 0,1,,100)
处的近似值。
20
21
MATLAB中的数值微分
数值微分的实现
在MATLAB中，没有直接提供求数值导数的函数，只有
计算向前差分的函数diff，其调用格式为：
DX=diff(X)：计算向量X的向前差分， DX(i)=X(i+1)-X(i)，i=1,2,…,n-1。 DX=diff(X,n)：计算X的n阶向前差分。 diff(X,2)=diff(diff(X))。
其中 在x0 , x1 ,, xn 之间，上式第一项是 f ( xi ) 的近 似值。
12
300
Seebeck Coefficient a, mV/K
数值分析-第六章-数值积分
k 0
而对应的误差为
b
b f (n1) ( )
I In
(
a
f
(
x)
Ln
(
x))dx
a (n 1)! wn1(x)dx
Newton-Cotes公式
当节点为等距节点时，对应的插值型求积公式称为 Newton-Cotes 公式。
梯形公式：最简单的 Newton-Cotes 公式
a
2
梯形公式的误差
梯形公式的误差为：
b f ( )
E I T a 2 (x a)(x b)dx
注意到对任意的 x [a,b] ，有 (x a)(x b) 0，根据积分中值定理，
若 f "(x) C[a,b] ，有
E f ()
b
(x a)(x b)dx
第六章 数值积分
数值积分的基本概念 数值积分的基本思想 代数精度 插值型求积公式
Newton-Cotes 求积公式 梯形公式、辛普森公式、一般的 Newton-Cotes 公式 复化积分公式：复化梯形公式、复化辛普森公式 区间逐次分半法
Romberg（龙贝格）积分
高斯型求积公式
数值积分的基本概念
微积分中定积分的定义为：b Nhomakorabean
a
f
(x
)dx
lim
n m a xxk
k01
xk
f
k( ，)
n
b
n
可用 xk f (xk ) 作为原积分的近似： a f (x)dx xk f (xk ) 。
k 1
k 1
进一步推广得到更一般的公式：
数值分析
这里h为常数，称为步长， , 和 分别为向前, 向后和中心差分算子.
以向前差分为例 利用一阶差分可定义二阶差分
2 f k f k 1 f k f k 2 2 f k 1 f k .
还可以定义m阶差分 m f k m 1 f k 1 m 1 f k ;
xi ƒ(xi) x0 x1 x2 x3 xn ƒ(x0) ƒ(x1) ƒ(x2) ƒ(x3) ƒ(xn) 一阶 均差 二阶均差 三阶均差 n阶均差
ƒ[x0, x1] ƒ[x1, x2] ƒ[x0, x1, x2] ƒ[x2, x3] ƒ[x1, x2, x3] ƒ[x0, x1, x2, x3] ƒ[xn-1, xn] ƒ[xn-2, xn-1, xn] ƒ[xn-3, xn-2, xn-1, xn] ƒ[x0, x1,…, xn]
0 ( x) 1 1 ( x ) x x0 2 ( x ) ( x x0 )( x x1 )
LL n ( x ) ( x x0 )( x x1 ) L ( x xn1 )
当增加一个节点xn+1时，只需加上基函数
n1 ( x xi ) 即可.
n n n j
f k ( I E ) f k (1)
1 2
n1
………… f [ x, x0 , ... , xn1 ] f [ x0 , ... , xn ] ( x xn ) f [ x, x0 , ... , xn ]
1 + (x x0) 2 + … … + (x x0)…(x xn1)
n1
f [ x, x0 ] f [ x0 , x1 ] ( x x1 ) f [ x, x0 , x1 ]
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
中心差商
f ( x0 h) f ( x0 h) f ' ( x0 ) 2h
由Taylor展开
h2 h3 f ( x0 h) f ( x0 ) hf '( x0 ) f ''( x0 ) f '''(1 ), x0 1 x0 h 2! 3! h2 h3 f ( x0 h) f ( x0 ) hf '( x0 ) f ''( x0 ) f '''( 2 ), x0 h 2 x0 2! 3!
[ x0 , x2 ]
16
类似地，可以得到n=2的带余项的二阶三点公式
1 h 2 (4) f ( x0 ) h 2 [ f ( x0 ) 2 f ( x1 ) f (x2 )] [hf (1 ) 6 f ( 2 )] 1 h 2 (4) f ( ) f ( x1 ) 2 [ f ( x0 ) 2 f ( x1 ) f (x2 )] h 12 1 h 2 (4) f ( 2 )] f ( x2 ) 2 [ f ( x0 ) 2 f ( x1 ) f ( x2 )] [hf (1 ) h 6
26
数值积分雏形 数格子
例：计算 R1 R R2 的 颗粒数量(分数)
分布密度，f
F
R2 R1
f ( R)dR
F
R1
R2
颗粒尺寸，R
原始方法：曲线下面面积所占方格数目
27
定积分的计算
当 f (x) 是闭区间上的连续函数时，如何 计算下列定积分？
I ( f ) f ( x)dx
23
例 用不同的方法求函数f(x)的数值导数，并在同一个坐标系中 做出f'(x)的图像。 程序如下： f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2'); g=inline('(3*x.^2+4*x-1)./sqrt(x.^3+2*x.^2x+12)/2+1/6./(x+5).^(5/6)+5'); x=-3:0.01:3; p=polyfit(x,f(x),5); %用5次多项式p拟合f(x) dp=polyder(p); %对拟合多项式p求导数dp dpx=polyval(dp,x); %求dp在假设点的函数值 dx=diff(f([x,3.01]))/0.01; %直接对f(x)求数值导数 gx=g(x); %求函数f的导函数g在假设点的导数 plot(x,dpx,x,dx,'.',x,gx,'-'); %作图
7
向后差商
f ( x0 ) f ( x0 h) f ' ( x0 ) h
由Taylor展开
h2 f ( x0 h) f ( x0 ) hf ' ( x0 ) f ' ' ( ), x0 x0 h 2!
误差
f ( x0 ) f ( x0 h) h R( x) f ' ( x0 ) f ' ' ( ) O(h) h 2!
试计算美国20世纪的年增长率
若记t时刻的人口为 x(t ),则人口的增长率为
dx dt r (t ) x(t )
dx 如何求 dt
2
在实际问题中，往往会遇到某函数f(x)是用表格 表示的,用通常的导数定义无法求导,因此要寻求其他 方法近似求导。常用的数值微分方法有: 一. 运用差商求数值微分 二. 运用插值函数求数值微分 三. 运用样条插值函数求数值微分
10
插值型求导公式
主要方法：插值多项式，讨论函数f(x)导数的近似值求法。
定义：若函数f(x)在节点 xi (i 1,2,, n 1) 处的函数 值已知，作f(x)的n次插值多项式 Pn ( x) ，并用 Pn ( x) 近似代 替f(x)，即
f ( x) Pn ( x)
由于Pn ( x) 是多项式，容易求其导数，故对应于f(x) 的每一个插值多项式 Pn ( x) ，建立一个数值微分公式
f ( x) Pn( x)
这样建立起来的数值微分公式，称为插值型数值微分公式。 11
带有余项的数值微分公式
即使 Pn ( x) 与f(x)的函数值处处相差甚微，但两个 函数的导数在某些点上的值仍可能有很大的差异 所以要对误差进行分析
n1 ( x) df (n1) ( ) f ( n1) ( ) 1 ( x) f ( x) Pn( x) n (n 1)! (n 1)! dx
10
0
1
2
annealing time t, 1000 s
5
10
20
50
100
13
如果只计算结点处的导数值
f ( xi ) Pn( xi )
数值微分公式的余项为
f ( ) 1 ( xi ) f ( xi ) Pn( xi ) n (n 1)!
( n 1)
14
带余项的两点数值微分公式
22
例 设x由[0,2π]间均匀分布的10个点组成，求sinx的 1~3阶差分。
命令如下： X=linspace(0,2*pi,10); Y=sin(X); DY=diff(Y); %计算Y的一阶差分 D2Y=diff(Y,2); %计算Y的二阶差分，也可用命令diff(DY) 计算 D3Y=diff(Y,3); %计算Y的三阶差分，也可用diff(D2Y)或 diff(DY,2)
xi 1 0.02k (k 0,1,,100)
处的近似值。
20
21
MATLAB中的数值微分
数值微分的实现
在MATLAB中，没有直接提供求数值导数的函数，只有
计算向前差分的函数diff，其调用格式为：
DX=diff(X)：计算向量X的向前差分， DX(i)=X(i+1)-X(i)，i=1,2,…,n-1。 DX=diff(X,n)：计算X的n阶向前差分。 diff(X,2)=diff(diff(X))。
Chapter 5 数值微分与数值积分
§1 数值微分
先看一个实例: 已知20世纪美国人口的统计数据为(单位:百万)
年份 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 人口 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4
a b
积分中值定理
I ( f ) f ( x)dx (b a) f ( )
a
28
b
上述公式的几何意义
f ( ) 的矩形面积恰好等于 所求曲边梯形的面积 I ( f ) 。 y
当h充分小时 , 可用差商来逼近导数
5
差商公式
f ( a h) f ( a ) 向前差商近似计算： f (a ) h
f ( a ) f ( a h) 向后差商近似计算： f (a ) h
f ( a h) f ( a h) 中心差商近似计算： f (a) 2h
误差
f ( x0 h) f ( x0 h) R( x) f ' ( x0 ) 2h h2 h2 [ f ' ' ' (1 ) f ' ' ' ( 2 )] f ' ' ' ( ) O(h 2 ) 12 6
9
分析
1、步长过大则截断误差显著，但如果步长太小又会导 致舍入误差的增长。 2、在实际计算时，在保证截断误差满足精度要求的前 提下选取尽可能大的步长 3、事先给出一个合适的步长往往是很困难的。通常在 在变步长的过程中实现步长的自动选择。
其中 在x0 , x1 ,, xn 之间，上式第一项是 f ( xi ) 的近 似值。
12
300
Seebeck Coefficient a, ) A B exp( 30
250 200
误差放大 !!!
150 20
da dt
100 50 0
da dt
注意
S "( xi ) M i
19
1 例： 利用函数 f ( x) 2 在节点 1 25 x
S’(-1)=0.0740, S’(1)=-0.0740
xi 1 0.1i (i 0,1, ,20) 上的函数值和边界条件
构造三次样条插值函数S(x)，并用它来计算f(x)和 f’(x)在下列点
（n=1）
x x0 x x1 p1 ( x) f ( x0 ) f ( x1 ) x0 x1 x1 x0
f ( x1 ) f ( x0 ) f ( x0 ) h f ( x1 ) f ( x0 ) f ( x1 ) h f ( x1 ) f ( x0 ) h f ( x0 ) f ( ) h 2 , [ x0 , x1 ] f ( x1 ) f ( x0 ) h f ( x1 ) f ( ) h 2
15
类似地，可以得到n=2的带余项的三点公式
1 h2 f ( x0 ) 2h [3 f ( x0 ) 4 f ( x1 ) f (x2 )] 3 f ( ) 1 h2 f ( x1 ) [ f ( x0 ) f ( x2 )] f ( ) 2h 6 1 h2 f ( x2 ) [ f ( x0 ) 4 f ( x1 ) 3 f ( x2 )] f ( ) 2h 3
四. 运用数值积分求数值微分