基于 MATLAB的数值分析以软件 MATLAB作为辅助工具介绍数值分析（科学与工程计算）的基本内容,注重讲授一些求解方程以及结果可视化的知识和技巧,使同学们能够有效地解决问题并处理计算结果。
内容包括：
1,MATLAB编程和绘图
2、数值分析的数学基础
3、数值算法在工程、科学中的应用第一章 MATLAB入门
1,MATLAB的命令窗口工作都在此处完成。
2、怎样进行计算运算对象：矩阵算术运算符,+ - * ^,*,/,^
倒除,\ a\b=b/a
变量与变量名,变量名和变量名类型不需声明 。
3、数据显示格式默认格式,5位（ format short)
format long 16位
format e 短的浮点格式
format long e 长的浮点格式
4、清除命令
clear:清除所有使用过的变量或某个（些）变量
clc,清除命令窗口
5、程序结构分支,if __else__end; if__elseif__end; if__break__end
循环,for__end; while__end
6、读写输入数据,z=input(‘type youe input:’) 键盘输入格式化输出,fprintf(‘e_format %12.5e \n’,vol)
7、数学函数
8、功能函数
sort(x) sum(x) max(x) min(x) mod(x,y) rand(n) eval(s)
9、编程（编写 M文件）
10、绘图第二章 数值代数内容：数值代数就是研究有关 矩阵计算 的问题。主要包括：
1,线性代数方程组的求解；
2,矩阵特征值问题要求,1、掌握用 MATLAB求解的方法
2、知道那些问题是困难的，那些问题是不可解的。
A=zeros(m,n) m行 n列的零矩阵
I=eye(n) n阶单位矩阵
A=ones(m,n) 元素均为 1
A’ A的转置
A(:,k) A(k,:) A(m1:m2,n1:n2)
inv(A) A的逆
size(A) A的大小
hilb(n) Hilbert矩阵
2.1 矩阵
bAx?
情况 1,m=n(正规方程），最常见；
情况 2,m<n(不定方程）；
情况 3,m>n(超定方程）；
本节只介绍情况 1。
倍计算时间大约是上面的效率最高
50 *)(
\
bAi n vx
bAx
MATLAB命令：
2.2 解线性代数方程组的 MATLAB命令线性代数方程组并不总是数值可解的。只有当矩阵 A
的行列式不为零时才行！矩阵 A的行列式即使不为零，但当很小或很大时，解的误差可能很大。
计算矩阵行列式的 MATLAB命令：
)d e t( AD?
2.4 病态问题有许多线性代数方程组理论上是可解的，但实际计算中由于受到舍入误差的影响而无法得到精确解。此类问题成为病态问题。
病态问题的计算过程中，小的舍入误差或系数矩阵的微小变化都可能使解产生很大误差。（例子 P97)
2.3 不可解问题病态矩阵的一个重要标志是条件数：
1)( AAAc o n d
MATLAB命令,)( Ac ond
当矩阵是病态时，其条件数一定很大，但它并不能直接说明解的误差。
线性方程组解的误差程度也取决于计算环境的精度。
条件数和行列式与计算环境是相互独立的。所以大条件数或小行列式未必意味无法直接精确求得线性方程组的解，
它只意味着有很大误差可能。而实际上如果采用更高精度的计算环境则很可能得到非常满意的解。
Hilbert矩阵是非常著名的病态矩阵（ hilb(n)),它经常用来检验算法的数值稳定性的好坏。
两种原因使我们想了解求解线性代数方程组的算法。
一是实际工作中要用其它计算机语言（ Fortran&C等）编写应用程序；二是 MATLAB处理大型稀疏矩阵方程组显得很笨拙或无能为力。
nnn
n
n
ba
baa
baaa
bABbAx
2222
111211
() ( 初等行变换）
2.5 线性代数方程组的求解方法（算法）
由线性代数的理论：
下面讨论实现过程：
第一步：消元。进行到第 k步时必有
nnnknnk
knkkkkk
kknkkkk
nkk
baaa
baaa
baaa
baaaa
1,
1,11,1,1
1,
111,1111
a(k,k)作为主元，第 k行依次乘 a(i,k)/a(k,k)加到第 i行,i=k+1:n。
总共 n-1步完成,k=1:n-1.
一,Gauss消去法
e nd
e nd
e nd
e nd
jkakiajiajia
:nkjf o r
kbkiaibib
kkakiakia
nkif o r
kkaif
nkf o r
);,(),(),(),(
1
);(),()()(
);,(/),(),(
:1
e n d ;b r e a k,,0),(
1:1
当 a(k,k)=0，则上述消去法无法进行；或当其绝对值相对太小可能会出现大的计算误差。选主元法可避免这种情况。下面介绍常用的 按列选主元的 Gauss法 。
e nd
e nd
e nd
jkakiajiajia
:nkjf or
kbkiaibib
kkakiakia
nkif or
tbkk btankkka
kkbkbnkkkankka
b( k ) tba( k,k,n) ;ta
kkkkk
knkaabskkq
nkf or
);,(),(),(),(
1
);(),()()(
);,(/),(),(
:1;)(;):,(
);()();:,():,(;1
) ) ) ;,:((m a x (],[
1:1
列主元 Gauss法运算量 (只考虑乘除运算）：
第 k步 =n-k次除法
+ n-k次乘法
+（ n-k)*(n-k)次乘法 ;
)(
3
1
)(2)(
3
23
1
1
1
1
2
no
nnn
knkn
n
k
n
k
总的乘除运算量 =
第二步：回代求解
1:1);,(/)),((
);,(/
1
niiiaxjiabx
nnabx
n
ij
jii
nn
e n d
iiadibib
e n d
jbjiadd
nijf o r
d
nif o r
nnanbnb
);,(/))(()(
);(),(
:1;0
1:1:1
);,(/)()(
%%%%%%%%
二,LU分解法
LU分解的目的是将矩阵 A转换为两个矩阵的乘积，即
。方程为上三角矩阵。此时，是单位下三角矩阵，
yUxbLybLU xbAx
LLUA
,
U,
好处是：对于线性方程组，如果需要多次求解不同的非齐次项，此时 LU分解的效率将大大超过高斯消去法。
LU分解的 MATLAB命令,[l,u]=lu(A) 和 [l,u,p]=lu(A)
前面讲到的不选主元的高斯消去法和列主元高斯消去法将能实现 LU
分解。
不选主元的高斯消去法用于下面两类矩阵肯定能成，即 严格对角占优矩阵或对称正定矩阵,其他矩阵就难说了！
列主元高斯消去法是解决一般中小型稠密矩阵方程组最有效的方法之一。下面讲解列主元高斯消去法实现 LU分解的算法。
.1:1
) ) ),,:((m a x (],[,1)(
.
1
),(
),(
1
),(
),1(
1
1
i)jP( i,
nk
knkaa b skkqkkkkkkkk
kka
kna
kka
kka
L
j
k
为其中行交换后的矩阵，行与第为单位矩阵第记
1,LU分解的代数理论
) ),(,(
,
333
22322
1131211
11222211
3
333
22322
1131211
1122
2
332
222
11211
11
kkkkPP
U
a
aa
aaa
aaaa
APLPLPLPL
aa
aa
aaa
aaaa
APLPL
aa
aa
aa
aaa
APL
k
nn
n
n
n
nnnn
nnn
n
n
n
nnn
n
n
n
其中，
于是：
,
.
~
,
,
~
,,
~
,
~;)(
)(
)(
))()((
121
121211
121121
121211
121
12121
121121
123211211
11222211
PPPP
PPLPPL
PPPLPPPL
PLPLLL
UAPPP
PPLPP
PPPLPPP
PPLPPPLPL
UAPLPLPLPL
n
nn
nkkkkknk
nnnnnn
nn
nn
nkkkkkn
nnnnnnnnn
nnnn
令：
.
~~~
,
.
~~~
1
1
1
2
1
1
121
为单位下三角矩阵。可以证明其中
L
LLLLLUPA
UPALLL
nn
nn
现在我们只要将列主元高斯消去法稍加改造即是 LU分解的算法。
列主元高斯消去法的矩阵表示：
.),,,
),:1(
U
))(,()(
.
'
,2,1
11222211
knkkkkk
k
nnnn
Laaa
knkAA
A
kkkkPPkkk
UAPLPLPLPL
唯一确定了（
部分，即的下三角的每一输出的矩阵；的上三角为输出的矩阵；唯一确定了
2,LU分解算法
e nd
e nd
e nd
jkakiajiajia
:nkjf or
kkakiakia
nkif or
tankkkka
nkkkkankka
a( k,k,n) ;ta
mkkkk
knkaabsmq
nkf or
as i z enm
aL G aus sakkf unc t i on
);,(),(),(),(
1
);,(/),(),(
:1;):),((
);:),(():,(;1)(
) ) ) ;,:((m a x (],[
1:1
);(],[
)(],[
运算量和高斯消去法一样
e nd
e nd
t;,k ) a( k k ( i );a( k k ( i ),k ) a( i,k )
a( i,k ) ; t
,n-k i f or
,n-f or k
e nd
t;k ( k ),,) p( k;p( k k ( k ),,)p( k,,)
kpt
nkf or
neyep
aL G aus sakk
as i z enm
aj i es anj i aof e npulf unc t i on
11
21
,) ;,(
1:1
);(
);(],[
);(],[
)(],,[
e nd
e nd
a( i,j ) ; u( i,j )
i,n j f or
:nf o r i
e nd
e n d
e nd
a( i,j ) ;i,j ) l (
e l s e;i,j ) l (
j i f i
:i j f or
:nf o r i
1
1
1
1
3、三对角矩阵方程的追赶法三对角且主对角严格占优矩阵方程是一类来源丰富的问题。比如，微分方程数值解或样条插值等问题中的正规方程组。解这种问题必须考虑其矩阵稀疏的特征，减少算法的计算量。
nn
n
ab
c
ab
ca
T
1
22
11
三对角矩阵形如：
fTx?问题：求解
T的 LU分解具有形式：
n
n
n u
c
u
cu
U
l
l
L
1
2
11
2
,
1
1
1
由 T=LU推得：
).:2(,/ ; 1111 niclauublau iiiiiii
1:1 ; /)(,/
:2 ;,
1
111
niuxcyxuyx
niylfyfy
yUx
fLy
fTx
iiiiinnn
iiii
e n d
1) ) / a ( i) ;f ( i*c ( i)-( f ( i)f ( i)
1:-1:1-nif or
f ( n ) / a ( n ) ;f ( n )
e n d
1) ;-f ( i*b( i)-f ( i)f ( i)
n:2if or
e n d
1) ;-c ( i*b( i)-a ( i)a ( i)
1) ;-b( i) / a ( ib( i)
n:2if or
le n g t h ( a ) ;n
f)c,b,a,z h u ig a n f a (ff u n c t io n
最终解为 f.
乘除运算量：
5n=O(n).
4、对称正定矩阵方程的 cholesky分解法对称正定矩阵方程的来源比较丰富，比如线性回归、拟合等问题。解这类问题必须考虑其矩阵对称的特征，减少算法的计算量。
为对称正定矩阵。其中问题：求解 A,bAx?
由于 A为对称正定矩阵,A必有 cholesky分解：
nnn
T
d
d
d
ll
l
L
L D LA
2
1
21
21
D,
1
1
1
,其中
j
j
k
jkikkijij
j
k
kjkjjj
ii
T
dlldal
dlad
nidalad
LD LA
/)(
n:1jin,:2j
:2,/,
1
1
1
1
2
111111
得到关系：
的两边，比较
e n d
e n d
j ) ;t ) / a ( j,-j)( a ( i,j)a ( i,
e n d
k) ;a ( j,*k)a ( i,*k)a ( k,t t
1-j:1kf or
0; t
n:1jif or
t;-j)a ( j,j)a ( j,
e n d
);,(*2)^,(t t
1-j:1kf or
0; t
n:2jf or
e n d
);a ( i,1) / d( 1a ( i,1)
n:2if or
s i z e ( a ) ;n][ m,
)c h ol e s ky ( aD][ L,
kkakja
f unc t i on
高斯消去法的一半。
，是很大时，相当于当运算量：乘除为上三角矩阵说明：
命令：
3
23
6
1
)29(
6
1
L,
)(
nn
nnn
LLA
Ac holLM A T L A B
T
三、线性代数方程组的迭代法线性代数方程组的迭代法并不适用于所有问题，但它对一些特定类型的问题非常有效。当问题是大型稀疏矩阵方程时，高斯消去法的效率会变得非常低，而且有时还会超出内存要求。对于这样的问题就需要使用迭代法。
大系统问题的求解最终归结为大型稀疏矩阵方程。比如，
电网络、场方程的数值计算、运筹问题等。
尽管迭代法的种类很多，这里只介绍其中的三种：
1,Jacobi迭代；
2,Gauss-Seidel迭代 ;
3、超松弛迭代法（ SOR)。
.x
,3;,2,1,0,2;1
)(*
*)(
)()1(
)0(
的精度为近似解为方程组的真解，其中收敛条件：、
迭代：、
、选定初始值：
算法：
k
k
kk
n
x
xx
kdBxx
Rx
dBxxbAx
1,线性代数方程组的迭代法的一般理论迭代算法简单，但问题是：
1、能否保证算法收敛？
2、能否充分利用矩阵的稀疏性，使运算量和存贮量尽量少。
[定理一 ] 迭代法收敛的充分必要条件是 0lim
1
k
k B[定理二 ] 迭代法收敛的充分必要条件是
1)(?B?
[定理三 ] 迭代法收敛的 充分条件 是
)1()(*( k )
)0(*( k )
1
x
*
1
x,1
kk
k
xx
G
G
x
xx
G
G
xB
和且谱范数行和范数列和范数
)(B
m a x; m a x
2
1
1
1
11
T
n
j
ij
ni
n
i
ij
nj
BB
bBbB
收敛性定理对于任意的初值 nRx?)0(
迭代矩阵 B的构造充分利用矩阵的稀疏性，使运算量和存贮量尽量少的办法就是要求迭代矩阵 B与原矩阵 A有相同的稀疏结构。具体就是：
0
a
0
aa0
U,
0aa
0a
0
L
);a,,a,d i a g ( aDU,LDA
n1,-n
1n12
1-nnn1
21
nn2211
)( bxULDbAx
n,1i,)x-(1) / axaa-(b
b)x )D 1-(-U( L ) x-(D
0,)
1
1(
1
SO R 3
n:1i,) / axaa-(b
) Se i d e l-G a u s s2
n,1i,) / axa-(bx
)( 1
( k )
ii
( k )
j
n
1ij
ij
)1(
1-i
1j
iji
)1(
( k )1)(k
ii
( k )
j
n
1ij
ij
)1(
1-i
1j
iji
)1(
)()1(
ii
( k )
j
n
ij1,j
iji
1)(k
i
)()1(
k
j
k
i
k
j
k
i
kk
kk
xx
DDD
xx
bUxxLD
bxULDxJ a c o b i
松弛因子。迭代：、
（迭代：、
迭代：、
常用迭代法及其收敛性
nRx )0(
定理 4,SOR收敛的必要条件是 0<ω<2.
定理 5：如果 A是严格对角占优矩阵或不可分弱对角占优矩阵，则
1,Jacobi收敛；
2,gauss-seidel收敛；
3、当 0<ω<=1 时,SOR必收敛。
定理 6：如果 A是对称正定矩阵，则
1、当 2D-A正定时,Jacobi收敛；
2,gauss-seidel收敛；
3、当 0<ω<2时,SOR必收敛。
其他矩阵不能保证收敛四、线性代数方程组共轭梯度法
)(m i n
0)()(
2
1
f ( x )
*
***
*
xfbAx
bAxxfxf
xA
xbAxx
TT
，因此由于
。唯一的极小点是对称正定矩阵时，有当二次函数：
即，解对称正定矩阵方程组的问题等价于求二次函数极小点的问题。
如何寻找 f(x)的极小点呢？基本思想就是 下降法 。
1、下降法的基本原理是：的极小值。具体的算法到以“最快的速度”下降使逐步产生一串点出发从某一点
)(
)()()f ( x
,,,,,,
)()1(( 0 )
)()1()0()0(
xf
xfxf
xxxx
k
k
为止。
或
、直到即的极小点上求、在射线；下降的方向处确定一个使、在
。现以求得
x
)f ( x 3
);(m i n)f ( x
,)(2
)(1
)(1)(k
( k )
)()(
0
1)(k
)1()()(
)()(
)(
k
kk
t
kkk
kk
k
x
tpxf
xxftpxx
pxfx
x
2、最速下降法我们知道，函数的负梯度方向是函数值下降最快的方向。
称为最速下降法作为下降方向的算法取 ( 1 ) )( )()()(?kkk Axbxfp
b)( A x)(ptAp)(p
pb)A ( p)tp(x
2
1
)tpA ( x)(p
2
1
)tpf ( x
dt
d
( k )T( k )( k )T( k )
( k )T( k )T( k )( k )( k )( k )T( k )
( k )( k )
( 3 ),x
( 2 )
)(
)(
,0)(
)()(1)(k
)()(
)()(
)()(
k
k
k
kTk
kTk
k
kk
ptx
App
pp
ttpxf
dt
d
得到时即当
.;
)(
)(;,,;
.0,1,2,k,
)()()1(
)()(
)()(
)(
)()(
)0(
k
k
kk
kTk
kTk
k
k
kk
n
rtxx
Arr
rr
t
e n db r e a krif
Axbr
Rx
最速下降算法：
function x=zuisu(a,b,x,eyta)
n=length(b);
r=b-a*x;
while max(abs(r))>eyta
t=r'*r/(r'*a*r);
x=x+t*r;
r=b-a*x;
end
这一算法有简单易行，可充分利用 A的稀疏性等特点，
但当 A的最大特征值远远大于最小特征值时收敛速度变得非常慢，以至于完全不适用。最速下降法并非最速！下面的共轭梯度法可有效解决这一问题。
3、共轭梯度法能否选出比负梯度方向更好的下降方向吗？能！
.0r4;1,,2,1,0,0)( 3
,,,,p2
,,,,1;
)(
)(
,;
)(
)(
,;,2,1,0,;,,
1)-(n
)()(
)()1(( 0 )
)()1()0(
)()(
)1()(
1
)(
1
)1()1(
)1()1(
)()(
)()(
)()()1(
)()(
)0()0()0()0()0(
、
、
共轭的；是两两关于、
是两两正交的；、可以证明：
其中；其中以求得，下降方向
kirp
App
rrr
App
rp
spsrp
Axbr
App
pp
tptxx
kpx
rpAxbrRx
kTi
k
k
kTk
kTk
k
k
k
kk
kk
kTk
kTk
k
k
k
kk
kk
n
).1(
10
*)(
)1(
nkxx
nr
k
n 步就能得到精确解计算说明上述计算过程最多这一算法称为共轭梯度法，简称 CG(conjugate gradient)
实际计算中由于误差的影响,)(kr 之间的正交性很快就会消失，以至于有限步内不能得到精确解，所以 CG实际是一种迭代算法。
CG是保稀疏且便于并行处理的算法，它是求解对称正定矩阵方程最优秀的算法之一。而且目前它已推广到一般矩阵方程，称为预优共轭梯度算法。
实用共轭梯度算法如下：
function x=cg(a,b,x,eyta)
n=length(b);
r=b-a*x;
p=r;
q0=r'*r;
while q0>eyta
w=a*p;
t=q0/(p'*w);
x=x+t*p;
r=r-t*w;
q=r'*r;
s=q/q0;
p=r+s*p;
q0=q;
end
运算量：每一步迭代的乘除运算量为
nn 22?
§ 2.6 矩阵特征值问题
).,(0)(,xxIARA nn 的解求问题,
一、幂法
MATLAB命令,D=eig(A)和 [X,D]=eig(A)
.,,,
,0,
21
21
n
n
xxx
nA
线性无关的特征向量个并且对应的特征值设
)1(
111
1
11111
)(
1111
)(
)1()()0()0(
)()(;0),(
,2,1,;0,
kkkk
kk
kkn
vxxv
xvk
kAvvvRv
显然常数充分大后，则当；?
.
),
(k
12
)(
)1()1(
)1(
)(
1
收敛因子值和特征向量。的按模最大的近似特征分别为分量的按模最大的是充分大时，结论：当
Av
vv
v
v
k
kk
lk
l
k
l
实用算法
：
时终止。当；
)1()(
)()(
)(
)(
)1()(
)0(;/
*))((
));(m a x (],[;
,2,1
];1,,1,1[
kk
k
kk
k
k
k
k
k
kk
vv
mvv
mlvs i gnm
vabslm
Avv
k
v
function [lmta,x]=mifa(a,epsl)
[m,n]=size(a);
x0=ones(n,1);
while 1
x=a*x0;
[lmta,m]=max(abs(x));
lmta=sign(x(m))*lmta;
x=x/lmta;
if abs(x-x0)<epsl,break,end;
x0=x;
end
一、幂法二、反幂法及其原点位移反幂法用来求 A的按模最小的特征值。
.,,,
,0,
21
121
n
nn
xxx
nA
线性无关的特征向量个并且对应的特征值设
n
k
k
nn
k
k
nn
m
xxv
A
xxAxAx
1
l i m
)m a x (/l i m
.
,
1
)(
1
1
1
幂法即可！收敛因子执行上述所以只要对因时终止。当；
)1()(
)()(
)(
)(
)1(1)(
)0(;/
*))((
));(m a x (],[;
,2,1
];1,,1,1[
kk
k
kk
k
k
k
k
k
kk
vv
mvv
mlvs i gnm
vabslm
vAv
k
v
实用算法
：
function [lmta,x]=fanmifa(a,epsl)
[m,n]=size(a);
x0=ones(n,1);
while 1
x=a\x0;
[lmta,m]=max(abs(x));
lmta=sign(x(m))*lmta;
x=x/lmta;
if max(abs(x-x0))<epsl,break,end;
x0=x;
end
lmta=1/lmta;
1、反幂法
2、带原点位移的反幂法反幂法与“原点位移”相配合，求指定点附近的某个特征值和特征向量，并可用于加速幂法的收敛性。
.
1
)(
sA,,0
)()(
,
s
m
sIA
liss
xsxsIAxAx
AsIAs
k
l
lli
用反幂法求得特征值，可对最近的的离点是即可知，若的原点位移。由为称对于给定的常数另一方面,s离 A的某个特征值 λ越近，收敛越快。因此不论用幂法求 A的按模最大特征值，还是利用反幂法求 A的按模最小特征值，为了加快收敛，均可以用迭代 m步后的近似值 lmta作为最初始的位移值,实行 动态位移迭代 。
动态位移幂法
function [lmta,x]=dongtamifa(a,epsl0,epsl)
[lmta,x]=mifa(a,epsl0);x0=x;lmta0=lmta;
n=length(x0);
while 1
x=(a-lmta0*eye(n))\x0;
[lmta,m]=max(abs(x));
lmta=sign(x(m))*lmta;
x=x/lmta;
lmta=1/lmta+lmta0;
if max(abs(x-x0))<epsl,break,end;
x0=x;
lmta0=lmta
end
动态位移反幂法
function [lmta,x]=dongtaifanmifa(a,epsl0,epsl)
[lmta,x]=fanmifa(a,epsl0);x0=x;lmta0=lmta;
n=length(x0);
while 1
x=(a-lmta0*eye(n))\x0;
[lmta,m]=max(abs(x));
lmta=sign(x(m))*lmta;
x=x/lmta;
lmta=1/lmta+lmta0;
if max(abs(x-x0))<epsl,break,end;
x0=x;
lmta0=lmta;
end
三、幂法综述
1、幂法和反幂法只能用于求解 可对角化矩阵的实数特征值和特征向量，不能求解复数特征值 ；
2、幂法和反幂法均为线性收敛，收敛速度由收敛因子决定，效率不高。
3,动态位移可以大大减小收敛因子加速收敛 ；
4,不适于求解全部特征值。
5、对称矩阵自然适用于幂法，此时采用 2-范数标准化向量的算法至少平方收敛。
求解一般矩阵全部特征值的办法是下面的相似变换法四、矩阵全部特征值的 QR迭代算法求解一般矩阵全部特征值的办法是相似变换法。其理论根据是：
。均为一阶或二阶实方阵其中使正交矩阵
ii
ss
s
s
nn
A
A
AA
AAA
TAPP
PRA
,
,,
222
11211
1
.l i m
,,)Q (A
).(],[
,
.,2,1,A,A,
)(
m
)()1)()()(1-( m )1)(m)(
)()(1)(m)()(( m ))0(
TA
AAARQQ
AqrRQM A TL A B
RQQRA
mQRRQAA
m
mmmmmm
mmmm
且相似于即命令：交三角化分解。
解为正为上三角矩阵，称此分为正交矩阵，其中，
实现过程：
（
1、矩阵的两种正交变换平面旋转变换变换。为平面旋转矩阵，或则矩阵记
G i v e n s
k i
1
cs-
sc
1
)k,J ( i,
,1,s i n,c o s
nn
R
nkisc
kijxy
y
xx
xx
x
s
xx
x
c
kijxy
cxsxy
sxcx
JxyRx
IJJkiJJ
jj
k
ki
ki
k
ki
i
jj
kik
ki
n
T
,,
0
y
,,0 x
,,
y
,
.),,(
22
i
2222
k
i
则时，取一个非零元素。即它的意义是消去向量的的分量表示正交矩阵，即?
镜面反射矩阵
.)0,,0,1(.
2
),(
2
1
),,,(
s i g n ( x,/) ),(m a x (
,,,),,,(2
,1
2,1,
11
1
22
1
2
2
211
2
1
221
1
2
T
T
T
T
n
n
T
n
T
Tn
eeHxH
uuI
u
u
u
u
IH
txu
xxxexu
xtxxxa bst
xxxxxx
HHH
H
rH o us e ho l d e
IHR
使确定的镜面反射矩阵
）
不全为零，则由）（
对称正交；）（
具有下列性质：矩阵变换矩阵。或称为镜面反射矩阵，矩阵镜面反射矩阵的意义是“成批”消去向量的非零元素。
function [a,b,x]=householder(x)
% x gived,seek a househouder convert H=I-bxx',make Hx=-a.
n=length(x);
t=max(abs(x));
if t==0
a=0;
b=0;
else
x=x/t;
a=sqrt(x'*x);
if x(1)<0,a=-a;end
x(1)=x(1)+a;
b=1/(a*x(1));
a=t*a;
end
2、方阵的正交三角化分解方阵的正交三角化分解,即为上三角矩阵。为正交矩阵，其中 RQQRA,?
;*
**;~
)1(;
*;
~
) ) ;1(,,(],,[;
*;) ) ;1(,,(],,[
2
2
1
12
2
2
2
2
12212
1
1
111
A
a
a
AHH
H
e y e
H
A
a
AHb x xIHArh ou s e ho l d exba
A
a
AHb x xIHArh ou s e ho l d exba
T
T
分解方法：
function [arfa,bata,a]=qrfenjie(a)
[m,n]=size(a);
for k=1:n-1
[arfa(k),bata(k),a(k:n,k)]=householder(a(k:n,k));
for j=k+1:n
a(k:n,j)=a(k:n,j)-bata(k)*a(k:n,k)'*a(k:n,j)*a(k:n,k);
end
end
function [q,r]=qrgouzao(a)
%实现 a=qr,q是正交矩阵,r是上三角矩阵。
[m,n]=size(a);
[arfa,bata,a]=qrfenjie(a);
q=eye(n)-bata(1)*a(1:n,1)*a(1:n,1)';
for k=2:n-1
q(1:k-1,k:n)=q(1:k-1,k:n)-bata(k)*q(1:k-1,k:n)*a(k:n,k)*a(k:n,k)';
q(k:n,k:n)=q(k:n,k:n)-bata(k)*q(k:n,k:n)*a(k:n,k)*a(k:n,k)';
end
for i=1:n-1
for j=i+1:n
r(i,j)=a(i,j);
end
end
for i=1:n-1
r(i,i)=-arfa(i);
end
r(n,n)=a(n,n);
3、化矩阵为 Hessenberg型矩阵。不可约时矩阵。所有称为（上）
形如：
H e s s e nb e r g
hH e s s e nb e r g
hh
hhh
hhhh
hhhh
H
ii
nnnn
nn
nn
nn
,0
,
,1
1,
31,332
21,22221
11,11211
矩阵。为使得正交矩阵定理：
H e s s e n b e r g
AQQBQRA Tnn,,
function [q,h]=hessenberghua(a);
% h=q'aq,h为 Hessenberg矩阵,q为正交矩阵。
[m,n]=size(a);
q=eye(n);
for k=1:n-2
[arfa(k),bata(k),a(k+1:n,k)]=householder(a(k+1:n,k));
p=eye(n-k)-bata(k)*a(k+1:n,k)*a(k+1:n,k)';
a(1:k,k+1:n)=a(1:k,k+1:n)*p;
a(k+1:n,k+1:n)=p*a(k+1:n,k+1:n)*p;
q(1:n,k+1:n)=q(1:n,k+1:n)*p;
end
h=a;
for i=1:n-2
h(i+1,i)=-arfa(i);
end
for j=1:n-2
for i=j+2:n
h(i,j)=0;
end
end
注：当 A为对称矩阵时，其
Hessenberg 形为对称三对角矩阵。
4,Hessenberg矩阵在 QR分解下的不变性变换下的不变性。矩阵这种性质称为变换，的到矩阵。称此为由均为和下，
正交相似变换，、
分解，、
矩阵，则在为设
QRH e s s e n be r g
QRHHH e s s e n be r gHQ
RQHQQH
QRH
H e s s e n be r gH
T
~~
~
2
QR 1
显然，实现 Hessenberg矩阵的 QR分解用 Givens变换合适。即执行 n-1次 Givens变换,
运算量。一过程需完成这然后计算为上三角形使
2
121
121
4
.
~
,
,1:1),,1,(
n
JJRJH
RHJJJniiiJJ
T
n
TT
nnii
function [q,h]=givensqr(h)
%h为不可约 hessenberg矩阵，用 givens变换进行 QR分解。
[m,n]=size(h);
q=eye(n);
for k=1:n-1
if h(k+1,k)~=0
c=h(k,k)/sqrt(h(k,k)^2+h(k+1,k)^2);
s=h(k+1,k)/sqrt(h(k,k)^2+h(k+1,k)^2);
d=[c,s;-s,c];
h(k:k+1,k:n)=d*h(k:k+1,k:n);
if k==1
q=[d' zeros(2,n-2);zeros(2,n-2)' eye(n-2)];
else
q=q*[eye(k-1) zeros(k-1,2) zeros(k-1,n-k-1)...;zeros(k-1,2)' d' zeros(2,n-k-1)...;zeros(k-1,n-k-1)' zeros(2,n-k-1)' eye(n-k-1)];
end
end
end
5、基本 QR算法
QR( 2 ) A
)(A 1)
,2,1,0,2
)(],[P
1
kk1k
121k
00
0
T
n
TT
kkk
nn
JJJQRQ
kQR
Ahu aH e s s e n be r gA
AH e s s e n be r gA
RA
注：（
迭代、
形为、约化算法：
一般 QR 迭代算法的收敛性比较复杂，这里不再介绍。仅指出，在一定条件下由基本 QR算法生成的系列
kA 收敛为准上三角矩阵，对角块按特征值的模从大到小排列。
function [q,h]=jibenQR(a,epsl)
[m,n]=size(a);
[q0,h]=hessenberghua(a);
t=zeros(1,n-1);
while 1
for i=1:n-1
if abs(h(i+1,i))<=(abs(h(i,i))+abs(h(i+1,i+1)))*epsl
h(i+1,i)=0;
t(i)=i+1;
end
end
k=1;
while k<=n-2
if t(k)==t(k+1)
yes=1;
k=n-1;
else
k=k+1;
yes=0;
end
end
if yes==0,break,end;
[q,h]=givensqr(h);
h=h*q;
q=q0*q;
q0=q;
end
6,QR 算法的改善基本 QR算法计算量和存储量都很大，且若收敛是线性的。
因此，需要从这两方面加以改善。
（ 1）划分和收缩
（ 2）原点位移第三章 数值逼近内容：
数值逼近也称数据模拟，是为离散数据（不论何种途径得到）建立简单连续模型的方法。包括：
1、插值方法
2、数据拟合方法
3、最佳逼近要求：
1、了解 MATLAB功能函数
2、掌握数据模拟的各种方法及应用
§ 3.1 问题的提出例一,在某化学反应里，侧的生成物的质量浓度 y与时间
t(min)的关系入表。为了研究该化学反应的性质，如反映速率等,欲求 y与 t之间的关系 y=f(t)。
1 2 3 4 6 8 10 12 14 16
4.00 6.41 8.01 8.79 9.53 9.86 10.33 10.42 0.53 10.61
例二,逢山开道（山区修建公路的路线选择问题 1994）
在某山区修建公路,已知该山区的地形高度见表一 （ 单位
m） 。 雨季在山谷形成一溪流,雨量最大时溪流水面宽度 W
与 x（ 溪流最深处的 ） 坐标的关系可以表示为
4 00 02 40 0,522 40 0
4/3
xxxw
要求：从山脚开始经居民点至矿区修一条公路 （ 一般公路,
桥,隧道 ） 给出路线设计方案使工程成本最小 。
3200 1430 1450 1460 1550 1500 1600
2800 950 1190 1370 1500 1200 1100
2400 910 1090 1270 1500 1200 1100
2000 880 1060 1230 1390 1500 1500
1600 830 980 1180 1320 1450 1420
1200 740 880 1080 1130 1250 1280
800 *650 760 880 970 1020 1050
400 510 620 730 800 850 870
0 370 470 550 660 670 690
Y/x 0 400 800 1200 1600 2000
部分数据表：
工程种类 一般路段 桥梁 隧道工程成本
（ 元 /m）
300 1500 （ 长度 <=300m ）,3000 （ 长度
>300m）
α=上升高度 /
水平距离
<=0.125 =0 <=0.100
2000
工程要求：
例三,水道测量模型 （ 1989.美国 ）
问题 水平方向的坐标 x,y以 Td（ =0.914m） 为单位,水深方向 Z 以 Ft(=30.48cm)为单位,下表给出了水面直角坐标 (x,y)处的水深 Z,这是在低潮时测得的 。 如果船的吃水深度为 5Ft,试问在矩形城 (75,200)(-50,50)中行船应避免进入那些区域?
工作,1,要根据数据做出 ① 地貌图 （ 三维 ）,② 等高线 （ 等势图 ） 二维 。 2,由坡度限制,沿给定的网格点选择路线
x 129.0 140.0 108.5 88.0 185.5 195.5 105.5
y 7.5 141.5 28.0 147.0 22.5 137.5 85.5
z 4 8 6 8 6 8 8
157.5 107.5 77.0 81.0 162.0 162.0 117.5
-6.5 -81.0 3.0 56.5 -66.5 84.0 -38.5
9 9 8 8 9 4 9
x
y
z
工作：将海底曲面图绘制出来 。
总之：根据给定采样数据而定出其实际分布图,这就是数值模拟问题,采用的方法就是插值和拟合 。
§ 3.2 数值模拟概论在科学与工程等实际问题中，实际数学模型或者难于建立或者难于求解，但其数据模型（由实验或测量所得到的一批离散数据）容易得到。能否通过处理这些数据来建立实际模型呢？ 这里我们仅以一维问题来说明。
给定,y=f(x)数据
n10 x xx?
n10 y y y?
通过这批数据希望找到对象 y=f(x)的更多的信息（或全部信息）。通常只能找到 f(x)的近似表达式 φ(x),φ(x)是根据研究对象的特征确定的。
0
s i nc o s
:
n
nnnn xbxax
，则研究对象是周期振动的比如对象是指数增（减）并最终稳定到某个值,则
xeccx 10
对象是直线运动、匀加速运动，则
xaax 10
2210 xaxaax
一般地，我们都要通过对研究对象实际的了解，做以下工作：
1、确定它的基本特征函数
);(,),(,21 xxx m
2、对这些特征函数（基函数）进行一种恰当的构造，得到
f(x)近似的数学模型,
非线性的,.,,,,,,; 21 mcccx?
最常用也就是最简单的一种构造方式就是,
线性
1
xcx jj
m
j
.)(,,,21 的线性组合即 xm
3、按某种寻优策略，由已知数据确定未知参数
mccc,.,,,,,,21寻优策略的不同，产生了不同的数值逼近方法
§ 3.3 插值方法寻优策略 就是 过点,即
niyx jj,,2,1,
根据需要可附加补充原则,
np xxcx,0
nxx,0
1,p阶光滑度
2、边界条件,处光滑性,周期性等
（ 整体） 误差,
),(,0 nxxxxxfxR
插值方法中最常用的一类就是 代数插值 。
代数插值：即多项式的插值 m
m xaxaax10?
一、代数插值
（一）拉格朗日插值公式（ Lagrange)
n10 x xx?
n10 y y y?
给定,y=f(x)数据确定次数不超过 n的多项式 P(x),使
niyx ii,,2,1,0,( ）
).(
)!1(
)(
)()()(
)()(
L a g r a n g e;,,2,1,0,
)()(
)(
)(;)()(
)1(
0
0
x
n
f
xxfxR
xlyx
nj
xxx
x
xl
xxx
n
n
j
jj
jj
j
n
j
j
整体误差：
插值公式为：则令
拉格朗日插值公式的讨论
。其中 ))((m a x
8
)(;)(
],[
2
2
01
2
1
01
0
0
10
1
1
10
xfM
xx
M
xR
y
xx
xx
y
xx
xx
x
xxx
线性插值公式（ n=1的情况）
显然，两插值节点越近，越精确。
2、二次插值（抛物线插值,n=2的情况）
).)()((
!3
)(
)(;
))((
))((
))((
))((
))((
))((
)(
2102
2
1202
10
1
2101
20
0
2010
21
2
xxxxxx
f
xR
y
xxxx
xxxx
y
xxxx
xxxx
y
xxxx
xxxx
x
n越大精度越高！是这样吗？下面我们看一实例。
实例演示,分别用 n=1,2,4,8,16,32时的拉格朗日插值逼近
f(x)。
]1,1[,251 1)( 2 xxxf
function f=lagelangri(t,y,x)
n=length(t);
for j=1:n
l(j)=1;
for i=1:n
if i~=j
l(j)=l(j)*(x-t(i))/(t(j)-t(i));
end
end
end
f=y*l';
定义 lagrange函数
function f=shiyanhs(x)
f=1./(1+25*x.^2); 定义实验函数
function g=lagelangritu(a,b,n)
h=(b-a)/n;
t=a:h:b;
f0=shiyanhs(t);
x=a:0.02:b;
m=length(x);
for i=1:m
p(i)=lagelangri(t,f0,x(i));
end
f=shiyanhs(x);
plot(x,f,'r',x,p);
gtext('f(x)');
gtext('p(x,n)');
axis([-1,1,-1,1]);
画图,f(x) 和 p(x)
总之，我们看到利用拉格朗日公式逼近函数 f(x),n小不行，
n大也不行。这是为什么呢？
，看出从误差：
)())(()!1( )()( 10
)1(
n
n
xxxxxxnfxR
当然大。内取的插值节点少，小代表在区间、
而增加。
的增加附近的值随、在边界的均值附近比较小，而的值，在（、
跑得快！（
比的增加呈指数级增长，的值，常常随、
)(],[3
n
,,,)())(()2
)!1
)(1
0
1010
)1(
xRban
xx
xxxxxxxxxx
n
nxf
n
nn
n
差。计算量和计算稳定性变的增加，随着
（（（（（
（（（（（
插值公式来看：、从
)(
)()(;,,2,1,0
,
)))))
)))))
)(
4
0
1100
1100
xln
xlyx
nj
xxxxxxxxxx
xxxxxxxxxx
xl
Lag r an ge
j
n
j
jj
njjjjjjj
njj
j
综合以上分析我们可以得到如下结论：
在小区间上应用低次拉格朗日差值公式（ n<=6)有效。
问题是对大的区间怎么办？那就是用 分段低次插值。
（二）分段低次插值分段线性插值,
n10 x xx?
n10 y y y?
给定,y=f(x)数据确定 P(x),使
niyx ii,,2,1,0,(2 ）、
上是线性函数。在每一小区间,],[)(1 1 ii xxxP?
根据要求立即得到：
即为所求！
.
)(
.,,2,1],,[
1
1
1
1
1
1
1
1
iii
i
i
i
i
i
i
i
ii
i
i
ii
i
ii
xxh
y
h
xx
y
h
xx
y
xx
xx
y
xx
xx
xP
nixxx?
若用基函数的形式表示，则：
0
)(
11
10
10
1
0,令
xxx
xxx
xx
xx
xN
即为所求！);(yP ( x )
0
)(N
1.-n,1,2,k
]x,[x x 0
)(N
n
0k
k
10
1
1
1
n
1k1-k
1
1
1
1
1
1
k
xN
xxx
xxx
xx
xx
x
xxx
xx
xx
xxx
xx
xx
x
k
n
nn
nn
n
kk
kk
k
kk
kk
k
function f=fenduanxianxing(t,y,x)
% x在 [t(1),t(n)]
n=length(t);
i=1;
while 1
if (x>=t(i))&(x<=t(i+1)),break,end;
i=i+1;
end
f=y(i)*(x-t(i+1))/(t(i)-t(i+1))+y(i+1)*(x-t(i))/(t(i+1)-t(i));
定义分段线性插值函数
]1,1[,251 1)( 2 xxxf
实验：用分段线性插值逼近函数
function g=fenduanxianxingtu(a,b,n)
h=(b-a)/n;
t=a:h:b;
f0=shiyanhs(t);
x=a:0.02:b;
m=length(x);
for i=1:m
f(i)=fenduanxianxing(t,f0,x(i));
end
f1=shiyanhs(x);
plot(x,f1,'.-r',x,f);
gtext('f(x)');
gtext('p(x,n)');
axis([-1,1,-1,1]);
g=fenduanxianxingtu(-1,1,50)
计算：
分段线性插值的讨论：
误差分析
2
0
)2(
],[
)2(
],[
2
1
8
R ( x )
],,[
,)(m a x,)m a x (
)(m a x,
8
R ( x )
],[
0
1
h
M
xxx
xfMhh
xfMh
M
xxx
n
xxx
i
i
xxx
ii
i
ii
n
ii
有则令时，有更多的插值节点。、逼近精度提高，需要只是连续，光滑性差；、
缺点：
好。计算简单，数值稳定性、
时，、
优点：
2
)(1
)(2
);()(01
xP
xP
xfxPh
下面的实例演示可以说明
g=fenduanxianxingtu(-1,1,n)
分别用 n=4,8,16,32的分段线性差值逼近函数（ h=2/n)
]1,1[,251 1)( 2 xxxf
样条函数与样条插值,
n10 x xx?
n10 y y y?
给定,y=f(x)数据确定 s(x),使
niyxs
xxCxs
ii
n
k
,,2,1,0,(3
),()(2 01
）、
、
)2
],[)(1 1
kk
xxxs ii
的多项式。（
上是次数不超过在每一小区间、
若 s(x)存在，则称 s(x)为对应于给定数据的一个 k次样条插值。满足条件 1和 2的函数 S(x)称为 k次样条函数。
样条函数的这种数学定义来源于 绘图员沿着受压铁约束的样条（一根具有很好弹性的木条），画出各种光滑曲线（样条曲线）。这种样条曲线可以看成弹性细梁受集中载荷作用而生成的 挠度曲线。
S(x)的存在性讨论：
).()(
.,2,1,0
0 0
0 !/
)(
1
m
RCx
m
x
xmx
xx
m
m
m
m
m
次半截幂函数定义：
是待定常数。其中次样条函数可表示为：
jj
n
j
k
jj
k
j
j
jk
ba
kxxbjxaxs
k
,
!)(!)(
1
10
实际中常用 k=2,3的情况，即二次样条和三次样条个待定常数。共有
：二次样条函数可表示为
2
!2)(
!2
)(
1
1
2
2
2102
n
xxb
x
axaaxs
n
j
jj
个待定常数。共有
：三次样条函数可表示为
3
!3)(!)(
1
1
3
3
0
3
n
xxbjxaxs
n
j
jj
j
j
j
二次样条插值公式因为二次样条函数的确定需要 n+2个条件，而
niyxs ii,0 )(2
给定了 n+1个，所以还需补充一个条件（自然是边界条件）。通常有
).()(,)()( 2
))(()(1
202
'
20
'
2
''
2
01
01'
00
'
2
nn
nn
xsxsxsxs
yxs
xx
yy
yxs
、周期性条件或、
补充条件类型：
补充条件 1下的二次样条插值公式
( 2 ) /)(
y)x-(x
)(],,[
( 1 ) 1:0
,/)(,
)(,)(
0
'
0
0
01
0
'
01000
'
00
'
210
11i
1122
hy
h
yy
c
c
yxsxxx
ni
xxhhyyy
yxsyxs
iiiiiiii
iiii
得到由得到由
1:1),0()0(),()(2
1:0),)(()()(
],,[ ],[)(1
'
2
'
210
1
2
12
112
nixsxsxxCxs
nixxxxcxxxs
xxxxxxs
ii
iiiiii
iiii
、
内为二次多项式在、
采用如下方法来确定：
( 3) 1:1
/)(
1:1),0(s)0(s
1
11
1
1
111-i
'
2
'
2
ni
h
h
yy
h
yy
c
h
h
c
hchc
nixx
i
i
ii
i
ii
i
i
i
i
iiiii
ii
得到由联立式（ 1）、（ 2）、（ 3）确定了二次样条插值公式：
))(()()(
],[,],,[;1:0,
12
10
i
iiiiii
iin
ii
xxxxcxxxs
xxxixxx
nic
第三步：；即确定第二步：
、、第一步：计算算法个待定常数。的确定立求解：当然，我们也可采用联
2
!2)(
!2
)(
)0(
:0,)(
1
1
2
2
2102
'
00
'
2
2
n
xxb
x
axaaxs
yxs
niyxs
n
j
jj
ii
这样不如前面的方法简单！
实验：用二次样条插值逼近 f(x),h=(b-a)/n.比较 n=4,8,16,32的效果
]1,1[,251 1)( 2 xxxf
二次样条插值的误差分析
).(
12
5
3
4
),m a x (
,
)()()(
1],,[)(
)4(
2
3
2
4
abffchh
chR
chR
chR
xsxfxR
baCxf
i
其中及其导数的界的估计为的余项的二次样条插值则补充条件定理：设三次样条插值公式因为三次样条函数的确定需要 n+3个条件，而
niyxs ii,0 )(3
给定了 n+1个，所以还需补充二个条件（自然是边界条件）。通常有
).0()0(
)0()0(
)()( 3
)0()0(2
)0()0(,1
202
'
20
'
2
202
3003
''
3
'
00
'
3
n
n
n
nn
nn
xsxs
xsxs
xsxs
yxsyxs
yxsyxs
）（自然有：周期性条件问题
，：问题
，问题补充条件类型：
解由插值条件和补充条件构成的线性方程组，可以得到三次样条插值公式。但这样比较复杂！
下面只就问题 2来讨论，并采用另一种途径称之为三 弯矩法。
)(,)(;)(
],,[
.:0,)(],[)(
),()(],[)(1
1133
1
1
3
1
313
10
2
313
得并使其满足对此式做两次不定积分则内是线性函数。设在故
，内为三次多项式且在、
iiii
i
i
i
i
i
i
ii
iiii
ii
yxsyxs
h
xx
m
h
xx
mxs
xxx
nimxsxxxs
xxCxsxxxs
( 1 )
h
)
6
m
-(
h-
)
6
m
-(
6
)(
6
)(
)(
i
2
1i
1
i
1
2
i
3
1
3
1
3
ii
i
ii
i
i
i
i
i
i
i
xxh
y
xxh
y
h
xx
m
h
xx
mxs
],[
1:1
)0()0(
),,()(2
1
33
0
2
3
时当于是、
ii
ii
n
xxx
ni
xsxs
xxCxs
1
1
1113
1
13
1
1
2
1
2
1
3
6
-)/-(
2
1
)0(
6
-)/-(
2
1
)0(
6
-)/-(
2
)(
2
)(
)(
i
ii
iiiiii
i
ii
iiiiii
i
ii
iii
i
i
i
i
i
i
h
mm
hyyhmxs
h
mm
hyyhmxs
h
mm
hyy
h
xx
m
h
xx
mxs
( 2 ) 1-n:1i
2
1:1 ),0()0(
11
33
iiiiii
ii
dmmm
nixsxs
)(
6
1
1
11
1
1
1
i
ii
i
ii
ii
i
ii
ii
i
i
h
yy
h
yy
hh
d
μ λ
hh
h
其中，
( 3 )
)0(,)0(3
00
3003
nn
nn
ym
ym
yxsyxs
、
综合上述，构造三次样条插值公式的算法为：
第一步：式（ 2）、（ 3）构成如下方程组
0112
2
21
1
1
1
1
1
11
0110
0
01
1
12
1
1
00
1
2
2
1
1111
)/()(6
2:2),/()(6
)/()(6
1:1,1,
,
2
2
2;),,( ;),,(
( 4 )
mhh
h
yy
h
yy
d
nihh
h
yy
h
yy
d
mhh
h
yy
h
yy
d
ni
hh
h
μ
ymym
A
dddmmm
dAm
nnn
n
nn
n
nn
n
ii
i
ii
i
ii
i
ii
ii
i
i
nn
n
n
T
n
T
n
注：
其中
],,[
.,,,,4
1
110
则若
）得到弯矩量方程组（第二步：利用追赶法解
ii
nn
xxx
mmmm?
( 1 )
h
)
6
m
-(
h-
)
6
m
-(
6
)(
6
)(
)(
i
2
1i
1
i
1
2
i
3
1
3
1
3
ii
i
ii
i
i
i
i
i
i
i
xxh
y
xxh
y
h
xx
m
h
xx
mxs
.
)m i n (
)m a x (
,2/)(),m a x (
,
8
3
24
1
3 8 4
5
)()()(
2],,[)(
1
2
3
4
2
4
i
i
i
h
h
chh
chR
hR
hR
hR
xsxfxR
baCxf
其中及其导数的界的估计为的余项的三次样条插值则问题定理：设三次样条插值的误差分析实验：用三次样条插值逼近 f(x),h=(b-a)/n.比较 n=4,8,16,32的效果
]1,1[,251 1)( 2 xxxf
均匀分划下的样条函数与样条插值一、基本 B样条函数从单位阶跃函数谈起
00
02/1
01
)(
x
x
x
xu
0阶 B样条函数：
)2/1()2/1()(0 xuxux
1阶 B样条函数：
2/1
2/1 01 10
11
)()(
x
x x
xx
dttx
2阶 B样条函数：
2/1
2/1
2
2
12
2/30
2/32/18/9*2/3*2/1
2/14/3
)()(
x
x
x
xxx
xx
dttx
3阶 B样条函数：
2/1
2/1
23
23
23
20
213/42*6/1
13/2*2/1
)()(
x
x
x
xxxx
xxx
dttx
一般的有 k次 B样条函数：
).()(
,2,1,)()(
1
2/1
2/1 1
RCx
kdxxx
k
k
x
x kk
且
实际上常用 k=1,2,3。而 k=3最用实用价值。
二、均匀分划下的样条函数逼近已知
n10 x x?x
n10 y y?y
njjhxx j,0,0
niyx
k
h
xx
cxs
iik
n
j
j
kjk
:0,)(s
B),()(
),3,k ( 1,2
0
使样条函数次确定对
B样条函数插值问题,
问题的解,
k=1时,一次 B样条插值公式,
)()(
0
11?
n
j
j
j h
xx
yxs
它就是分段线性插值。
k=2时,二次 B样条插值公式,
nn y
y
y
c
c
c
1
0
1
0
8
61
1
61
16
解
)()(
0
22?
n
j
j
j h
xx
cxs
得到误差分析结果与前面一般二次样条插值一样！但在边界点的光滑性不能保证。
k=3时,三次 B样条插值公式,
nn y
y
y
c
c
c
1
0
1
0
6
41
1
41
14
解得到
)()(
0
33?
n
j
j
j h
xx
cxs
误差分析结果与前面一般三次样条插值一样！
但在边界点的光滑性不能保证。
注：
为使 B样条插值在边界点也光滑，可以在两边延拓。以三次 B样条插值为例。即
nnii
n
j
j
j
yxyxniyx
h
xx
cxs
)(s)(s,:0,)(s
B3),()(
30033
1
1
33
，且使样条函数次确定解：
)1()(2)1(
)2/1()2/1()(
);2/1()2/1()(
211
223
223
xxx
xxx
xxx
( 1)
)2(
1
)2(
1
)(
1
)(
)(
1
)(
);(
1
)(
11
2
0101
2
3
1
1
2
3
3
1
1
2
03
3
1
1
2
3
nnnn
n
j
jn
n
j
j
j
n
j
j
yccc
h
yccc
h
jnc
h
xs
jc
h
xs
h
xx
c
h
xs
.2;23
6
6
6
41
1
41
14
2;6/,6/1
21
)2(
:0
64
:0,)(
2
110
2
101
1
2
01
1
2
1
2
0
2
00
11
3
nnnn
nnn
nnn
iiii
ii
yhcccyhccc
cy
y
cy
c
c
c
yhycyhyc
ni
yccc
niyxs
、
、
、
）得到）和（联立式（
实验,用一、二、三阶 B样条插值逼近 f(x),h=(b-a)/n.比较
n=4,8,16,32的效果
]1,1[,251 1)( 2 xxxf
n
为
4
的情况
n=8的情况
n=16的情况
n=32的情况阶梯函数与折线函数的磨光
).(,
)()(s
0,,,
0
0
00
00
xf
h
xx
yx
Bihxxyx
n
j
j
j
i
n
iii
简记即为阶梯（逼近）函数样条插值公式阶则对给定的型值点
样条插值函数。显然他就是一阶得到一次磨光函数进行一次磨光对阶梯函数
b
.)(
)(
1
)(
1
)(
)(
1
)(
0
1
0
2/
2/
0
2/
2/
01
2/
2/
0
0
n
j
j
j
n
j
hx
hx
j
j
hx
hx
hx
hx
h
xx
y
dt
h
xt
y
h
dttf
h
xf
dttf
h
xf
n
j
j
kjk h
xx
yxf
kxf
0
0
).()(
)( 次磨光（逼近）函数为的以此类推，
的逼近性质。、下面着重分析 )()( 32 xfxf
定理 1：
2
2/112/12
2
2/112/12
2
)4(
11
2
2
22
112
2/12
24
)(
)()(
1
)(
8
)(
)()(
2
1
)(
)3
8
)(
)()2(
1
)(
)2
)(
8
)(y
)2(
8
1
)(
1
)(
h
y
xyyy
h
xf
h
y
xyyyxf
h
y
xyyyy
h
xf
hohyyyyxf
xxxf
iiii
iiii
iiiii
i
iiiii
ii
、相切性质
、凹凸性质
）、逼近性质处有：及相邻节点的中点在节点
)33(
8
1
)(
1
)(
)(
48
1
)(
2
1
)(
)3
8
)(
)()2(
1
)(
)2
)(
6
)(y
)2(
6
1
)(
1
)(2
11212/13
21112/13
2
)4(
11
2
3
22
113
2/13
iiiiiii
iiiiiii
iiiii
i
iiiii
ii
yyyy
h
yy
h
xf
yyyyyyxf
h
y
xyyyy
h
xf
hohyyyyxf
xxxf
、相切性质
、凹凸性质
）、逼近性质处有：及相邻节点的中点在节点：定理
磨光法在外形设计问题重的应用例 某外形设计问题所给的型值 如下表。要求寻找这样的曲线，满足：
1）通过两个端点，
2）在 x(i)处逼近原型值的绝对误差不超过 1厘米，
3）在 x(i)处的二阶导数符号与原型值点的二阶差符号相同。
1112 2 ), iiiiii yyyyyx 的二阶差注：型值（
X(i)
cm
30 80 130 180 230 280 330 380 430 480
Y(i)
cm
80 110 132 148.75 163 175 185.5 195 204 212.75
02
02
1
,,,
)()(
).,(),,(
3)(
3)1:1()(2
.9,50,30
)
h
xx
(Ωy( x )f
3
11
101
11101
1
1
33
1111
03
3
0
n
0i
i
3i3
nnn
nnn
n
i
i
i
nn
n
i
yyy
yy y
yyhxxhxx
h
xx
yxf
yxyx
xxxf
nixxf
nhx
）得到：由条件待定。
此时三次磨光函数为两侧延拓，即增加右），必须将型值点向左处满足条件、在边界点要使
）。满足条件在内节点保证了显然，定理磨光函数）设计曲线可采用三次解：由条件
）。不满足条件直接计算表明由
）：现在我们来验证条件即：
2
)2(
6
1
)(
2
.2
,2
113
11
101
iiiii
nnn
yyyyxf
yyy
yyy
下面我们采用盈亏修正法来提高精度。其算法如下：
)(?)(?
.2?,2?,:0,?
)0
3
1
1
3
11101
h
xx
yxf
yyyyyyniyy
i
n
i
i
nnnii
开始（初始化）
）。转步形成新的三次磨光函数修正型值下一步。，停止计算；否则执行
）计算误差
1)4
)(?)(
)3
.2?
,2?
,:0,
)2
:0
)4?(
6
1
)(
1
3
1
1
3
11
101
113
h
xx
yxf
yyy
yyy
niryy
rif
ni
yyyyxfyr
i
n
i
i
nnn
iii
iiiiiii
实例计算结果
§ 3.4 数据拟合的最小二乘法问题的引入先讨论两个变量 x,y的情况。通过观测变量 x,y积累了一组资料（ x(i),y(i)),i=1:n，一般来说 n比较大。我们的任务是从这批实验数据出发，寻求一近似函数 φ(x)曲逼近 y。
由于观测数据都带有观测误差，数目又较大，对于这类问题运用插值函数取描述 y往往是不适当的。现在，用开始引入的例 1来说明建立近似函数的一种方法，即最小二乘法。
研究 y与 t的关系的最小二乘法，通常步骤如下：
1）先描 t,y的草图
2）凭视觉猜想一数学模型
t
t
t
tt
etb
tctccta
t
1
)(
1
ln))(ln ()(
)(:)(
)(:)(
/
2
210
，即本例可用两种试验模型
3）计算总误差
n
i
ii tycccI
a
1
2
210 ))((),,(
)(
来说是就试验模型
4）希望总误差最小
m in),,( 210?cccI
.2,1,0,0
m in),,( 210
j
c
I
cccI
j
m i n)),,,((),,,(
,),,,(,
:1),,,,,(
m,,,
1
2
2121
1
21
21
21
m
j
jnjjjn
n
i
iin
jjnjj
n
xxxycccI
xcxxx
mjyxxx
xxxy
则有若选择模型次观测得到观测值：有关系，并通过与多个变量如果数据拟合的最小二乘法
.:1,)()()(
:1,0 m i n)(
,,,
)(
.,,,m i n,)(
,)()(,:1,,
11 1
1
2
21
1
2
21
1
2
1
mkxycxx
mk
c
I
xyI
ccc
xy
cccxyI
xcxniyx
n
i
ikii
m
j
j
n
i
ikiji
k
n
i
iii
m
n
i
iii
m
n
i
iii
m
j
jjii
采用：确定参数越小，模型越好。
为均方误差。均方误差当参数确定后，称确定参数按最小二乘优化准则确定试验模型：根据观测数据
mmmmmm
m
m
n
i
ikiik
T
kk
n
i
ikijij
T
kjkkj
T
n
T
njjjj
b
b
b
c
c
c
aaa
aaa
aaa
mjk
xyyyb
xxa
yy
xxx
2
1
2
1
21
22221
11211
1
1
21
21
:1,
,)(,
,)()(),(
,,,,yy
,)(,),(),(
则令
对称矩阵方程例 1的数据拟合结果其中解：转化为采用平方误差：
从而拟合函数为解出为：正规方程采用
,
)()2(
9486.3
.04832 0.01436.11490.4)(
),04832 0.0,1436.1,1490.4(
01.8530
59.757
49.88
14043 410396826
1039682676
8267610
)()1(
/
2
2
2
1
0
2
210
bxaz
et
ttt
c
c
c
c
bAc
tctcct
t
.,ln,ln,1 baztx
);3 6 1 8.2,3 5 4 2.2,3 4 3 7.2,3 3 5 1.2,2 8 8 5.2,2 5 4 4.2,1 7 3 6.2,0 8 0 7.2,8 5 7 9.1,3 8 6 3.1(;0 6 2 5.0,0 7 1 4.0,0 8 3 3.0,1 0 0.0,1 2 5 0.0,1 6 6 7.0,2 5 0 0.0,3 3 3 3.0,5 0 0.0,0 0 0.1
,
z
x
zx 的数据对应：根据已知数据得到
1109.0
3411.11)(
.0579.1,3411.11,0579.1,4284.2
9586.4
4362.21
4930.16923.2
6923.210
2
/0 5 7 9.1
平方误差所求拟合函数为从而解出为正规方程
t
a
et
eba
b
a
bAc
显然，比第一种模型要好得多！
）的分布。和（拟合模型下图显示出原始数据,2)1(
内容包括：
1,MATLAB编程和绘图
2、数值分析的数学基础
3、数值算法在工程、科学中的应用第一章 MATLAB入门
1,MATLAB的命令窗口工作都在此处完成。
2、怎样进行计算运算对象：矩阵算术运算符,+ - * ^,*,/,^
倒除,\ a\b=b/a
变量与变量名,变量名和变量名类型不需声明 。
3、数据显示格式默认格式,5位（ format short)
format long 16位
format e 短的浮点格式
format long e 长的浮点格式
4、清除命令
clear:清除所有使用过的变量或某个（些）变量
clc,清除命令窗口
5、程序结构分支,if __else__end; if__elseif__end; if__break__end
循环,for__end; while__end
6、读写输入数据,z=input(‘type youe input:’) 键盘输入格式化输出,fprintf(‘e_format %12.5e \n’,vol)
7、数学函数
8、功能函数
sort(x) sum(x) max(x) min(x) mod(x,y) rand(n) eval(s)
9、编程（编写 M文件）
10、绘图第二章 数值代数内容：数值代数就是研究有关 矩阵计算 的问题。主要包括：
1,线性代数方程组的求解；
2,矩阵特征值问题要求,1、掌握用 MATLAB求解的方法
2、知道那些问题是困难的，那些问题是不可解的。
A=zeros(m,n) m行 n列的零矩阵
I=eye(n) n阶单位矩阵
A=ones(m,n) 元素均为 1
A’ A的转置
A(:,k) A(k,:) A(m1:m2,n1:n2)
inv(A) A的逆
size(A) A的大小
hilb(n) Hilbert矩阵
2.1 矩阵
bAx?
情况 1,m=n(正规方程），最常见；
情况 2,m<n(不定方程）；
情况 3,m>n(超定方程）；
本节只介绍情况 1。
倍计算时间大约是上面的效率最高
50 *)(
\
bAi n vx
bAx
MATLAB命令：
2.2 解线性代数方程组的 MATLAB命令线性代数方程组并不总是数值可解的。只有当矩阵 A
的行列式不为零时才行！矩阵 A的行列式即使不为零，但当很小或很大时，解的误差可能很大。
计算矩阵行列式的 MATLAB命令：
)d e t( AD?
2.4 病态问题有许多线性代数方程组理论上是可解的，但实际计算中由于受到舍入误差的影响而无法得到精确解。此类问题成为病态问题。
病态问题的计算过程中，小的舍入误差或系数矩阵的微小变化都可能使解产生很大误差。（例子 P97)
2.3 不可解问题病态矩阵的一个重要标志是条件数：
1)( AAAc o n d
MATLAB命令,)( Ac ond
当矩阵是病态时，其条件数一定很大，但它并不能直接说明解的误差。
线性方程组解的误差程度也取决于计算环境的精度。
条件数和行列式与计算环境是相互独立的。所以大条件数或小行列式未必意味无法直接精确求得线性方程组的解，
它只意味着有很大误差可能。而实际上如果采用更高精度的计算环境则很可能得到非常满意的解。
Hilbert矩阵是非常著名的病态矩阵（ hilb(n)),它经常用来检验算法的数值稳定性的好坏。
两种原因使我们想了解求解线性代数方程组的算法。
一是实际工作中要用其它计算机语言（ Fortran&C等）编写应用程序；二是 MATLAB处理大型稀疏矩阵方程组显得很笨拙或无能为力。
nnn
n
n
ba
baa
baaa
bABbAx
2222
111211
() ( 初等行变换）
2.5 线性代数方程组的求解方法（算法）
由线性代数的理论：
下面讨论实现过程：
第一步：消元。进行到第 k步时必有
nnnknnk
knkkkkk
kknkkkk
nkk
baaa
baaa
baaa
baaaa
1,
1,11,1,1
1,
111,1111
a(k,k)作为主元，第 k行依次乘 a(i,k)/a(k,k)加到第 i行,i=k+1:n。
总共 n-1步完成,k=1:n-1.
一,Gauss消去法
e nd
e nd
e nd
e nd
jkakiajiajia
:nkjf o r
kbkiaibib
kkakiakia
nkif o r
kkaif
nkf o r
);,(),(),(),(
1
);(),()()(
);,(/),(),(
:1
e n d ;b r e a k,,0),(
1:1
当 a(k,k)=0，则上述消去法无法进行；或当其绝对值相对太小可能会出现大的计算误差。选主元法可避免这种情况。下面介绍常用的 按列选主元的 Gauss法 。
e nd
e nd
e nd
jkakiajiajia
:nkjf or
kbkiaibib
kkakiakia
nkif or
tbkk btankkka
kkbkbnkkkankka
b( k ) tba( k,k,n) ;ta
kkkkk
knkaabskkq
nkf or
);,(),(),(),(
1
);(),()()(
);,(/),(),(
:1;)(;):,(
);()();:,():,(;1
) ) ) ;,:((m a x (],[
1:1
列主元 Gauss法运算量 (只考虑乘除运算）：
第 k步 =n-k次除法
+ n-k次乘法
+（ n-k)*(n-k)次乘法 ;
)(
3
1
)(2)(
3
23
1
1
1
1
2
no
nnn
knkn
n
k
n
k
总的乘除运算量 =
第二步：回代求解
1:1);,(/)),((
);,(/
1
niiiaxjiabx
nnabx
n
ij
jii
nn
e n d
iiadibib
e n d
jbjiadd
nijf o r
d
nif o r
nnanbnb
);,(/))(()(
);(),(
:1;0
1:1:1
);,(/)()(
%%%%%%%%
二,LU分解法
LU分解的目的是将矩阵 A转换为两个矩阵的乘积，即
。方程为上三角矩阵。此时，是单位下三角矩阵，
yUxbLybLU xbAx
LLUA
,
U,
好处是：对于线性方程组，如果需要多次求解不同的非齐次项，此时 LU分解的效率将大大超过高斯消去法。
LU分解的 MATLAB命令,[l,u]=lu(A) 和 [l,u,p]=lu(A)
前面讲到的不选主元的高斯消去法和列主元高斯消去法将能实现 LU
分解。
不选主元的高斯消去法用于下面两类矩阵肯定能成，即 严格对角占优矩阵或对称正定矩阵,其他矩阵就难说了！
列主元高斯消去法是解决一般中小型稠密矩阵方程组最有效的方法之一。下面讲解列主元高斯消去法实现 LU分解的算法。
.1:1
) ) ),,:((m a x (],[,1)(
.
1
),(
),(
1
),(
),1(
1
1
i)jP( i,
nk
knkaa b skkqkkkkkkkk
kka
kna
kka
kka
L
j
k
为其中行交换后的矩阵，行与第为单位矩阵第记
1,LU分解的代数理论
) ),(,(
,
333
22322
1131211
11222211
3
333
22322
1131211
1122
2
332
222
11211
11
kkkkPP
U
a
aa
aaa
aaaa
APLPLPLPL
aa
aa
aaa
aaaa
APLPL
aa
aa
aa
aaa
APL
k
nn
n
n
n
nnnn
nnn
n
n
n
nnn
n
n
n
其中，
于是：
,
.
~
,
,
~
,,
~
,
~;)(
)(
)(
))()((
121
121211
121121
121211
121
12121
121121
123211211
11222211
PPPP
PPLPPL
PPPLPPPL
PLPLLL
UAPPP
PPLPP
PPPLPPP
PPLPPPLPL
UAPLPLPLPL
n
nn
nkkkkknk
nnnnnn
nn
nn
nkkkkkn
nnnnnnnnn
nnnn
令：
.
~~~
,
.
~~~
1
1
1
2
1
1
121
为单位下三角矩阵。可以证明其中
L
LLLLLUPA
UPALLL
nn
nn
现在我们只要将列主元高斯消去法稍加改造即是 LU分解的算法。
列主元高斯消去法的矩阵表示：
.),,,
),:1(
U
))(,()(
.
'
,2,1
11222211
knkkkkk
k
nnnn
Laaa
knkAA
A
kkkkPPkkk
UAPLPLPLPL
唯一确定了（
部分，即的下三角的每一输出的矩阵；的上三角为输出的矩阵；唯一确定了
2,LU分解算法
e nd
e nd
e nd
jkakiajiajia
:nkjf or
kkakiakia
nkif or
tankkkka
nkkkkankka
a( k,k,n) ;ta
mkkkk
knkaabsmq
nkf or
as i z enm
aL G aus sakkf unc t i on
);,(),(),(),(
1
);,(/),(),(
:1;):),((
);:),(():,(;1)(
) ) ) ;,:((m a x (],[
1:1
);(],[
)(],[
运算量和高斯消去法一样
e nd
e nd
t;,k ) a( k k ( i );a( k k ( i ),k ) a( i,k )
a( i,k ) ; t
,n-k i f or
,n-f or k
e nd
t;k ( k ),,) p( k;p( k k ( k ),,)p( k,,)
kpt
nkf or
neyep
aL G aus sakk
as i z enm
aj i es anj i aof e npulf unc t i on
11
21
,) ;,(
1:1
);(
);(],[
);(],[
)(],,[
e nd
e nd
a( i,j ) ; u( i,j )
i,n j f or
:nf o r i
e nd
e n d
e nd
a( i,j ) ;i,j ) l (
e l s e;i,j ) l (
j i f i
:i j f or
:nf o r i
1
1
1
1
3、三对角矩阵方程的追赶法三对角且主对角严格占优矩阵方程是一类来源丰富的问题。比如，微分方程数值解或样条插值等问题中的正规方程组。解这种问题必须考虑其矩阵稀疏的特征，减少算法的计算量。
nn
n
ab
c
ab
ca
T
1
22
11
三对角矩阵形如：
fTx?问题：求解
T的 LU分解具有形式：
n
n
n u
c
u
cu
U
l
l
L
1
2
11
2
,
1
1
1
由 T=LU推得：
).:2(,/ ; 1111 niclauublau iiiiiii
1:1 ; /)(,/
:2 ;,
1
111
niuxcyxuyx
niylfyfy
yUx
fLy
fTx
iiiiinnn
iiii
e n d
1) ) / a ( i) ;f ( i*c ( i)-( f ( i)f ( i)
1:-1:1-nif or
f ( n ) / a ( n ) ;f ( n )
e n d
1) ;-f ( i*b( i)-f ( i)f ( i)
n:2if or
e n d
1) ;-c ( i*b( i)-a ( i)a ( i)
1) ;-b( i) / a ( ib( i)
n:2if or
le n g t h ( a ) ;n
f)c,b,a,z h u ig a n f a (ff u n c t io n
最终解为 f.
乘除运算量：
5n=O(n).
4、对称正定矩阵方程的 cholesky分解法对称正定矩阵方程的来源比较丰富，比如线性回归、拟合等问题。解这类问题必须考虑其矩阵对称的特征，减少算法的计算量。
为对称正定矩阵。其中问题：求解 A,bAx?
由于 A为对称正定矩阵,A必有 cholesky分解：
nnn
T
d
d
d
ll
l
L
L D LA
2
1
21
21
D,
1
1
1
,其中
j
j
k
jkikkijij
j
k
kjkjjj
ii
T
dlldal
dlad
nidalad
LD LA
/)(
n:1jin,:2j
:2,/,
1
1
1
1
2
111111
得到关系：
的两边，比较
e n d
e n d
j ) ;t ) / a ( j,-j)( a ( i,j)a ( i,
e n d
k) ;a ( j,*k)a ( i,*k)a ( k,t t
1-j:1kf or
0; t
n:1jif or
t;-j)a ( j,j)a ( j,
e n d
);,(*2)^,(t t
1-j:1kf or
0; t
n:2jf or
e n d
);a ( i,1) / d( 1a ( i,1)
n:2if or
s i z e ( a ) ;n][ m,
)c h ol e s ky ( aD][ L,
kkakja
f unc t i on
高斯消去法的一半。
，是很大时，相当于当运算量：乘除为上三角矩阵说明：
命令：
3
23
6
1
)29(
6
1
L,
)(
nn
nnn
LLA
Ac holLM A T L A B
T
三、线性代数方程组的迭代法线性代数方程组的迭代法并不适用于所有问题，但它对一些特定类型的问题非常有效。当问题是大型稀疏矩阵方程时，高斯消去法的效率会变得非常低，而且有时还会超出内存要求。对于这样的问题就需要使用迭代法。
大系统问题的求解最终归结为大型稀疏矩阵方程。比如，
电网络、场方程的数值计算、运筹问题等。
尽管迭代法的种类很多，这里只介绍其中的三种：
1,Jacobi迭代；
2,Gauss-Seidel迭代 ;
3、超松弛迭代法（ SOR)。
.x
,3;,2,1,0,2;1
)(*
*)(
)()1(
)0(
的精度为近似解为方程组的真解，其中收敛条件：、
迭代：、
、选定初始值：
算法：
k
k
kk
n
x
xx
kdBxx
Rx
dBxxbAx
1,线性代数方程组的迭代法的一般理论迭代算法简单，但问题是：
1、能否保证算法收敛？
2、能否充分利用矩阵的稀疏性，使运算量和存贮量尽量少。
[定理一 ] 迭代法收敛的充分必要条件是 0lim
1
k
k B[定理二 ] 迭代法收敛的充分必要条件是
1)(?B?
[定理三 ] 迭代法收敛的 充分条件 是
)1()(*( k )
)0(*( k )
1
x
*
1
x,1
kk
k
xx
G
G
x
xx
G
G
xB
和且谱范数行和范数列和范数
)(B
m a x; m a x
2
1
1
1
11
T
n
j
ij
ni
n
i
ij
nj
BB
bBbB
收敛性定理对于任意的初值 nRx?)0(
迭代矩阵 B的构造充分利用矩阵的稀疏性，使运算量和存贮量尽量少的办法就是要求迭代矩阵 B与原矩阵 A有相同的稀疏结构。具体就是：
0
a
0
aa0
U,
0aa
0a
0
L
);a,,a,d i a g ( aDU,LDA
n1,-n
1n12
1-nnn1
21
nn2211
)( bxULDbAx
n,1i,)x-(1) / axaa-(b
b)x )D 1-(-U( L ) x-(D
0,)
1
1(
1
SO R 3
n:1i,) / axaa-(b
) Se i d e l-G a u s s2
n,1i,) / axa-(bx
)( 1
( k )
ii
( k )
j
n
1ij
ij
)1(
1-i
1j
iji
)1(
( k )1)(k
ii
( k )
j
n
1ij
ij
)1(
1-i
1j
iji
)1(
)()1(
ii
( k )
j
n
ij1,j
iji
1)(k
i
)()1(
k
j
k
i
k
j
k
i
kk
kk
xx
DDD
xx
bUxxLD
bxULDxJ a c o b i
松弛因子。迭代：、
（迭代：、
迭代：、
常用迭代法及其收敛性
nRx )0(
定理 4,SOR收敛的必要条件是 0<ω<2.
定理 5：如果 A是严格对角占优矩阵或不可分弱对角占优矩阵，则
1,Jacobi收敛；
2,gauss-seidel收敛；
3、当 0<ω<=1 时,SOR必收敛。
定理 6：如果 A是对称正定矩阵，则
1、当 2D-A正定时,Jacobi收敛；
2,gauss-seidel收敛；
3、当 0<ω<2时,SOR必收敛。
其他矩阵不能保证收敛四、线性代数方程组共轭梯度法
)(m i n
0)()(
2
1
f ( x )
*
***
*
xfbAx
bAxxfxf
xA
xbAxx
TT
，因此由于
。唯一的极小点是对称正定矩阵时，有当二次函数：
即，解对称正定矩阵方程组的问题等价于求二次函数极小点的问题。
如何寻找 f(x)的极小点呢？基本思想就是 下降法 。
1、下降法的基本原理是：的极小值。具体的算法到以“最快的速度”下降使逐步产生一串点出发从某一点
)(
)()()f ( x
,,,,,,
)()1(( 0 )
)()1()0()0(
xf
xfxf
xxxx
k
k
为止。
或
、直到即的极小点上求、在射线；下降的方向处确定一个使、在
。现以求得
x
)f ( x 3
);(m i n)f ( x
,)(2
)(1
)(1)(k
( k )
)()(
0
1)(k
)1()()(
)()(
)(
k
kk
t
kkk
kk
k
x
tpxf
xxftpxx
pxfx
x
2、最速下降法我们知道，函数的负梯度方向是函数值下降最快的方向。
称为最速下降法作为下降方向的算法取 ( 1 ) )( )()()(?kkk Axbxfp
b)( A x)(ptAp)(p
pb)A ( p)tp(x
2
1
)tpA ( x)(p
2
1
)tpf ( x
dt
d
( k )T( k )( k )T( k )
( k )T( k )T( k )( k )( k )( k )T( k )
( k )( k )
( 3 ),x
( 2 )
)(
)(
,0)(
)()(1)(k
)()(
)()(
)()(
k
k
k
kTk
kTk
k
kk
ptx
App
pp
ttpxf
dt
d
得到时即当
.;
)(
)(;,,;
.0,1,2,k,
)()()1(
)()(
)()(
)(
)()(
)0(
k
k
kk
kTk
kTk
k
k
kk
n
rtxx
Arr
rr
t
e n db r e a krif
Axbr
Rx
最速下降算法：
function x=zuisu(a,b,x,eyta)
n=length(b);
r=b-a*x;
while max(abs(r))>eyta
t=r'*r/(r'*a*r);
x=x+t*r;
r=b-a*x;
end
这一算法有简单易行，可充分利用 A的稀疏性等特点，
但当 A的最大特征值远远大于最小特征值时收敛速度变得非常慢，以至于完全不适用。最速下降法并非最速！下面的共轭梯度法可有效解决这一问题。
3、共轭梯度法能否选出比负梯度方向更好的下降方向吗？能！
.0r4;1,,2,1,0,0)( 3
,,,,p2
,,,,1;
)(
)(
,;
)(
)(
,;,2,1,0,;,,
1)-(n
)()(
)()1(( 0 )
)()1()0(
)()(
)1()(
1
)(
1
)1()1(
)1()1(
)()(
)()(
)()()1(
)()(
)0()0()0()0()0(
、
、
共轭的；是两两关于、
是两两正交的；、可以证明：
其中；其中以求得，下降方向
kirp
App
rrr
App
rp
spsrp
Axbr
App
pp
tptxx
kpx
rpAxbrRx
kTi
k
k
kTk
kTk
k
k
k
kk
kk
kTk
kTk
k
k
k
kk
kk
n
).1(
10
*)(
)1(
nkxx
nr
k
n 步就能得到精确解计算说明上述计算过程最多这一算法称为共轭梯度法，简称 CG(conjugate gradient)
实际计算中由于误差的影响,)(kr 之间的正交性很快就会消失，以至于有限步内不能得到精确解，所以 CG实际是一种迭代算法。
CG是保稀疏且便于并行处理的算法，它是求解对称正定矩阵方程最优秀的算法之一。而且目前它已推广到一般矩阵方程，称为预优共轭梯度算法。
实用共轭梯度算法如下：
function x=cg(a,b,x,eyta)
n=length(b);
r=b-a*x;
p=r;
q0=r'*r;
while q0>eyta
w=a*p;
t=q0/(p'*w);
x=x+t*p;
r=r-t*w;
q=r'*r;
s=q/q0;
p=r+s*p;
q0=q;
end
运算量：每一步迭代的乘除运算量为
nn 22?
§ 2.6 矩阵特征值问题
).,(0)(,xxIARA nn 的解求问题,
一、幂法
MATLAB命令,D=eig(A)和 [X,D]=eig(A)
.,,,
,0,
21
21
n
n
xxx
nA
线性无关的特征向量个并且对应的特征值设
)1(
111
1
11111
)(
1111
)(
)1()()0()0(
)()(;0),(
,2,1,;0,
kkkk
kk
kkn
vxxv
xvk
kAvvvRv
显然常数充分大后，则当；?
.
),
(k
12
)(
)1()1(
)1(
)(
1
收敛因子值和特征向量。的按模最大的近似特征分别为分量的按模最大的是充分大时，结论：当
Av
vv
v
v
k
kk
lk
l
k
l
实用算法
：
时终止。当；
)1()(
)()(
)(
)(
)1()(
)0(;/
*))((
));(m a x (],[;
,2,1
];1,,1,1[
kk
k
kk
k
k
k
k
k
kk
vv
mvv
mlvs i gnm
vabslm
Avv
k
v
function [lmta,x]=mifa(a,epsl)
[m,n]=size(a);
x0=ones(n,1);
while 1
x=a*x0;
[lmta,m]=max(abs(x));
lmta=sign(x(m))*lmta;
x=x/lmta;
if abs(x-x0)<epsl,break,end;
x0=x;
end
一、幂法二、反幂法及其原点位移反幂法用来求 A的按模最小的特征值。
.,,,
,0,
21
121
n
nn
xxx
nA
线性无关的特征向量个并且对应的特征值设
n
k
k
nn
k
k
nn
m
xxv
A
xxAxAx
1
l i m
)m a x (/l i m
.
,
1
)(
1
1
1
幂法即可！收敛因子执行上述所以只要对因时终止。当；
)1()(
)()(
)(
)(
)1(1)(
)0(;/
*))((
));(m a x (],[;
,2,1
];1,,1,1[
kk
k
kk
k
k
k
k
k
kk
vv
mvv
mlvs i gnm
vabslm
vAv
k
v
实用算法
：
function [lmta,x]=fanmifa(a,epsl)
[m,n]=size(a);
x0=ones(n,1);
while 1
x=a\x0;
[lmta,m]=max(abs(x));
lmta=sign(x(m))*lmta;
x=x/lmta;
if max(abs(x-x0))<epsl,break,end;
x0=x;
end
lmta=1/lmta;
1、反幂法
2、带原点位移的反幂法反幂法与“原点位移”相配合，求指定点附近的某个特征值和特征向量，并可用于加速幂法的收敛性。
.
1
)(
sA,,0
)()(
,
s
m
sIA
liss
xsxsIAxAx
AsIAs
k
l
lli
用反幂法求得特征值，可对最近的的离点是即可知，若的原点位移。由为称对于给定的常数另一方面,s离 A的某个特征值 λ越近，收敛越快。因此不论用幂法求 A的按模最大特征值，还是利用反幂法求 A的按模最小特征值，为了加快收敛，均可以用迭代 m步后的近似值 lmta作为最初始的位移值,实行 动态位移迭代 。
动态位移幂法
function [lmta,x]=dongtamifa(a,epsl0,epsl)
[lmta,x]=mifa(a,epsl0);x0=x;lmta0=lmta;
n=length(x0);
while 1
x=(a-lmta0*eye(n))\x0;
[lmta,m]=max(abs(x));
lmta=sign(x(m))*lmta;
x=x/lmta;
lmta=1/lmta+lmta0;
if max(abs(x-x0))<epsl,break,end;
x0=x;
lmta0=lmta
end
动态位移反幂法
function [lmta,x]=dongtaifanmifa(a,epsl0,epsl)
[lmta,x]=fanmifa(a,epsl0);x0=x;lmta0=lmta;
n=length(x0);
while 1
x=(a-lmta0*eye(n))\x0;
[lmta,m]=max(abs(x));
lmta=sign(x(m))*lmta;
x=x/lmta;
lmta=1/lmta+lmta0;
if max(abs(x-x0))<epsl,break,end;
x0=x;
lmta0=lmta;
end
三、幂法综述
1、幂法和反幂法只能用于求解 可对角化矩阵的实数特征值和特征向量，不能求解复数特征值 ；
2、幂法和反幂法均为线性收敛，收敛速度由收敛因子决定，效率不高。
3,动态位移可以大大减小收敛因子加速收敛 ；
4,不适于求解全部特征值。
5、对称矩阵自然适用于幂法，此时采用 2-范数标准化向量的算法至少平方收敛。
求解一般矩阵全部特征值的办法是下面的相似变换法四、矩阵全部特征值的 QR迭代算法求解一般矩阵全部特征值的办法是相似变换法。其理论根据是：
。均为一阶或二阶实方阵其中使正交矩阵
ii
ss
s
s
nn
A
A
AA
AAA
TAPP
PRA
,
,,
222
11211
1
.l i m
,,)Q (A
).(],[
,
.,2,1,A,A,
)(
m
)()1)()()(1-( m )1)(m)(
)()(1)(m)()(( m ))0(
TA
AAARQQ
AqrRQM A TL A B
RQQRA
mQRRQAA
m
mmmmmm
mmmm
且相似于即命令：交三角化分解。
解为正为上三角矩阵，称此分为正交矩阵，其中，
实现过程：
（
1、矩阵的两种正交变换平面旋转变换变换。为平面旋转矩阵，或则矩阵记
G i v e n s
k i
1
cs-
sc
1
)k,J ( i,
,1,s i n,c o s
nn
R
nkisc
kijxy
y
xx
xx
x
s
xx
x
c
kijxy
cxsxy
sxcx
JxyRx
IJJkiJJ
jj
k
ki
ki
k
ki
i
jj
kik
ki
n
T
,,
0
y
,,0 x
,,
y
,
.),,(
22
i
2222
k
i
则时，取一个非零元素。即它的意义是消去向量的的分量表示正交矩阵，即?
镜面反射矩阵
.)0,,0,1(.
2
),(
2
1
),,,(
s i g n ( x,/) ),(m a x (
,,,),,,(2
,1
2,1,
11
1
22
1
2
2
211
2
1
221
1
2
T
T
T
T
n
n
T
n
T
Tn
eeHxH
uuI
u
u
u
u
IH
txu
xxxexu
xtxxxa bst
xxxxxx
HHH
H
rH o us e ho l d e
IHR
使确定的镜面反射矩阵
）
不全为零，则由）（
对称正交；）（
具有下列性质：矩阵变换矩阵。或称为镜面反射矩阵，矩阵镜面反射矩阵的意义是“成批”消去向量的非零元素。
function [a,b,x]=householder(x)
% x gived,seek a househouder convert H=I-bxx',make Hx=-a.
n=length(x);
t=max(abs(x));
if t==0
a=0;
b=0;
else
x=x/t;
a=sqrt(x'*x);
if x(1)<0,a=-a;end
x(1)=x(1)+a;
b=1/(a*x(1));
a=t*a;
end
2、方阵的正交三角化分解方阵的正交三角化分解,即为上三角矩阵。为正交矩阵，其中 RQQRA,?
;*
**;~
)1(;
*;
~
) ) ;1(,,(],,[;
*;) ) ;1(,,(],,[
2
2
1
12
2
2
2
2
12212
1
1
111
A
a
a
AHH
H
e y e
H
A
a
AHb x xIHArh ou s e ho l d exba
A
a
AHb x xIHArh ou s e ho l d exba
T
T
分解方法：
function [arfa,bata,a]=qrfenjie(a)
[m,n]=size(a);
for k=1:n-1
[arfa(k),bata(k),a(k:n,k)]=householder(a(k:n,k));
for j=k+1:n
a(k:n,j)=a(k:n,j)-bata(k)*a(k:n,k)'*a(k:n,j)*a(k:n,k);
end
end
function [q,r]=qrgouzao(a)
%实现 a=qr,q是正交矩阵,r是上三角矩阵。
[m,n]=size(a);
[arfa,bata,a]=qrfenjie(a);
q=eye(n)-bata(1)*a(1:n,1)*a(1:n,1)';
for k=2:n-1
q(1:k-1,k:n)=q(1:k-1,k:n)-bata(k)*q(1:k-1,k:n)*a(k:n,k)*a(k:n,k)';
q(k:n,k:n)=q(k:n,k:n)-bata(k)*q(k:n,k:n)*a(k:n,k)*a(k:n,k)';
end
for i=1:n-1
for j=i+1:n
r(i,j)=a(i,j);
end
end
for i=1:n-1
r(i,i)=-arfa(i);
end
r(n,n)=a(n,n);
3、化矩阵为 Hessenberg型矩阵。不可约时矩阵。所有称为（上）
形如：
H e s s e nb e r g
hH e s s e nb e r g
hh
hhh
hhhh
hhhh
H
ii
nnnn
nn
nn
nn
,0
,
,1
1,
31,332
21,22221
11,11211
矩阵。为使得正交矩阵定理：
H e s s e n b e r g
AQQBQRA Tnn,,
function [q,h]=hessenberghua(a);
% h=q'aq,h为 Hessenberg矩阵,q为正交矩阵。
[m,n]=size(a);
q=eye(n);
for k=1:n-2
[arfa(k),bata(k),a(k+1:n,k)]=householder(a(k+1:n,k));
p=eye(n-k)-bata(k)*a(k+1:n,k)*a(k+1:n,k)';
a(1:k,k+1:n)=a(1:k,k+1:n)*p;
a(k+1:n,k+1:n)=p*a(k+1:n,k+1:n)*p;
q(1:n,k+1:n)=q(1:n,k+1:n)*p;
end
h=a;
for i=1:n-2
h(i+1,i)=-arfa(i);
end
for j=1:n-2
for i=j+2:n
h(i,j)=0;
end
end
注：当 A为对称矩阵时，其
Hessenberg 形为对称三对角矩阵。
4,Hessenberg矩阵在 QR分解下的不变性变换下的不变性。矩阵这种性质称为变换，的到矩阵。称此为由均为和下，
正交相似变换，、
分解，、
矩阵，则在为设
QRH e s s e n be r g
QRHHH e s s e n be r gHQ
RQHQQH
QRH
H e s s e n be r gH
T
~~
~
2
QR 1
显然，实现 Hessenberg矩阵的 QR分解用 Givens变换合适。即执行 n-1次 Givens变换,
运算量。一过程需完成这然后计算为上三角形使
2
121
121
4
.
~
,
,1:1),,1,(
n
JJRJH
RHJJJniiiJJ
T
n
TT
nnii
function [q,h]=givensqr(h)
%h为不可约 hessenberg矩阵，用 givens变换进行 QR分解。
[m,n]=size(h);
q=eye(n);
for k=1:n-1
if h(k+1,k)~=0
c=h(k,k)/sqrt(h(k,k)^2+h(k+1,k)^2);
s=h(k+1,k)/sqrt(h(k,k)^2+h(k+1,k)^2);
d=[c,s;-s,c];
h(k:k+1,k:n)=d*h(k:k+1,k:n);
if k==1
q=[d' zeros(2,n-2);zeros(2,n-2)' eye(n-2)];
else
q=q*[eye(k-1) zeros(k-1,2) zeros(k-1,n-k-1)...;zeros(k-1,2)' d' zeros(2,n-k-1)...;zeros(k-1,n-k-1)' zeros(2,n-k-1)' eye(n-k-1)];
end
end
end
5、基本 QR算法
QR( 2 ) A
)(A 1)
,2,1,0,2
)(],[P
1
kk1k
121k
00
0
T
n
TT
kkk
nn
JJJQRQ
kQR
Ahu aH e s s e n be r gA
AH e s s e n be r gA
RA
注：（
迭代、
形为、约化算法：
一般 QR 迭代算法的收敛性比较复杂，这里不再介绍。仅指出，在一定条件下由基本 QR算法生成的系列
kA 收敛为准上三角矩阵，对角块按特征值的模从大到小排列。
function [q,h]=jibenQR(a,epsl)
[m,n]=size(a);
[q0,h]=hessenberghua(a);
t=zeros(1,n-1);
while 1
for i=1:n-1
if abs(h(i+1,i))<=(abs(h(i,i))+abs(h(i+1,i+1)))*epsl
h(i+1,i)=0;
t(i)=i+1;
end
end
k=1;
while k<=n-2
if t(k)==t(k+1)
yes=1;
k=n-1;
else
k=k+1;
yes=0;
end
end
if yes==0,break,end;
[q,h]=givensqr(h);
h=h*q;
q=q0*q;
q0=q;
end
6,QR 算法的改善基本 QR算法计算量和存储量都很大，且若收敛是线性的。
因此，需要从这两方面加以改善。
（ 1）划分和收缩
（ 2）原点位移第三章 数值逼近内容：
数值逼近也称数据模拟，是为离散数据（不论何种途径得到）建立简单连续模型的方法。包括：
1、插值方法
2、数据拟合方法
3、最佳逼近要求：
1、了解 MATLAB功能函数
2、掌握数据模拟的各种方法及应用
§ 3.1 问题的提出例一,在某化学反应里，侧的生成物的质量浓度 y与时间
t(min)的关系入表。为了研究该化学反应的性质，如反映速率等,欲求 y与 t之间的关系 y=f(t)。
1 2 3 4 6 8 10 12 14 16
4.00 6.41 8.01 8.79 9.53 9.86 10.33 10.42 0.53 10.61
例二,逢山开道（山区修建公路的路线选择问题 1994）
在某山区修建公路,已知该山区的地形高度见表一 （ 单位
m） 。 雨季在山谷形成一溪流,雨量最大时溪流水面宽度 W
与 x（ 溪流最深处的 ） 坐标的关系可以表示为
4 00 02 40 0,522 40 0
4/3
xxxw
要求：从山脚开始经居民点至矿区修一条公路 （ 一般公路,
桥,隧道 ） 给出路线设计方案使工程成本最小 。
3200 1430 1450 1460 1550 1500 1600
2800 950 1190 1370 1500 1200 1100
2400 910 1090 1270 1500 1200 1100
2000 880 1060 1230 1390 1500 1500
1600 830 980 1180 1320 1450 1420
1200 740 880 1080 1130 1250 1280
800 *650 760 880 970 1020 1050
400 510 620 730 800 850 870
0 370 470 550 660 670 690
Y/x 0 400 800 1200 1600 2000
部分数据表：
工程种类 一般路段 桥梁 隧道工程成本
（ 元 /m）
300 1500 （ 长度 <=300m ）,3000 （ 长度
>300m）
α=上升高度 /
水平距离
<=0.125 =0 <=0.100
2000
工程要求：
例三,水道测量模型 （ 1989.美国 ）
问题 水平方向的坐标 x,y以 Td（ =0.914m） 为单位,水深方向 Z 以 Ft(=30.48cm)为单位,下表给出了水面直角坐标 (x,y)处的水深 Z,这是在低潮时测得的 。 如果船的吃水深度为 5Ft,试问在矩形城 (75,200)(-50,50)中行船应避免进入那些区域?
工作,1,要根据数据做出 ① 地貌图 （ 三维 ）,② 等高线 （ 等势图 ） 二维 。 2,由坡度限制,沿给定的网格点选择路线
x 129.0 140.0 108.5 88.0 185.5 195.5 105.5
y 7.5 141.5 28.0 147.0 22.5 137.5 85.5
z 4 8 6 8 6 8 8
157.5 107.5 77.0 81.0 162.0 162.0 117.5
-6.5 -81.0 3.0 56.5 -66.5 84.0 -38.5
9 9 8 8 9 4 9
x
y
z
工作：将海底曲面图绘制出来 。
总之：根据给定采样数据而定出其实际分布图,这就是数值模拟问题,采用的方法就是插值和拟合 。
§ 3.2 数值模拟概论在科学与工程等实际问题中，实际数学模型或者难于建立或者难于求解，但其数据模型（由实验或测量所得到的一批离散数据）容易得到。能否通过处理这些数据来建立实际模型呢？ 这里我们仅以一维问题来说明。
给定,y=f(x)数据
n10 x xx?
n10 y y y?
通过这批数据希望找到对象 y=f(x)的更多的信息（或全部信息）。通常只能找到 f(x)的近似表达式 φ(x),φ(x)是根据研究对象的特征确定的。
0
s i nc o s
:
n
nnnn xbxax
，则研究对象是周期振动的比如对象是指数增（减）并最终稳定到某个值,则
xeccx 10
对象是直线运动、匀加速运动，则
xaax 10
2210 xaxaax
一般地，我们都要通过对研究对象实际的了解，做以下工作：
1、确定它的基本特征函数
);(,),(,21 xxx m
2、对这些特征函数（基函数）进行一种恰当的构造，得到
f(x)近似的数学模型,
非线性的,.,,,,,,; 21 mcccx?
最常用也就是最简单的一种构造方式就是,
线性
1
xcx jj
m
j
.)(,,,21 的线性组合即 xm
3、按某种寻优策略，由已知数据确定未知参数
mccc,.,,,,,,21寻优策略的不同，产生了不同的数值逼近方法
§ 3.3 插值方法寻优策略 就是 过点,即
niyx jj,,2,1,
根据需要可附加补充原则,
np xxcx,0
nxx,0
1,p阶光滑度
2、边界条件,处光滑性,周期性等
（ 整体） 误差,
),(,0 nxxxxxfxR
插值方法中最常用的一类就是 代数插值 。
代数插值：即多项式的插值 m
m xaxaax10?
一、代数插值
（一）拉格朗日插值公式（ Lagrange)
n10 x xx?
n10 y y y?
给定,y=f(x)数据确定次数不超过 n的多项式 P(x),使
niyx ii,,2,1,0,( ）
).(
)!1(
)(
)()()(
)()(
L a g r a n g e;,,2,1,0,
)()(
)(
)(;)()(
)1(
0
0
x
n
f
xxfxR
xlyx
nj
xxx
x
xl
xxx
n
n
j
jj
jj
j
n
j
j
整体误差：
插值公式为：则令
拉格朗日插值公式的讨论
。其中 ))((m a x
8
)(;)(
],[
2
2
01
2
1
01
0
0
10
1
1
10
xfM
xx
M
xR
y
xx
xx
y
xx
xx
x
xxx
线性插值公式（ n=1的情况）
显然，两插值节点越近，越精确。
2、二次插值（抛物线插值,n=2的情况）
).)()((
!3
)(
)(;
))((
))((
))((
))((
))((
))((
)(
2102
2
1202
10
1
2101
20
0
2010
21
2
xxxxxx
f
xR
y
xxxx
xxxx
y
xxxx
xxxx
y
xxxx
xxxx
x
n越大精度越高！是这样吗？下面我们看一实例。
实例演示,分别用 n=1,2,4,8,16,32时的拉格朗日插值逼近
f(x)。
]1,1[,251 1)( 2 xxxf
function f=lagelangri(t,y,x)
n=length(t);
for j=1:n
l(j)=1;
for i=1:n
if i~=j
l(j)=l(j)*(x-t(i))/(t(j)-t(i));
end
end
end
f=y*l';
定义 lagrange函数
function f=shiyanhs(x)
f=1./(1+25*x.^2); 定义实验函数
function g=lagelangritu(a,b,n)
h=(b-a)/n;
t=a:h:b;
f0=shiyanhs(t);
x=a:0.02:b;
m=length(x);
for i=1:m
p(i)=lagelangri(t,f0,x(i));
end
f=shiyanhs(x);
plot(x,f,'r',x,p);
gtext('f(x)');
gtext('p(x,n)');
axis([-1,1,-1,1]);
画图,f(x) 和 p(x)
总之，我们看到利用拉格朗日公式逼近函数 f(x),n小不行，
n大也不行。这是为什么呢？
，看出从误差：
)())(()!1( )()( 10
)1(
n
n
xxxxxxnfxR
当然大。内取的插值节点少，小代表在区间、
而增加。
的增加附近的值随、在边界的均值附近比较小，而的值，在（、
跑得快！（
比的增加呈指数级增长，的值，常常随、
)(],[3
n
,,,)())(()2
)!1
)(1
0
1010
)1(
xRban
xx
xxxxxxxxxx
n
nxf
n
nn
n
差。计算量和计算稳定性变的增加，随着
（（（（（
（（（（（
插值公式来看：、从
)(
)()(;,,2,1,0
,
)))))
)))))
)(
4
0
1100
1100
xln
xlyx
nj
xxxxxxxxxx
xxxxxxxxxx
xl
Lag r an ge
j
n
j
jj
njjjjjjj
njj
j
综合以上分析我们可以得到如下结论：
在小区间上应用低次拉格朗日差值公式（ n<=6)有效。
问题是对大的区间怎么办？那就是用 分段低次插值。
（二）分段低次插值分段线性插值,
n10 x xx?
n10 y y y?
给定,y=f(x)数据确定 P(x),使
niyx ii,,2,1,0,(2 ）、
上是线性函数。在每一小区间,],[)(1 1 ii xxxP?
根据要求立即得到：
即为所求！
.
)(
.,,2,1],,[
1
1
1
1
1
1
1
1
iii
i
i
i
i
i
i
i
ii
i
i
ii
i
ii
xxh
y
h
xx
y
h
xx
y
xx
xx
y
xx
xx
xP
nixxx?
若用基函数的形式表示，则：
0
)(
11
10
10
1
0,令
xxx
xxx
xx
xx
xN
即为所求！);(yP ( x )
0
)(N
1.-n,1,2,k
]x,[x x 0
)(N
n
0k
k
10
1
1
1
n
1k1-k
1
1
1
1
1
1
k
xN
xxx
xxx
xx
xx
x
xxx
xx
xx
xxx
xx
xx
x
k
n
nn
nn
n
kk
kk
k
kk
kk
k
function f=fenduanxianxing(t,y,x)
% x在 [t(1),t(n)]
n=length(t);
i=1;
while 1
if (x>=t(i))&(x<=t(i+1)),break,end;
i=i+1;
end
f=y(i)*(x-t(i+1))/(t(i)-t(i+1))+y(i+1)*(x-t(i))/(t(i+1)-t(i));
定义分段线性插值函数
]1,1[,251 1)( 2 xxxf
实验：用分段线性插值逼近函数
function g=fenduanxianxingtu(a,b,n)
h=(b-a)/n;
t=a:h:b;
f0=shiyanhs(t);
x=a:0.02:b;
m=length(x);
for i=1:m
f(i)=fenduanxianxing(t,f0,x(i));
end
f1=shiyanhs(x);
plot(x,f1,'.-r',x,f);
gtext('f(x)');
gtext('p(x,n)');
axis([-1,1,-1,1]);
g=fenduanxianxingtu(-1,1,50)
计算：
分段线性插值的讨论：
误差分析
2
0
)2(
],[
)2(
],[
2
1
8
R ( x )
],,[
,)(m a x,)m a x (
)(m a x,
8
R ( x )
],[
0
1
h
M
xxx
xfMhh
xfMh
M
xxx
n
xxx
i
i
xxx
ii
i
ii
n
ii
有则令时，有更多的插值节点。、逼近精度提高，需要只是连续，光滑性差；、
缺点：
好。计算简单，数值稳定性、
时，、
优点：
2
)(1
)(2
);()(01
xP
xP
xfxPh
下面的实例演示可以说明
g=fenduanxianxingtu(-1,1,n)
分别用 n=4,8,16,32的分段线性差值逼近函数（ h=2/n)
]1,1[,251 1)( 2 xxxf
样条函数与样条插值,
n10 x xx?
n10 y y y?
给定,y=f(x)数据确定 s(x),使
niyxs
xxCxs
ii
n
k
,,2,1,0,(3
),()(2 01
）、
、
)2
],[)(1 1
kk
xxxs ii
的多项式。（
上是次数不超过在每一小区间、
若 s(x)存在，则称 s(x)为对应于给定数据的一个 k次样条插值。满足条件 1和 2的函数 S(x)称为 k次样条函数。
样条函数的这种数学定义来源于 绘图员沿着受压铁约束的样条（一根具有很好弹性的木条），画出各种光滑曲线（样条曲线）。这种样条曲线可以看成弹性细梁受集中载荷作用而生成的 挠度曲线。
S(x)的存在性讨论：
).()(
.,2,1,0
0 0
0 !/
)(
1
m
RCx
m
x
xmx
xx
m
m
m
m
m
次半截幂函数定义：
是待定常数。其中次样条函数可表示为：
jj
n
j
k
jj
k
j
j
jk
ba
kxxbjxaxs
k
,
!)(!)(
1
10
实际中常用 k=2,3的情况，即二次样条和三次样条个待定常数。共有
：二次样条函数可表示为
2
!2)(
!2
)(
1
1
2
2
2102
n
xxb
x
axaaxs
n
j
jj
个待定常数。共有
：三次样条函数可表示为
3
!3)(!)(
1
1
3
3
0
3
n
xxbjxaxs
n
j
jj
j
j
j
二次样条插值公式因为二次样条函数的确定需要 n+2个条件，而
niyxs ii,0 )(2
给定了 n+1个，所以还需补充一个条件（自然是边界条件）。通常有
).()(,)()( 2
))(()(1
202
'
20
'
2
''
2
01
01'
00
'
2
nn
nn
xsxsxsxs
yxs
xx
yy
yxs
、周期性条件或、
补充条件类型：
补充条件 1下的二次样条插值公式
( 2 ) /)(
y)x-(x
)(],,[
( 1 ) 1:0
,/)(,
)(,)(
0
'
0
0
01
0
'
01000
'
00
'
210
11i
1122
hy
h
yy
c
c
yxsxxx
ni
xxhhyyy
yxsyxs
iiiiiiii
iiii
得到由得到由
1:1),0()0(),()(2
1:0),)(()()(
],,[ ],[)(1
'
2
'
210
1
2
12
112
nixsxsxxCxs
nixxxxcxxxs
xxxxxxs
ii
iiiiii
iiii
、
内为二次多项式在、
采用如下方法来确定：
( 3) 1:1
/)(
1:1),0(s)0(s
1
11
1
1
111-i
'
2
'
2
ni
h
h
yy
h
yy
c
h
h
c
hchc
nixx
i
i
ii
i
ii
i
i
i
i
iiiii
ii
得到由联立式（ 1）、（ 2）、（ 3）确定了二次样条插值公式：
))(()()(
],[,],,[;1:0,
12
10
i
iiiiii
iin
ii
xxxxcxxxs
xxxixxx
nic
第三步：；即确定第二步：
、、第一步：计算算法个待定常数。的确定立求解：当然，我们也可采用联
2
!2)(
!2
)(
)0(
:0,)(
1
1
2
2
2102
'
00
'
2
2
n
xxb
x
axaaxs
yxs
niyxs
n
j
jj
ii
这样不如前面的方法简单！
实验：用二次样条插值逼近 f(x),h=(b-a)/n.比较 n=4,8,16,32的效果
]1,1[,251 1)( 2 xxxf
二次样条插值的误差分析
).(
12
5
3
4
),m a x (
,
)()()(
1],,[)(
)4(
2
3
2
4
abffchh
chR
chR
chR
xsxfxR
baCxf
i
其中及其导数的界的估计为的余项的二次样条插值则补充条件定理：设三次样条插值公式因为三次样条函数的确定需要 n+3个条件，而
niyxs ii,0 )(3
给定了 n+1个，所以还需补充二个条件（自然是边界条件）。通常有
).0()0(
)0()0(
)()( 3
)0()0(2
)0()0(,1
202
'
20
'
2
202
3003
''
3
'
00
'
3
n
n
n
nn
nn
xsxs
xsxs
xsxs
yxsyxs
yxsyxs
）（自然有：周期性条件问题
，：问题
，问题补充条件类型：
解由插值条件和补充条件构成的线性方程组，可以得到三次样条插值公式。但这样比较复杂！
下面只就问题 2来讨论，并采用另一种途径称之为三 弯矩法。
)(,)(;)(
],,[
.:0,)(],[)(
),()(],[)(1
1133
1
1
3
1
313
10
2
313
得并使其满足对此式做两次不定积分则内是线性函数。设在故
，内为三次多项式且在、
iiii
i
i
i
i
i
i
ii
iiii
ii
yxsyxs
h
xx
m
h
xx
mxs
xxx
nimxsxxxs
xxCxsxxxs
( 1 )
h
)
6
m
-(
h-
)
6
m
-(
6
)(
6
)(
)(
i
2
1i
1
i
1
2
i
3
1
3
1
3
ii
i
ii
i
i
i
i
i
i
i
xxh
y
xxh
y
h
xx
m
h
xx
mxs
],[
1:1
)0()0(
),,()(2
1
33
0
2
3
时当于是、
ii
ii
n
xxx
ni
xsxs
xxCxs
1
1
1113
1
13
1
1
2
1
2
1
3
6
-)/-(
2
1
)0(
6
-)/-(
2
1
)0(
6
-)/-(
2
)(
2
)(
)(
i
ii
iiiiii
i
ii
iiiiii
i
ii
iii
i
i
i
i
i
i
h
mm
hyyhmxs
h
mm
hyyhmxs
h
mm
hyy
h
xx
m
h
xx
mxs
( 2 ) 1-n:1i
2
1:1 ),0()0(
11
33
iiiiii
ii
dmmm
nixsxs
)(
6
1
1
11
1
1
1
i
ii
i
ii
ii
i
ii
ii
i
i
h
yy
h
yy
hh
d
μ λ
hh
h
其中，
( 3 )
)0(,)0(3
00
3003
nn
nn
ym
ym
yxsyxs
、
综合上述，构造三次样条插值公式的算法为：
第一步：式（ 2）、（ 3）构成如下方程组
0112
2
21
1
1
1
1
1
11
0110
0
01
1
12
1
1
00
1
2
2
1
1111
)/()(6
2:2),/()(6
)/()(6
1:1,1,
,
2
2
2;),,( ;),,(
( 4 )
mhh
h
yy
h
yy
d
nihh
h
yy
h
yy
d
mhh
h
yy
h
yy
d
ni
hh
h
μ
ymym
A
dddmmm
dAm
nnn
n
nn
n
nn
n
ii
i
ii
i
ii
i
ii
ii
i
i
nn
n
n
T
n
T
n
注：
其中
],,[
.,,,,4
1
110
则若
）得到弯矩量方程组（第二步：利用追赶法解
ii
nn
xxx
mmmm?
( 1 )
h
)
6
m
-(
h-
)
6
m
-(
6
)(
6
)(
)(
i
2
1i
1
i
1
2
i
3
1
3
1
3
ii
i
ii
i
i
i
i
i
i
i
xxh
y
xxh
y
h
xx
m
h
xx
mxs
.
)m i n (
)m a x (
,2/)(),m a x (
,
8
3
24
1
3 8 4
5
)()()(
2],,[)(
1
2
3
4
2
4
i
i
i
h
h
chh
chR
hR
hR
hR
xsxfxR
baCxf
其中及其导数的界的估计为的余项的三次样条插值则问题定理：设三次样条插值的误差分析实验：用三次样条插值逼近 f(x),h=(b-a)/n.比较 n=4,8,16,32的效果
]1,1[,251 1)( 2 xxxf
均匀分划下的样条函数与样条插值一、基本 B样条函数从单位阶跃函数谈起
00
02/1
01
)(
x
x
x
xu
0阶 B样条函数：
)2/1()2/1()(0 xuxux
1阶 B样条函数：
2/1
2/1 01 10
11
)()(
x
x x
xx
dttx
2阶 B样条函数：
2/1
2/1
2
2
12
2/30
2/32/18/9*2/3*2/1
2/14/3
)()(
x
x
x
xxx
xx
dttx
3阶 B样条函数：
2/1
2/1
23
23
23
20
213/42*6/1
13/2*2/1
)()(
x
x
x
xxxx
xxx
dttx
一般的有 k次 B样条函数：
).()(
,2,1,)()(
1
2/1
2/1 1
RCx
kdxxx
k
k
x
x kk
且
实际上常用 k=1,2,3。而 k=3最用实用价值。
二、均匀分划下的样条函数逼近已知
n10 x x?x
n10 y y?y
njjhxx j,0,0
niyx
k
h
xx
cxs
iik
n
j
j
kjk
:0,)(s
B),()(
),3,k ( 1,2
0
使样条函数次确定对
B样条函数插值问题,
问题的解,
k=1时,一次 B样条插值公式,
)()(
0
11?
n
j
j
j h
xx
yxs
它就是分段线性插值。
k=2时,二次 B样条插值公式,
nn y
y
y
c
c
c
1
0
1
0
8
61
1
61
16
解
)()(
0
22?
n
j
j
j h
xx
cxs
得到误差分析结果与前面一般二次样条插值一样！但在边界点的光滑性不能保证。
k=3时,三次 B样条插值公式,
nn y
y
y
c
c
c
1
0
1
0
6
41
1
41
14
解得到
)()(
0
33?
n
j
j
j h
xx
cxs
误差分析结果与前面一般三次样条插值一样！
但在边界点的光滑性不能保证。
注：
为使 B样条插值在边界点也光滑，可以在两边延拓。以三次 B样条插值为例。即
nnii
n
j
j
j
yxyxniyx
h
xx
cxs
)(s)(s,:0,)(s
B3),()(
30033
1
1
33
，且使样条函数次确定解：
)1()(2)1(
)2/1()2/1()(
);2/1()2/1()(
211
223
223
xxx
xxx
xxx
( 1)
)2(
1
)2(
1
)(
1
)(
)(
1
)(
);(
1
)(
11
2
0101
2
3
1
1
2
3
3
1
1
2
03
3
1
1
2
3
nnnn
n
j
jn
n
j
j
j
n
j
j
yccc
h
yccc
h
jnc
h
xs
jc
h
xs
h
xx
c
h
xs
.2;23
6
6
6
41
1
41
14
2;6/,6/1
21
)2(
:0
64
:0,)(
2
110
2
101
1
2
01
1
2
1
2
0
2
00
11
3
nnnn
nnn
nnn
iiii
ii
yhcccyhccc
cy
y
cy
c
c
c
yhycyhyc
ni
yccc
niyxs
、
、
、
）得到）和（联立式（
实验,用一、二、三阶 B样条插值逼近 f(x),h=(b-a)/n.比较
n=4,8,16,32的效果
]1,1[,251 1)( 2 xxxf
n
为
4
的情况
n=8的情况
n=16的情况
n=32的情况阶梯函数与折线函数的磨光
).(,
)()(s
0,,,
0
0
00
00
xf
h
xx
yx
Bihxxyx
n
j
j
j
i
n
iii
简记即为阶梯（逼近）函数样条插值公式阶则对给定的型值点
样条插值函数。显然他就是一阶得到一次磨光函数进行一次磨光对阶梯函数
b
.)(
)(
1
)(
1
)(
)(
1
)(
0
1
0
2/
2/
0
2/
2/
01
2/
2/
0
0
n
j
j
j
n
j
hx
hx
j
j
hx
hx
hx
hx
h
xx
y
dt
h
xt
y
h
dttf
h
xf
dttf
h
xf
n
j
j
kjk h
xx
yxf
kxf
0
0
).()(
)( 次磨光（逼近）函数为的以此类推，
的逼近性质。、下面着重分析 )()( 32 xfxf
定理 1：
2
2/112/12
2
2/112/12
2
)4(
11
2
2
22
112
2/12
24
)(
)()(
1
)(
8
)(
)()(
2
1
)(
)3
8
)(
)()2(
1
)(
)2
)(
8
)(y
)2(
8
1
)(
1
)(
h
y
xyyy
h
xf
h
y
xyyyxf
h
y
xyyyy
h
xf
hohyyyyxf
xxxf
iiii
iiii
iiiii
i
iiiii
ii
、相切性质
、凹凸性质
）、逼近性质处有：及相邻节点的中点在节点
)33(
8
1
)(
1
)(
)(
48
1
)(
2
1
)(
)3
8
)(
)()2(
1
)(
)2
)(
6
)(y
)2(
6
1
)(
1
)(2
11212/13
21112/13
2
)4(
11
2
3
22
113
2/13
iiiiiii
iiiiiii
iiiii
i
iiiii
ii
yyyy
h
yy
h
xf
yyyyyyxf
h
y
xyyyy
h
xf
hohyyyyxf
xxxf
、相切性质
、凹凸性质
）、逼近性质处有：及相邻节点的中点在节点：定理
磨光法在外形设计问题重的应用例 某外形设计问题所给的型值 如下表。要求寻找这样的曲线，满足：
1）通过两个端点，
2）在 x(i)处逼近原型值的绝对误差不超过 1厘米，
3）在 x(i)处的二阶导数符号与原型值点的二阶差符号相同。
1112 2 ), iiiiii yyyyyx 的二阶差注：型值（
X(i)
cm
30 80 130 180 230 280 330 380 430 480
Y(i)
cm
80 110 132 148.75 163 175 185.5 195 204 212.75
02
02
1
,,,
)()(
).,(),,(
3)(
3)1:1()(2
.9,50,30
)
h
xx
(Ωy( x )f
3
11
101
11101
1
1
33
1111
03
3
0
n
0i
i
3i3
nnn
nnn
n
i
i
i
nn
n
i
yyy
yy y
yyhxxhxx
h
xx
yxf
yxyx
xxxf
nixxf
nhx
）得到：由条件待定。
此时三次磨光函数为两侧延拓，即增加右），必须将型值点向左处满足条件、在边界点要使
）。满足条件在内节点保证了显然，定理磨光函数）设计曲线可采用三次解：由条件
）。不满足条件直接计算表明由
）：现在我们来验证条件即：
2
)2(
6
1
)(
2
.2
,2
113
11
101
iiiii
nnn
yyyyxf
yyy
yyy
下面我们采用盈亏修正法来提高精度。其算法如下：
)(?)(?
.2?,2?,:0,?
)0
3
1
1
3
11101
h
xx
yxf
yyyyyyniyy
i
n
i
i
nnnii
开始（初始化）
）。转步形成新的三次磨光函数修正型值下一步。，停止计算；否则执行
）计算误差
1)4
)(?)(
)3
.2?
,2?
,:0,
)2
:0
)4?(
6
1
)(
1
3
1
1
3
11
101
113
h
xx
yxf
yyy
yyy
niryy
rif
ni
yyyyxfyr
i
n
i
i
nnn
iii
iiiiiii
实例计算结果
§ 3.4 数据拟合的最小二乘法问题的引入先讨论两个变量 x,y的情况。通过观测变量 x,y积累了一组资料（ x(i),y(i)),i=1:n，一般来说 n比较大。我们的任务是从这批实验数据出发，寻求一近似函数 φ(x)曲逼近 y。
由于观测数据都带有观测误差，数目又较大，对于这类问题运用插值函数取描述 y往往是不适当的。现在，用开始引入的例 1来说明建立近似函数的一种方法，即最小二乘法。
研究 y与 t的关系的最小二乘法，通常步骤如下：
1）先描 t,y的草图
2）凭视觉猜想一数学模型
t
t
t
tt
etb
tctccta
t
1
)(
1
ln))(ln ()(
)(:)(
)(:)(
/
2
210
，即本例可用两种试验模型
3）计算总误差
n
i
ii tycccI
a
1
2
210 ))((),,(
)(
来说是就试验模型
4）希望总误差最小
m in),,( 210?cccI
.2,1,0,0
m in),,( 210
j
c
I
cccI
j
m i n)),,,((),,,(
,),,,(,
:1),,,,,(
m,,,
1
2
2121
1
21
21
21
m
j
jnjjjn
n
i
iin
jjnjj
n
xxxycccI
xcxxx
mjyxxx
xxxy
则有若选择模型次观测得到观测值：有关系，并通过与多个变量如果数据拟合的最小二乘法
.:1,)()()(
:1,0 m i n)(
,,,
)(
.,,,m i n,)(
,)()(,:1,,
11 1
1
2
21
1
2
21
1
2
1
mkxycxx
mk
c
I
xyI
ccc
xy
cccxyI
xcxniyx
n
i
ikii
m
j
j
n
i
ikiji
k
n
i
iii
m
n
i
iii
m
n
i
iii
m
j
jjii
采用：确定参数越小，模型越好。
为均方误差。均方误差当参数确定后，称确定参数按最小二乘优化准则确定试验模型：根据观测数据
mmmmmm
m
m
n
i
ikiik
T
kk
n
i
ikijij
T
kjkkj
T
n
T
njjjj
b
b
b
c
c
c
aaa
aaa
aaa
mjk
xyyyb
xxa
yy
xxx
2
1
2
1
21
22221
11211
1
1
21
21
:1,
,)(,
,)()(),(
,,,,yy
,)(,),(),(
则令
对称矩阵方程例 1的数据拟合结果其中解：转化为采用平方误差：
从而拟合函数为解出为：正规方程采用
,
)()2(
9486.3
.04832 0.01436.11490.4)(
),04832 0.0,1436.1,1490.4(
01.8530
59.757
49.88
14043 410396826
1039682676
8267610
)()1(
/
2
2
2
1
0
2
210
bxaz
et
ttt
c
c
c
c
bAc
tctcct
t
.,ln,ln,1 baztx
);3 6 1 8.2,3 5 4 2.2,3 4 3 7.2,3 3 5 1.2,2 8 8 5.2,2 5 4 4.2,1 7 3 6.2,0 8 0 7.2,8 5 7 9.1,3 8 6 3.1(;0 6 2 5.0,0 7 1 4.0,0 8 3 3.0,1 0 0.0,1 2 5 0.0,1 6 6 7.0,2 5 0 0.0,3 3 3 3.0,5 0 0.0,0 0 0.1
,
z
x
zx 的数据对应：根据已知数据得到
1109.0
3411.11)(
.0579.1,3411.11,0579.1,4284.2
9586.4
4362.21
4930.16923.2
6923.210
2
/0 5 7 9.1
平方误差所求拟合函数为从而解出为正规方程
t
a
et
eba
b
a
bAc
显然，比第一种模型要好得多！
）的分布。和（拟合模型下图显示出原始数据,2)1(
