矩阵和数值分析实验作业.docx
《矩阵和数值分析实验作业.docx》由会员分享，可在线阅读，更多相关《矩阵和数值分析实验作业.docx（36页珍藏版）》请在冰豆网上搜索。
矩阵和数值分析实验作业
一．求累加值
原题：
设
，分别编制从小到大和从大到小的顺序程序计算
并指出有效位数。
解：
1.从小到大
a.程序代码：
lejia1.m
formatlong
N=input（'请输入N='）;
s=0;
forj=2:
N
s=s+1/（j*j-1）;
end%求近似值
s1=s*10^（5）;
s2=fix（s1）;%取五位数字
S=10^（-5）*s2;
a=0.5*（1.5-1/N-1/（N+1））;
y=abs（S-a）;
fori=0:
100
ify>=1/2*10^（-i）
break
end%判断有效数字的位数
end
disp（'结果是:
'）;
S
disp（'有效数字的位数是:
'）;
i-1
b.运行结果：
2．从大到小
a.程序代码：
leijia2.m
formatlong
N=input（'pleaseinputN:
'）;
s=0;
forj=2:
N
s=s+1/[（N-j+2）^2-1];
end%求近似值
s1=s*10^（5）;
s2=fix（s1）;
S=10^（-5）*s2;%取五位数字
a=0.5*（1.5-1/N-1/（N+1））;
y=abs（S-a）;
fori=0:
100
ify>=1/2*10^（-i）
break
end%判断有效数字位数
end
disp（'结果是:
'）;
S
disp（'有效数字的位数是:
'）;
i-1
b.运行结果：
二．解线性方程组
1.原题：
分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组
迭代法计算停止条件为：
。
解：
1>.Jacobi迭代法：
a.程序代码：
Jacobi.m
A=input（'请输入矩阵A='）;
b=input（'请输入列向量b='）;
X=input（'请输入初始值X='）;
S=input（'请输入终止条件S='）;
max1=100;%最大循环次数
B=diag（A）;
G=diag（B）;%求对角阵
D=inv（G）;%求对角阵的逆矩阵
P=G-A;
fork=1:
max1
Y=D*P*X+D*b;
Z=abs（Y-X）;
[t]=max（Z）;
X=Y;
Iftbreak
end
end%计算方程的解
X
b.运行结果：
2>.G-S迭代法：
a.程序代码：
GS.m
A=input（'请输入矩阵A='）;
b=input（'请输入列向量b='）;
X=input（'请输入初始值X='）;
S=input（'请输入终止条件S='）;
max1=100;%最大循环次数
C=diag（A）;
P=diag（C）;%求A的对角阵
U=-triu（A）+P;%求A的下三角阵
G=A+U;
D=inv（G）;
fork=1:
max1
Y=D*U*X+D*b;
Z=abs（Y-X）;
[t]=max（Z）;
X=Y;
if（tbreak
end
end
X
b.运行结果：
2.原题：
用Gauss列主元消去法、QR方法求解如下方程组：
解：
1>.Gauss列主元消去法：
a.程序代码：
GAUSS.m
A=input（'请输入系数矩阵A='）;
b=input（'请输入列向量b='）;
G=[A,b];
B=[];%定义一个空矩阵
t=size（b）;%求系数矩阵的维数
k=t（1,1）;%求系数矩阵的维数
fori=1:
k
C=0;
fore=i:
k
C（e）=abs（G（e,i））;
end
l=max（C）;%挑选最大值
l2=find（C==l）;
l1=l2（1,1）;
T=G（l1,:
）;
G（l1,:
）=G（i,:
）;
G（i,:
）=T;%选主元，交换最大值所在的行
disp（'主元调整后的矩阵：
'）;
G
forj=i+1:
k
P=-G（j,i）/G（i,i）;
B（j,:
）=P*G（i,:
）+G（j,:
）;
G（j,:
）=B（j,:
）;
end%进行消元
disp（'消元后的矩阵：
'）;
G
end
b=G（:
k+1）;
G（:
k+1）=[];
X=[];
form=1:
k
X（k-m+1）=b（k-m+1）/G（k-m+1,k-m+1）;
b=b-X（k-m+1）*G（:
k-m+1）;
end%求方程组的解
disp（'方程组的解：
'）;
X
b.运行结果：
2>.QR法：
a.程序代码：
QR.m
A=input（'请输入系数矩阵A='）;
b=input（'请输入列向量b='）;
m=size（A）;
k=m（1,1）;%计算A的维数
B=A;
M=eye（size（A））;%定义一个与A同维数的单位矩阵
Q=M;
forr=1:
k-1;
I=eye（size（B））;
e=I（:
1）;
a=B（:
1）;%取出B的第一列
d=norm（a）;
w=a-d*e;
c=rot90（w）;
H=I-2*w*c/（c*w）;%求变换矩阵
P=H*B;
B=P;
B（1,:
）=[];
B（:
1）=[];
fori=1:
k-r+1
forj=1:
k-r+1
M（r+i-1,r+j-1）=H（i,j）;
end
end%将变换矩阵H变成k维
Z=M;
M=eye（size（A））;
Q=Z*Q;
end
Q
R=Q*A;
R
X=[];
b=Q*b;
foru=1:
k
X（k-u+1）=b（k-u+1）/R（k-u+1,k-u+1）;
b=b-X（k-u+1）*R（:
k-u+1）;
end%计算方程组的解
disp（'方程组的解：
'）;
X
b.运行结果：
三．非线性方程的迭代解法
1.原题：
用Newton迭代法求解方程
的根，计算停止的条件为：
。
解：
a.程序代码：
newton.m
symsx;
f=input（'请以x为变量输入函数f='）;
s=input（'请输入终止条件s='）;
X=input（'请输入初始值X='）;
max1=100;%最大循环次数
p=diff（f）;%求f的导函数
fork=1:
max1
y=X-subs（f,X）/subs（p,X）;%迭代公式
m=abs（y-X）;
X=y;
if（mbreak
end
end
disp（'方程的根：
'）;
X
a.运行结果：
2．原题：
利用Newton迭代法求多项式
的所有实零点，注意重根的问题。
解：
a.程序代码：
new.m
functionX=newton1（X,s）
symsx;
f=x^4-5.4*x^3+10.56*x^2-8.954*x+2.7951;
max1=20;
p=diff（f）;
fork=1:
max1
y=X-subs（f,X）/subs（p,X）;
m=abs（y-X）;
X=y;
if（mbreak
end
end%定义一个newton迭代法求线性方程组的的函数
symsx;
f=input（'请以x为变量输入多项式f='）;
c=input（'请输入多项式的次数c='）;
a=input（'请输入根所在的估计的最大区间的左端点a='）;
b=input（'请输入根所在的估计的最大区间的右端点b='）;
n=input（'请输入将[a,b]区间分成的个数n='）;
s=input（'请输入newton搜索时的终止条件s='）;
A=a:
（b-a）/n:
b;%将区间n等分
disp（'n等分后的区间为:
'）;
A
fori=1:
n
Y（i）=subs（f,A（i））;
Y（i+1）=subs（f,A（i+1））;
ifY（i）*Y（i+1）<0
M（i）=newton1（A（i）,s）;
elseifabs（Y（i））==0
M（i）=A（i）;%求多项式的根
end
end
B=find（M~=0）;
k1=size（B）;
k=k1（1,2）;
forkk=1:
k
X（kk）=M（B（kk））;%将根以矩阵形式表示
end
disp（'多项式的根为:
'）;
X
p=0;
c1=c-k;
forj=1:
k
q=f;
fore=1:
c1+1
q=diff（q,x）;
y=subs（q,X（j））;
ifabs（y）>=100*s
break
end%判断根的重数
end
ife>1
disp（'根的重数：
'）;
X（j）
disp（'的重数：
'）;
disp（e）;
p=p+e-1;%对根重数大于1的重数进行累加
end
ifp==k-c
break%判断重数是否已达到总重数
end
end
b.运行结果：
四．数值积分
原题：
分别用复化梯形公式和Romberg公式计算积分
要求误差不超过10-5，并给出节点个数。
解：
1>.复化梯形公式：
a.程序代码：
fuhuatixing.m
formatlong;
symsx;
f=input（'请以x为自变量输入函数f='）;
a=input（'请输入积分区域左端点a='）;
b=input（'请输入积分区域右端点b='）;
s=input（'请输入终止条件s='）;
n=100;
X=a:
（b-a）/n:
b;
t1=0;
fa=subs（f,a）;
fb=subs（f,b）;
fori=2:
n
t1=t1+subs（f,X（i））;
T（i）=（b-a）*（fa+fb+2*t1）/（2*n）;%用复化梯形公式计算积分结果
if（T（i）-T（i-1））<=s
break
end
end
disp（'积分结果:
'）;
T（i）
b.运行结果：
2>.Romberg公式：
a.程序代码：
romberg.m
symsx;
f=input（'请以x为自变量输入函数f='）;
a=input（'请输入积分区域左端点a='）;
b=input（'请输入积分区域右端点b='）;
s=input（'请输入终止条件s='）;
T（1,1）=1/2*（subs（f,a）+subs（f,b））;
n=100;
fori=2:
n
q1=0;
fore=1:
（2^（i-2））
q1=q1+subs（f,（a+（2*（e-1）+1）*（b-a）/2^（i-1）））;
end
q=（b-a）*q1/2^（i-1）;
T（i,1）=1/2*T（i-1,1）+q;%计算复化梯形得出的值
forj=2:
i
T（i,j）=（4^（j-1）*T（i,j-1）-T（i-1,j-1））/（4^（j-1）-1）;%计算n大于2时复化newton-cotes的值
end
ifabs（T（i,i）-T（i-1,i-1））break
end
end
disp（'T数表:
'）;
T
disp（'积分结果:
'）;
T（i,i）
b.运行结果：
升级会员 