大连理工大学 矩阵大作业Word文档格式.docx
《大连理工大学 矩阵大作业Word文档格式.docx》由会员分享，可在线阅读，更多相关《大连理工大学 矩阵大作业Word文档格式.docx（32页珍藏版）》请在冰豆网上搜索。
=i
p=p*x+a（j）;
j=j+1;
y=p;
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;
x（j）=s/a（j,j）;
j=j-1;
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,:
x=a;
y=b;
运行的结果为：
2 列主元消去法计算程序和结果：
[a,b]=lzy（A,B）;
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）;
ab（:
num_j+1）=b（:
1）;
[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<
forn=j:
num_j+1
ab（m,n）=ab（m,n）-（ab（m,j）/max）*ab（j,n）;
m=m+1;
a1=zeros（num_i,num_j）;
num_i
a1（:
k）=ab（:
b1=ab（:
num_i+1）;
function[b,c]=searmax（i,j,a）
k=i;
m=abs（a（k,j））;
c=k;
whilek<
ifm>
=abs（a（k,j））
continue
else
m=abs（a（k,j））;
c=k;
b=a（c,j）;
（1）LU分解的计算程序及结果：
function[l,u]=lufz（a）
l=eye（num_i）;
u=eye（num_i）;
forj=1:
u（1,j）=a（1,j）;
l（j,1）=a（j,1）/u（1,1）;
i=2;
j=i;
fork=1:
i-1
s=s+l（i,k）*u（k,j）;
u（i,j）=a（i,j）-s;
s=s+l（j+1,k）*u（k,i）;
l（j+1,i）=（a（j+1,i）-s）/u（i,i）;
s=s+l（i,k）*u（k,num_i）;
u（i,num_i）=a（i,num_i）-s;
输入要求解的矩阵后运行的结果为：
（2）带列主元的LU分解计算程序及结果
由于Matlab中自带带列主元的LU分解函数，故计算程序如下：
[L,U,P]=lu（A）;
运行结果如下：
为单位阵。
（3）求A的逆矩阵的计算程序及结果：
[num_i,num_j]=size（A）;
A_n=eye（num_i）;
E=eye（num_i）;
[L,U]=lufz（A）;
fori=1:
y=hd（1,L,E（:
i））;
A_n（:
i）=hd（0,U,y）;
disp（A_n）
functionx=hd（f,a,b）
x=zeros（num_i,1）;
switchf
case0
i=num_i-1;
x（num_i）=b（num_i）/a（num_i,num_j）;
whilei>
fork=i+1:
s=s+a（i,k）*x（k）;
x（i）=（b（i）-s）/a（i,i）;
x
（1）=b
（1）/a（1,1）;
j=2;
j-1
s=s+a（j,k）*x（k）;
x（j）=（b（j）-s）/a（j,j）;
error！
请输入正确的参数！
运行结果：
（4）求A的行列式的计算程序及结果：
s=1;
s=s*U（i,i）;
disp（['
矩阵的行列式值为：
num2str（s）]）
运行程序后结果与调用matlab中det（）函数结果如下：
（1）求cholesky分解程序及结果：
functionl=chole（a）
l=zeros（num_i）;
l（1,1）=sqrt（a（1,1））;
fork=2:
l（k,1）=a（k,1）/l（1,1）;
j=2;
=num_j
s1=0;
s1=s1+l（j,k）^2;
l（j,j）=sqrt（a（j,j）-s1）;
s2=0;
s2=s2+l（i,k）*l（j,k）;
l（i,j）=（a（i,j）-s2）/l（j,j）;
运行程序后的结果为：
（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）;
运行后的结果为:
程序和结果如下：
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（'
方程组的解为：
function[l,u]=sjfj（a）
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）;
l=tril（a）;
u=triu（a）;
num_i;
l（j,j）=1;
运行程序后结果为：
编程求解矩阵A的QR分解：
（1）QR分解计算程序：
function[Q,R]=hqr（A）
[n,n]=size（A）;
Q=eye（n）;
（n-1）
E=eye（n-i+1）;
e1=E（:
a=zeros（1,n-i+1）'
;
forj=1:
（n-i+1）
a（j）=A（j+i-1,i）;
av=sqrt（a'
*a）;
w=a-av*e1;
h=E-（2/（w'
*w））*（w*w'
q=eye（n）;
fork=i:
n
forj=i:
q（k,j）=h（k-i+1,j-i+1）;
A=q*A;
Q=q*Q;
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
a=a0;
b=b0;
c=c0;
方程组的解为x'
运行后的结果为：
（2）Gauss-Seidel迭代法求解的程序与结果：
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）;
x=[a,b,c];
（1）计算
在
上的根的程序：
symsx
f=2*x^3-5*x-1;
df=diff（f,x）;
g=x-（f/df）;
x0=1;
x1=subs（g,x0）;
dx=abs（x1-x0）;
ifdx<
disp（['
Newton法的解为：
x='
num2str（x1）]）
x0=x1;
x1=1.1;
x2=x1-（subs（f,x1）/（subs（f,x1）-subs（f,x0）））*（x1-x0）;
dx=abs（x2-x1）;
割线法的解为：
num2str（x2）]）
x1=x2;
运行后结果为：
Newton法
，割线法
（2）计算
f=exp（x）*sin（x）;
x0=-2;
x1=-2.1;
f=x*cos（x）-2;
x0=-4;
x2=-2;
x1=0.5*（x0+x2）;
f0=subs（f,x0）;
f1=subs（f,x1）;
c=f0*f1;
ifc<
x2=x1;
l=abs（x2-x0）;
ifl<
二分法的解为x='
num2str（x0）]）
运行后的结果为
（1）Lagrange插值程序：
functionyv=lagran（xd,yd,xv）
num=length（xd）;
%计算节点的个数
f=0;
num%创建lagrange插值函数
index=setdiff（1:
num,k）;
f=f+yd（k）*prod（（x-xd（index））./（xd（k）-xd（index）））;
yv=subs（f,xv）;
%计算xv点插值函数值
（2）绘制函数图程序：
f=1./（1+x.^2）;
xd2=-5:
2:
5;
yd2=subs（f,xd2）;
x0=-5:
0.1:
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;
num_xd-2
a_m（i）=2;
a_np=1:
num_xd-3;
num_xd-3
a_np（i）=0.5;
A=diag（a_m）+diag（a_np,-1）+diag（a_np,1）;
g=1:
n-1;
n-1
g（k）=3*（（0.5/h）*（yd（k+2）-yd（k+1））+（0.5/h）*（yd（k+1）-yd（k）））;
b=1:
b
（1）=g
（1）-0.5*mf;
b（n-1）=g（n-1）-0.5*ml;
n-2
b（k）=g（k）;
b=b.'
M=A\b;
m=1:
n+1;
m
（1）=mf;
m（n+1）=ml;
fori=2:
m（i）=M（i-1）;
num_xv=length（xv）;
yv=1:
num_xv;
num_xv
num_xd
ifxv（i）<
xd（j）
k1=j-1;
k2=k1+1;
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;
（2）绘图程序（改变程序中的h值即可改变步长）：
clearall
clc
h=0.5;
xd=-5:
h:
yd=1./（1+xd.^2）;
xv=-5:
y=1./（1+xv.^2）;
mf=0.014793;
ml=-0.014793;
yv=csi（xd,yd,xv,h,mf,ml）;
plot（xv,y）
plot（xv,yv）
（3）步长为2，1，
时原函数图与插值函数图的比较：
（1）复化梯形公式计算积分程序：
functiont=trape（f,a,b,n）
h=（b-a）/n;
index=（a+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）
index1=（a+0.5*h）:
（b-0.5*h）;
index2=（a+h）:
s1=sum（subs（f,index1））;
s2=sum（subs（
升级会员 