数值计算方法程序设计Word文档下载推荐.docx
《数值计算方法程序设计Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《数值计算方法程序设计Word文档下载推荐.docx（19页珍藏版）》请在冰豆网上搜索。
%首先说明：
追赶法是适用于三对角矩阵的线性方程组求解的方法，并不适用于其他类型矩阵。
%定义三对角矩阵A的各组成单元。
方程为Ax=d
%b为A的对角线元素（1~n），a为-1对角线元素（2~n），c为+1对角线元素（1~n-1）。
%A=[2-100
%-13-20
%0-24-3
%00-35]
a=[0-1-2-3];
c=[-1-2-3];
b=[2345];
d=[61-21];
n=length（b）;
u0=0;
y0=0;
a
（1）=0;
%“追”的过程
L
（1）=b
（1）-a
（1）*u0;
y
（1）=（d
（1）-y0*a
（1））/L
（1）;
u
（1）=c
（1）/L
（1）;
fori=2:
（n-1）
L（i）=b（i）-a（i）*u（i-1）;
y（i）=（d（i）-y（i-1）*a（i））/L（i）;
u（i）=c（i）/L（i）;
end
L（n）=b（n）-a（n）*u（n-1）;
y（n）=（d（n）-y（n-1）*a（n））/L（n）;
%“赶”的过程
x（n）=y（n）;
fori=（n-1）:
-1:
1
x（i）=y（i）-u（i）*x（i+1）;
3.特征向量的计算，幂法
5.2.2幂法的MATLAB程序
用幂法计算矩阵的主特征值和对应的特征向量的MATLAB主程序
function[k,lambda,Vk,Wc]=mifa（A,V0,jd,max1）
lambda=0;
k=1;
Wc=1;
jd=jd*;
state=1;
V=V0;
while（（k<
=max1）&
（state==1））
Vk=A*V;
[mj]=max（abs（Vk））;
mk=m;
tzw=abs（lambda-mk）;
Vk=（1/mk）*Vk;
Txw=norm（V-Vk）;
Wc=max（Txw,tzw）;
V=Vk;
lambda=mk;
state=0;
if（Wc>
jd）
k=k+1;
Wc=Wc;
if（Wc<
=jd）
请注意：
迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：
'
）
else
迭代次数k已经达到最大迭代次数max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc如下：
Vk=V;
k=k-1;
Wc;
例5.2.2用幂法计算下列矩阵的主特征值和对应的特征向量的近似向量，精度.并把
（1）和
（2）输出的结果与例中的结果进行比较.
（1）；
（2）；
（3）；
（4）.
解
（1）输入MATLAB程序
>
A=[1-1;
24];
V0=[1,1]'
;
[k,lambda,Vk,Wc]=mifa（A,V0,,100）,
[V,D]=eig（A）,Dzd=max（diag（D））,wuD=abs（Dzd-lambda）,wuV=V（:
2）./Vk,
运行后屏幕显示结果
k=lambda=Wc=
33
Vk=V=wuV=
Dzd=wuD=
3
由输出结果可看出，迭代33次，相邻两次迭代的误差Wc19e-007，矩阵的主特征值的近似值00和对应的特征向量的近似向量Vk（00，00,lambda与例5.1.1中的最大特征值近似相等，绝对误差约为37e-006，Vk与特征向量的第1个分量的绝对误差约等于0，第2个分量的绝对值相同.由wuV可以看出，的特征向量V（:
2）与Vk的对应分量的比值近似相等.因此，用程序计算的结果达到预先给定的精度.
（2）输入MATLAB程序
B=[123;
213;
336];
V0=[1,1,1]'
[k,lambda,Vk,Wc]=mifa（B,V0,,100）,[V,D]=eig（B）,
Dzd=max（diag（D））,wuD=abs（Dzd-lambda）,wuV=V（:
3）./Vk,
k=lambda=Wc=Dzd=wuD=
39090
Vk=wuV=
V=
0
（3）输入MATLAB程序
C=[122;
1-11;
4-121];
V0=[1,1,1]'
[k,lambda,Vk,Wc]=mifa（C,V0,,100）,[V,D]=eig（C）,
Dzd=max（diag（D））,wuD=abs（Dzd-lambda）,
Vzd=V（:
1）,wuV=V（:
1）./Vk,
运行后屏幕显示
k=lambda=Wc=
100
Dzd=wuD=
Vk=Vzd=wuV=
由输出结果可见，迭代次数k已经达到最大迭代次数max1=100，并且lambda的相邻两次迭代的误差58>
2,由wuV可以看出，lambda的特征向量Vk与真值Dzd的特征向量Vzd对应分量的比值相差较大，所以迭代序列发散.实际上，实数矩阵C的特征值的近似值为，并且对应的特征向量的近似向量分别为=（，，），
（,），
（,+,
+,是常数）.
（4）输入MATLAB程序
D=[-4140;
-5130;
-102];
[k,lambda,Vk,Wc]=mifa（D,V0,,100）,[V,Dt]=eig（D）,
Dtzd=max（diag（Dt））,wuDt=abs（Dtzd-lambda）,
Vzd=V（:
2）,wuV=V（:
2）./Vk,
19
Dtzd=wuDt=
Vk=Vzd=wuV=
（一）原点位移反幂法的MATLAB主程序1
用原点位移反幂法计算矩阵的特征值和对应的特征向量的MATLAB主程序1
function[k,lambdan,Vk,Wc]=ydwyfmf（A,V0,jlamb,jd,max1）
[n,n]=size（A）;
A1=A-jlamb*eye（n）;
jd=jd*;
RA1=det（A1）;
ifRA1==0
因为A-aE的n阶行列式hl等于零，所以A-aE不能进行LU分解.'
）
return
ifRA1~=0
forp=1:
n
h（p）=det（A1（1:
p,1:
p））;
hl=h（1:
n）;
fori=1:
ifh（1,i）==0
因为A-aE的r阶主子式等于零，所以A-aE不能进行LU分解.'
ifh（1,i）~=0
因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.'
Vk=V0;
[LU]=lu（A1）;
Yk=L\Vk;
Vk=U\Yk;
mk=m;
Vk1=Vk/mk;
Yk1=L\Vk1;
Vk1=U\Yk1;
[mj]=max（abs（Vk1））;
mk1=m;
Vk2=（1/mk1）*Vk1;
tzw1=abs（（mk-mk1）/mk1）;
tzw2=abs（mk1-mk）;
Txw1=norm（Vk）-norm（Vk1）;
Txw2=（norm（Vk）-norm（Vk1））/norm（Vk1）;
Txw=min（Txw1,Txw2）;
tzw=min（tzw1,tzw2）;
Vk=Vk2;
mk=mk1;
%Vk=Vk2,mk=mk1,
A-aE的秩R（A-aE）和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：
A-aE的秩R（A-aE）和各阶顺序主子式值hl、迭代次数k已经达到最大迭代次数max1,按模最小特征值的迭代值lambda,特征向量的迭代向量Vk,相邻两次迭代的误差Wc如下：
hl,RA1
[V,D]=eig（A,'
nobalance'
）,Vk;
lambdan=jlamb+1/mk1;
例5.3.2用原点位移反幂法的迭代公式（），根据给定的下列矩阵的特征值的初始值，计算与对应的特征向量的近似向量,精确到1.
（1）,；
（2），；
（3），.
A=[1-10;
-24-2;
0-12];
[k,lambda,Vk,Wc]=ydwyfmf（A,V0,,,10000）
因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.
k=lambda=Wc=hl=
3
Vk=V=D=
00
00
（2）输入MATLAB程序
A=[1-1;
V0=[20,1]'
[k,lambda,Vk,Wc]=ydwyfmf（A,V0,,,100）
因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU
升级会员 