大连理工大学矩阵大作业汇编.docx
《大连理工大学矩阵大作业汇编.docx》由会员分享，可在线阅读，更多相关《大连理工大学矩阵大作业汇编.docx（35页珍藏版）》请在冰点文库上搜索。
大连理工大学矩阵大作业汇编
（1）从大到小的顺序的计算程序：
functiony=snd（n）
formatlong
y=0;
ifn<2
disp（'请输入大于1的数！
'）
else
s=0;
i=2;
whilei<=n
s=single（s+（1/（i^2-1）））;
i=i+1;
end
y=s;
end
（2）从小到大的顺序的计算程序：
functiony=snx（n）
formatlong
y=0;
ifn<2
disp（'请输入大于1的数！
'）
else
s=0;
i=n;
while1
s=single（s+（1/（i^2-1）））;
i=i-1;
ifi==1
break
end
end
y=s;
end
（3）按两种顺序分别计算
并指出有效位数（编制程序时用单精度）
1
的计算结果：
2
的计算结果：
3
的计算结果：
计算时的有效位数为七位数。
1 秦九昭算法计算程序：
functiony=qjz（a,x）
j=3;
i=size（a,2）;
switchi
case1
y=a
（1）;
case2
y=a
（1）*x+a
（2）;
otherwise
p=a
（1）*x+a
（2）;
whilej<=i
p=p*x+a（j）;
j=j+1;
end
y=p;
end
2 计算
在点23的值。
计算结果如下：
当
时
。
1 Gauss法计算程序和结果：
程序：
A（1,:
）=[31,-13,0,0,0,-10,0,0,0];
A（2,:
）=[-13,35,-9,0,-11,0,0,0,0];
A（3,:
）=[0,-9,31,-10,0,0,0,0,0];
A（4,:
）=[0,0,-10,79,-30,0,0,0,-9];
A（5,:
）=[0,0,0,-30,57,-7,0,-5,0];
A（6,:
）=[0,0,0,0,-7,47,-30,0,0];
A（7,:
）=[0,0,0,0,0,-30,41,0,0];
A（8,:
）=[0,0,0,0,-5,0,0,27,-2];
A（9,:
）=[0,0,0,-9,0,0,0,-2,29];
B=[-15;27;-23;0;-20;12;-7;7;10];
[a,b]=gauss（A,B）;
j=size（a,2）;
whilej>=1
k=j+1;
s=b（j）;
whilek<=9
s=s-x（k）*a（j,k）;
k=k+1;
end
x（j）=s/a（j,j）;
j=j-1;
end
disp（x）
function[x,y]=gauss（a,b）
num_i=size（a,1）;
j=1;
whilej<=（num_i-1）
i=j+1;
whilei<=num_i
r=a（i,j）/a（j,j）;
a（i,:
）=a（i,:
）-r*a（j,:
）;
b（i,:
）=b（i,:
）-r*b（j,:
）;
i=i+1;
end
j=j+1;
end
x=a;
y=b;
运行的结果为：
。
2 列主元消去法计算程序和结果：
A（1,:
）=[31,-13,0,0,0,-10,0,0,0];
A（2,:
）=[-13,35,-9,0,-11,0,0,0,0];
A（3,:
）=[0,-9,31,-10,0,0,0,0,0];
A（4,:
）=[0,0,-10,79,-30,0,0,0,-9];
A（5,:
）=[0,0,0,-30,57,-7,0,-5,0];
A（6,:
）=[0,0,0,0,-7,47,-30,0,0];
A（7,:
）=[0,0,0,0,0,-30,41,0,0];
A（8,:
）=[0,0,0,0,-5,0,0,27,-2];
A（9,:
）=[0,0,0,-9,0,0,0,-2,29];
B=[-15;27;-23;0;-20;12;-7;7;10];
[a,b]=lzy（A,B）;
j=size（a,2）;
whilej>=1
k=j+1;
s=b（j）;
whilek<=9
s=s-x（k）*a（j,k）;
k=k+1;
end
x（j）=s/a（j,j）;
j=j-1;
end
disp（x）
function[a1,b1]=lzy（a,b）
[num_i,num_j]=size（a）;
ab=zeros（num_i,num_j+1）;
fork=1:
num_j
ab（:
k）=a（:
k）;
end
ab（:
num_j+1）=b（:
1）;
j=1;
whilej[max,max_i]=searmax（j,j,ab）;
i_nu=ab（j,:
）;
ab（j,:
）=ab（max_i,:
）;
ab（max_i,:
）=i_nu;
m=j+1;
whilem<=num_i
forn=j:
num_j+1
ab（m,n）=ab（m,n）-（ab（m,j）/max）*ab（j,n）;
end
m=m+1;
end
j=j+1;
end
a1=zeros（num_i,num_j）;
fork=1:
num_i
a1（:
k）=ab（:
k）;
end
b1=ab（:
num_i+1）;
function[b,c]=searmax（i,j,a）
num_i=size（a,1）;
k=i;
m=abs（a（k,j））;
c=k;
whilekk=k+1;
ifm>=abs（a（k,j））
continue
else
m=abs（a（k,j））;
c=k;
end
end
b=a（c,j）;
运行的结果为：
（1）LU分解的计算程序及结果：
function[l,u]=lufz（a）
[num_i,num_j]=size（a）;
l=eye（num_i）;
u=eye（num_i）;
forj=1:
num_j
u（1,j）=a（1,j）;
l（j,1）=a（j,1）/u（1,1）;
end
i=2;
whilei<=num_i
j=i;
whilejs=0;
fork=1:
i-1
s=s+l（i,k）*u（k,j）;
end
u（i,j）=a（i,j）-s;
s=0;
fork=1:
i-1
s=s+l（j+1,k）*u（k,i）;
end
l（j+1,i）=（a（j+1,i）-s）/u（i,i）;
j=j+1;
end
s=0;
fork=1:
i-1
s=s+l（i,k）*u（k,num_i）;
end
u（i,num_i）=a（i,num_i）-s;
i=i+1;
end
输入要求解的矩阵后运行的结果为：
（2）带列主元的LU分解计算程序及结果
由于Matlab中自带带列主元的LU分解函数，故计算程序如下：
A（1,:
）=[31,-13,0,0,0,-10,0,0,0];
A（2,:
）=[-13,35,-9,0,-11,0,0,0,0];
A（3,:
）=[0,-9,31,-10,0,0,0,0,0];
A（4,:
）=[0,0,-10,79,-30,0,0,0,-9];
A（5,:
）=[0,0,0,-30,57,-7,0,-5,0];
A（6,:
）=[0,0,0,0,-7,47,-30,0,0];
A（7,:
）=[0,0,0,0,0,-30,41,0,0];
A（8,:
）=[0,0,0,0,-5,0,0,27,-2];
A（9,:
）=[0,0,0,-9,0,0,0,-2,29];
[L,U,P]=lu（A）;
运行结果如下：
为单位阵。
（3）求A的逆矩阵的计算程序及结果：
A（1,:
）=[31,-13,0,0,0,-10,0,0,0];
A（2,:
）=[-13,35,-9,0,-11,0,0,0,0];
A（3,:
）=[0,-9,31,-10,0,0,0,0,0];
A（4,:
）=[0,0,-10,79,-30,0,0,0,-9];
A（5,:
）=[0,0,0,-30,57,-7,0,-5,0];
A（6,:
）=[0,0,0,0,-7,47,-30,0,0];
A（7,:
）=[0,0,0,0,0,-30,41,0,0];
A（8,:
）=[0,0,0,0,-5,0,0,27,-2];
A（9,:
）=[0,0,0,-9,0,0,0,-2,29];
[num_i,num_j]=size（A）;
A_n=eye（num_i）;
E=eye（num_i）;
[L,U]=lufz（A）;
fori=1:
num_i
y=hd（1,L,E（:
i））;
A_n（:
i）=hd（0,U,y）;
end
disp（A_n）
function[l,u]=lufz（a）
[num_i,num_j]=size（a）;
l=eye（num_i）;
u=eye（num_i）;
forj=1:
num_j
u（1,j）=a（1,j）;
l（j,1）=a（j,1）/u（1,1）;
end
i=2;
whilei<=num_i
j=i;
whilejs=0;
fork=1:
i-1
s=s+l（i,k）*u（k,j）;
end
u（i,j）=a（i,j）-s;
s=0;
fork=1:
i-1
s=s+l（j+1,k）*u（k,i）;
end
l（j+1,i）=（a（j+1,i）-s）/u（i,i）;
j=j+1;
end
s=0;
fork=1:
i-1
s=s+l（i,k）*u（k,num_i）;
end
u（i,num_i）=a（i,num_i）-s;
i=i+1;
end
functionx=hd（f,a,b）
[num_i,num_j]=size（a）;
x=zeros（num_i,1）;
switchf
case0
i=num_i-1;
x（num_i）=b（num_i）/a（num_i,num_j）;
whilei>=1
s=0;
fork=i+1:
num_i
s=s+a（i,k）*x（k）;
end
x（i）=（b（i）-s）/a（i,i）;
i=i-1;
end
case1
x
（1）=b
（1）/a（1,1）;
j=2;
whilej<=num_i
s=0;
fork=1:
j-1
s=s+a（j,k）*x（k）;
end
x（j）=（b（j）-s）/a（j,j）;
j=j+1;
end
otherwise
disp（'error！
请输入正确的参数！
'）
end
运行结果：
（4）求A的行列式的计算程序及结果：
A（1,:
）=[31,-13,0,0,0,-10,0,0,0];
A（2,:
）=[-13,35,-9,0,-11,0,0,0,0];
A（3,:
）=[0,-9,31,-10,0,0,0,0,0];
A（4,:
）=[0,0,-10,79,-30,0,0,0,-9];
A（5,:
）=[0,0,0,-30,57,-7,0,-5,0];
A（6,:
）=[0,0,0,0,-7,47,-30,0,0];
A（7,:
）=[0,0,0,0,0,-30,41,0,0];
A（8,:
）=[0,0,0,0,-5,0,0,27,-2];
A（9,:
）=[0,0,0,-9,0,0,0,-2,29];
[num_i,num_j]=size（A）;
[L,U]=lufz（A）;
s=1;
fori=1:
num_i
s=s*U（i,i）;
end
disp（['矩阵的行列式值为：
',num2str（s）]）
运行程序后结果与调用matlab中det（）函数结果如下：
（1）求cholesky分解程序及结果：
functionl=chole（a）
[num_i,num_j]=size（a）;
l=zeros（num_i）;
l（1,1）=sqrt（a（1,1））;
fork=2:
num_i
l（k,1）=a（k,1）/l（1,1）;
end
j=2;
whilej<=num_j
s1=0;
fork=1:
j-1
s1=s1+l（j,k）^2;
end
l（j,j）=sqrt（a（j,j）-s1）;
i=j+1;
whilei<=num_i
s2=0;
fork=1:
j-1
s2=s2+l（i,k）*l（j,k）;
end
l（i,j）=（a（i,j）-s2）/l（j,j）;
i=i+1;
end
j=j+1;
end
运行程序后的结果为：
（2）求解方程组程序及结果：
a=[2,1,-1,1;1,5,2,7;-1,2,10,1;1,7,1,11];
b=[13;-9;6;0];
L=chole（a）;
y=hd（1,L,b）;
x=hd（0,L',y）;
disp（x）
运行后的结果为:
程序和结果如下：
A=[4,2,0,0;3,-2,1,0;0,2,5,3;0,0,-1,6];
B=[6;2;10;5];
[L,U]=sjfj（A）;%
y=hd（1,L,B）;
x=hd（0,U,y）;
disp（'方程组的解为：
'）
disp（x）
function[l,u]=sjfj（a）
num_i=size（a,1）;
i=2;
whilei<=num_i
a（i,i-1）=a（i,i-1）/a（i-1,i-1）;
a（i,i）=a（i,i）-a（i,i-1）*a（i-1,i）;
i=i+1;
end
l=tril（a）;
u=triu（a）;
forj=1:
num_i;
l（j,j）=1;
end
运行程序后结果为：
编程求解矩阵A的QR分解：
（1）QR分解计算程序：
function[Q,R]=hqr（A）
[n,n]=size（A）;
Q=eye（n）;
fori=1:
（n-1）
E=eye（n-i+1）;
e1=E（:
1）;
a=zeros（1,n-i+1）';
forj=1:
（n-i+1）
a（j）=A（j+i-1,i）;
end
av=sqrt（a'*a）;
w=a-av*e1;
h=E-（2/（w'*w））*（w*w'）;
q=eye（n）;
fork=i:
n
forj=i:
n
q（k,j）=h（k-i+1,j-i+1）;
end
end
A=q*A;
Q=q*Q;
end
R=A;
Q=Q';
（2）输入矩阵A后的计算结果：
（1）Jacobi迭代法求解程序与结果：
a=0;
b=0;
c=0;
while1
a0=0.25*（7+b-c）;
b0=（1/8）*（-21-4*a-c）;
c0=（1/5）*（15+2*a-b）;
A=[abs（a-a0）,abs（b-b0）,abs（c-c0）];
f=max（A）;
iff<=0.001
x=[a0,b0,c0];
break
else
a=a0;
b=b0;
c=c0;
end
end
disp（'方程组的解为x'）
disp（x）
运行后的结果为：
（2）Gauss-Seidel迭代法求解的程序与结果：
a=0;
b=0;
c=0;
while1
a0=a;
b0=b;
c0=c;
a=（1/4）*（7+b-c）;
b=（1/8）*（-21-4*a-c）;
c=（1/5）*（15+2*a-b）;
A=[abs（a-a0）,abs（b-b0）,abs（c-c0）];
f=max（A）;
iff<=0.001
x=[a,b,c];
break
end
end
disp（'方程组的解为：
'）
disp（x）
运行后的结果为：
（1）计算
在
上的根的程序：
symsx
f=2*x^3-5*x-1;
df=diff（f,x）;
g=x-（f/df）;
x0=1;
while1
x1=subs（g,x0）;
dx=abs（x1-x0）;
ifdx<=0.001
disp（['Newton法的解为：
x=',num2str（x1）]）
break
else
x0=x1;
end
end
x0=1;
x1=1.1;
while1
x2=x1-（subs（f,x1）/（subs（f,x1）-subs（f,x0）））*（x1-x0）;
dx=abs（x2-x1）;
ifdx<=0.001
disp（['割线法的解为：
x=',num2str（x2）]）
break
else
x0=x1;
x1=x2;
end
end
运行后结果为：
Newton法
，割线法
。
（2）计算
在
上的根的程序：
symsx
f=exp（x）*sin（x）;
df=diff（f,x）;
g=x-（f/df）;
x0=-2;
while1
x1=subs（g,x0）;
dx=abs（x1-x0）;
ifdx<=0.001
disp（['Newton法的解为：
x=',num2str（x1）]）
break
else
x0=x1;
end
end
x0=-2;
x1=-2.1;
while1
x2=x1-（subs（f,x1）/（subs（f,x1）-subs（f,x0）））*（x1-x0）;
dx=abs（x2-x1）;
ifdx<=0.001
disp（['割线法的解为：
x=',num2str（x2）]）
break
else
x0=x1;
x1=x2;
end
end
运行后结果为：
Newton法
，割线法
。
symsx
f=x*cos（x）-2;
x0=-4;
x2=-2;
while1
x1=0.5*（x0+x2）;
f0=subs（f,x0）;
f1=subs（f,x1）;
c=f0*f1;
ifc<0
x2=x1;
else
x0=x1;
end
l=abs（x2-x0）;
ifl<=0.001
disp（['二分法的解为x=',num2str（x0）]）
break
end
end
运行后的结果为
。
（1）Lagrange插值程序：
functionyv=lagran（xd,yd,xv）
symsx
num=length（xd）;%计算节点的个数
f=0;
fork=1:
num%创建lagrange插值函数
index=setdiff（1:
num,k）;
f=f+yd（k）*prod（（x-xd（index））./（xd（k）-xd（index）））;
end
yv=subs（f,xv）;%计算xv点插值函数值
（2）绘制函数图程序：
symsx
f=1./（1+x.^2）;
xd2=-5:
2:
5;
yd2=subs（f,xd2）;
x0=-5:
0.1:
5;
y=subs（f,x0）;
plot（x0,y）
holdon
y2=lagran（xd2,yd2,x0）;
plot（x0,y2）
（3）运行程序后得到的步长分别为2，1，
的插值函数图与原图形的比较图如下：
（1）三次样条插值程序：
functionyv=csi（xd,yd,xv,h,mf,ml）
num_xd=length（xd）;
n=num_xd-1;
a_m=1:
num_xd-2;
fori=1:
num_xd-2
a_m（i）=2;
end
a_np=1:
num_xd-3;
fori=1:
num_xd-3
a_np（i）=0.5;
end
A=diag（a_m）+diag（a_np,-1）+diag（a_np,1）;
g=1:
n-1;
fork=1:
n-1
g（k）=3*（（0.5/h）*（yd（k+2）-yd（k+1））+（0.5/h）*（yd（k+1）-yd（k）））;
end
b=1:
n-1;
b
（1）=g
（1）-0.5*mf;
b（n-1）=g（n-1）-0.5*ml;
fork=2:
n-2
b（k）=g（k）;
end
b=b.';
M=A\b;
m=1:
n+1;
m
（1）=mf;
m（n+1）=ml;
fori=2:
n
m（i）=M（i-1）;
end
num_xv=length（xv）;
yv=1:
num_xv;
fori=1:
num_xv
forj=1:
num_xd
ifxv（i）k1=j-1;
k2=k1+1;
break
end
end
yv1=（h+2*（xv（i）-xd（k1）））*（xv（i）-xd（k2））^2*（yd（k1）/h^3）;
yv2=（h-2*（xv（i）-xd（k2）））*（xv（i）-xd（k1））^2*（yd（k2）/h^3）;
yv3=（xv（i）-xd（k1））*（xv（i）-xd（k2））^2*（m（k1）/h^2）;
yv4=（xv（i）-xd（k2））*（xv（i）-xd（k1））^2*（m（k2）/h^2）;
yv（i）=yv1+yv2+yv3+yv4;
end
（2）绘图程序（改变程序中的h值即可改变步长）：
clearall
clc
h=0.5;
xd=-5:
h:
5;
yd=1./（1+xd.^2）;
xv=-5:
0.1:
5;
y=1./（1+xv.^2）;
mf=0.014793;
ml=-0.014793;
yv=csi（xd,yd,xv,h,mf,ml）;
plot（xv,y）
holdon
plot（xv,yv）
（3）步长为2，1，
时原函数图与插值函数图的比较：
（1）复化梯形公式计算积分程序：
functiont=trape（f,a,b,n）
h=（b-a）/n;
index=（a+h）:
h:
（b-h）;
t1=sum（subs（f,index））;
t=（（b-a）/（2*n））*（subs（f,a）+2*t1+subs（f,b））;
（2）复化Simpson公式计算程序：
functions=simpson（f,a,b,n）
h=（b-a）/n;
index1=（a+0.5*h）:
h:
（b-0.5*h）;
index2=（a+h）:
h:
（b-h）;
s1=sum（subs（f,index1））;
s2=sum（subs
升级会员 