大连理工大学矩阵与数值分析上机作业汇编.docx
- 文档编号:23240063
- 上传时间:2023-05-15
- 格式:DOCX
- 页数:51
- 大小:307.04KB
大连理工大学矩阵与数值分析上机作业汇编.docx
《大连理工大学矩阵与数值分析上机作业汇编.docx》由会员分享,可在线阅读,更多相关《大连理工大学矩阵与数值分析上机作业汇编.docx(51页珍藏版)》请在冰豆网上搜索。
大连理工大学矩阵与数值分析上机作业汇编
矩阵与数值分析上机作业
学校:
大连理工大学
学院:
班级:
姓名:
学号:
授课老师:
注:
编程语言Matlab
1.考虑计算给定向量的范数:
输入向量力=(仏如…,切)T,输出||x||i9||x||29|k||oo«请编制一个通用程序,并用你编制的程序计算如下向童的范数:
对n=10,100,1000甚至更大的几计算其范数,你会发现什么结果?
你能否修改你的程序使得计算结果相对精确呢?
Norm,m函数
functions=Norm(x,m)
务求向量次的范数
务m取1,2,inf分别农示1,2,无穷范数n=length(x);
s=0;
switchm
case1务1-范数
fori=l:
n
s=s+abs(x(i));
case2磋-范数
fori=l:
n
S=s+x(i)A2;
end
s=sqrt(s);
caseinf%无穷-范数
s=ma:
-:
(abs(x));
end
计算向■X,y的范数
Testi.m
clearall;
clc;
nl=10;n2=100;n3=1000;
xl=1・/[l:
nl]f;x2=1・/[l:
n2]1;>:
3=1・/[l:
n3]1;
yl二Ty2=[l:
n2]Ty3=[l:
n3]T
disp(fn=10时');
dispf'x的1-范数「);disp(Norm(xl,1));dispfx的2-范数:
');disp(Norm(xl,2));
disp('x的无穷-范数:
■);disp(Norm(xlzinf));
dispfy的1-范数「);disp(Norm(ylf1));
disp{*y的2-范数:
1);disp(Norm(yl,2));
disp('y的无穷-范数:
;disp(Norm(ylzinf));disp('n=100时');
dispfx的1-范数:
');disp(Norm(x2f1));
dispfx的—范数「);disp(Norm(x2,2));disp('x的无穷-范数:
;disp(Norm(>:
2zinf));
的1-范数:
');disp(Norm(y2,1));disp(1丫的2■范数:
.);disp(Norm(y2,2));disp(1y的无穷-范数:
;disp(Norm(y2zinf));disp('11=1000时');
disp(fx的1-范数「);disp(Norm(x3,l));dispfx的2-范数「);disp(Norm(x3r2));disp('x的无穷-范数:
•);disp(Norm(>:
3,inf));disp(1y的2-范数:
');disp(Norm(y3z1));disp('y的2-范数:
1);disp(Norm(y3r2));disp('y的无方-范数:
’);disp(Norm(y3,inf));
y的1-范数:
500500;y的2-范数48272+004;y的无穷-范数:
1000
2.考虑匕=/(’)=芈也,其中定义/(0)=b此时/(巧是连续函数.用此•公式计算当[-10-叫10-15]时的函数值,画出图像。
另一方面,考虑下面算法:
〃=1+丄
if(1=1then
y=1
else
!
/=lnd/(d-l)
endif
用此算法计算xe[-10-15.10-15]时的函数值,画出图像•比较一下发生了什么?
酢
Test2.m
clearall;
clc;
n=100;%K|H]
h=2*10A(-15)/n;务步长
x=-10A(-15):
h:
10A(-15);
务第•种原函数
fl=zeros(1,n+1);
fork=l:
n+1
ifx(k)-=0
fl(k)=log(l+x(k))/x(k);
else
fl(k)=l;
end
end
subplot(2z1,1);
plot(x,flz--r*);
axis([-10A(-15)z10A(-15)z-lz2]);legend(f原图J;
$第二种算法
f2=zeros(1,n+1);
fork=l:
n+1
d=l+x(k);
if(d-=l)
f2(k)=log(d)/(d-1);
else
f2(k)=l;
end
end
subplot(2,1,2);
plot(xzf2z--r1);
axis([-10A(-15)z10A(-15)z-lz2]);legend('第二种算法1);
运行结果:
11111••d■11
第二种算法]
-
-
-
111
11111
0.8-0.6-0.4020020.40.60.8
x1严
显然第二种算法结果不准确,是因为汁算机中的舍入误差造成的,当xe[-10-15.10,s]时,d=l+x,计算机进行舍入造成d恒等于1,结果函数值恒为1。
3.首先编制一个利用秦九昶算法计算一个多项式在给定点的函数值的通用程序,你
的程序包括输入多项式的系数以及给定点,输出函数值.利用你编制的程序计算
p(x)=(x-2)9=x9-18r8+144x7-672护+2016式一4032r44-5376/-4608x2+230牡-512在工=2邻域附近的值。
画出p(t)在[1.95,20.5]上的图像.
酢:
秦九韶算法:
QinJS.m
functiony=QinJS(a,x)
务丫输出函数值
%a多项式系数,由高次到零次
显给定点、
n=length(a);
s=a(l);
fori=2:
n
s=s*x+a(i);
end
Y=s;
计算p(x):
test3.m
clearall;
clc;
・6:
0・2:
2・4;%>:
=2的邻域
dispCx=2的邻域:
f);x
a=[l-18144-6722016-40325376-46082304-512];
p=zeros(1,5);
fori=l:
5
p(i)=QinJS(a,x(i));
end
disp(•应多项式p值:
,);p
xk=l.95:
0.01:
20.5;
nk=length(xk);
pk=zeros(1,nk);
k=l;
fork=l:
nk
pk(k)=QinJS(azxk(k));
end
plot(xkzpk,*-r*);
xlabel('h');ylabel(fp(x)1);
运行结果:
x=2的邻域:
1.60001.80002.00002.20002.4000
相应多项式P值:
P=
1・0e-003*
-0.2621-0.000500.00050.2621
p(x)在•丫€[1・95,20・5]上的图像
4.编制计算给定矩阵4的厶U分解和PEU分解的通用程序,然后用你编制的程序完成下面两个计算任务:
(1)考虑
10
-1・•・
•••
•
•
•
0
■
■
■
1
■
•
•
A=
•■
••
■•
•
•
•
0
1
€Rnxn,
-1…
-1
1
1
-1・・・
-1
-1
1
自己取定并计算b=Ar•然后用你编制的不选主元和列主元的Gauss消去法求解该方程组,记你计算出的解为筑对72从5到30估计计算解的精度。
(2)对几从5到30计算其逆矩阵•
LU分解.LUDecom.m
function[L,U]=LUDecom(A)
务不带列主元的LU分解
N=size(A);
n=N(l);
L=eye(n);U=zeros(n);
fori=l:
n
U(lzi)=A(lzi);L(i,l)=A(izl)/U(lzl);
end
fori=2:
n
forj=i:
n
2=0;
fork=l:
i-1
z=z+L(izk)*U(kzj);
end
U(i,j)=A(i,j)-z;
end
forj=i+l:
n
z=0;
fork=l:
i-1
z=z+L(jzk)*U(k/i);
end
L(jzi)=(A(jzi)-z)/U(i,i);
end
end
PLU分解,PLUDecom.m
function[P,LzU]=PLUDecom(A)
务带列主元的LU分解
[mzm]=size(A);
U=A;
P=eye(m);
L=eye(m);
fori=l:
m
forj=i:
m
t(j)=U(jzi);
fork=l:
i-1
end
end
a=i;b=abs(t(i));
forj=i+l:
m
ifb b=abs(t(j)); a=j; end end if forj=1: m c=U(izj); U(izj)=U(azj); U(a,j)=c; end forj=1: m c=P(izj); P(i,j)=P(a,j); P(a,j)=c; end c=t(a); t(a)=t(i); t(i)=c; end U(i,i)=t(i); forj=i+l: m U(j,i)=t(j)/t(i); end forj=i+l: m fork=l: i-1 U(i,j)=U(izj)-U(i,k)*U(kzj); end end end L=tril(U,-1)+eye(m); U=triu(U,0); (1) (2)程序: Test4.m clearall; clc; forn=5: 30 x=zeros(nf1); A=-ones(n); A(: n)=ones(n,1); fori=l: n A(i,i)=l; forj=(i+1): (n-1) A(i,j)=0; end x(i)=l/i; end disp(1,Pin=');disp(n); disp('方程精确解: ; X b=A*x;b系数b disp「利用LU分解方程组的解: •); [LzU]=LUDecom(A);%LU分解 xLU=U\(L\b) disp「利用PLU分解方程组的解「); [P,L,U]=PLUDecom(A);%PLU分解 xPLU=U\(L\(P\b))电求解A的逆矩阵 disp('A的准确逆矩阵: ); InvA=inv(A) InvAL=zeros(n);$利用LU分解求A的逆矩阵I=eye(n); fori=l: n InvAL(: i)=U\(L\I(: i)); end disp「利用LU分解的A的逆矩阵: ; InvAL End 运行结果: (1)只列出n=5,6,7的结果 当n=5 方程精确解: X=1.0000 0.5000 0.3333 0.2500 0.2000 利用LU分解方程组的解: xLU= 1.0000 0.5000 0.3333 0.2500 0.2000 利用PLU分解方程组的解: xPLU= 1.0000 0.5000 0.3333 0.2500 0.2000 当n=6 方程精确解: X= 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 利用LU分解方程组的解: xLU= 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 利用PLU分解方程组的解: xPLU= 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 当n=7 方程精确解: X= 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 利用LU分解方程组的解: xLU= 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 利用PLU分解方程组的解: xPLU= 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 (2)只列出n=5f697时A的逆矩阵的结果当5 A的准确逆矩阵: InvA= 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0.5000 -0.2500 -0.2500 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0625 利用LU分解的A的逆矩阵: InvAL= 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0.5000 -0.2500 -0.2500 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0625 当沪6 A的准确逆矩阵: InvA= 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0313 利用LU分解的A的逆矩阵: InvAL= 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0313 当n=7 A的准确逆矩阵: InvA= 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0156 -0.0156 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0 0.5000 -0・5000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0156 0.0156 利用LU分解的A的逆矩阵: InvAL= 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0156 -0.0156 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0 0.5000 -0.2500 -0.1250 -0.0625 -0・0625 0 0 0 0.5000 -0.2500 -0.1250 -0・1250 0 0 0 0 0.5000 -0.2500 -0・2500 0 0 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0156 0.0156 5.编制计算对称正定阵的Cholesky^解的通用程序,并用你编制的程斤计算Ax=b.其中卫=(切)6Rnxn,5=订b可以由你自己取定,对“从10到20脸证程序的可至性。 丽: Choiesky分解: ChoIesky.mfunctionL=Cholesky(A)N=size(A); n=N(l); L=zeros(n); L(lzl)=sqrt(A(lzl));fori=2: n L(i/l)=A(i/l)/L(l/l); end forj=2: n sl=0; fork=l: j-1 sl=sl+L(jzk)A2; end L(j,j)=sqrt(A(jzj)-sl); fori=j+1: n s2=0; fork=l: j-1 s2=s2+L(izk)*L(jzk); end L(i,j)=(A(i,j)-s2)/L(j,j); end end 计算Ax=b;Test5.m clearall;clc; forn=10: 20 A=zeros(n,n); b=zeros(n,1); fori=l: n forj=1: n A(izj)=1/(i+j-1); end b(i,l)=i; end disp(rn=r);disp(n); disp(f方程组原始解1);xO=A\b disp(,利用Cholesky分解的方程组的解*); L=Cholesky(A) x=Lr\(L\b) end 运行结果: 只列岀了11的结果 n=10 方程组原始解 x0= 1・0e+008* -0.0000 0.0010 -0.0233 0.2330 -1.2108 3.5947 -6.3233 6.5114 -3.6233 0.8407 利用Choiesky分解的方程组的解x= 1・0e+008* -0.0000 0.0010 -0.0233 0.2330 -1.2105 3.5939 -6.3219 6.5100 -3.6225 0.8405 n= 11 方程组原始解 x0= 1.0e+009* 0.0000 -0.0002 0.0046 -0.0567 0.3687 -1.4039 3.2863 -4.7869 4.2260 -2.0685 0.4305 利用Choiesky分解的方程组的解 X= 1.0e+009* 0.0000 -0.0002 0.0046 -0.0563 0.3668 -1・3972 3.2716 -4.7669 4.2094 -2.0608 0.4290 6.⑴编制程序HoxLse{x)->其作用是对输入的向童丄,输出单位向童"使得(/—21ZU7)x=11创2妝 ⑵编制Householder变换阵H=J一2nuT€庫沁几乘以qG庫心尬的程序乩4,注意,你的程序并不显式的计算出H. (3)考虑矩阵 1 -1 2 3 3辺 4、 辺 A= -2 2 e 7T -\/10 2 -3 7 0 2 7 5/2丿 用你编制的程序计算H使得的第一列为a®的形式,并将H.4的结果显示。 (1)House,m functionu=House(x)n=length(x); el=eye(n,1); w=x-norm(xz2)*el; u=w/norm(w,2); (2)Hou_A.m functionHA=Hou__A(A) al=A(: 1); n=length(al); el=eye(n,1);w=al-norm(al,2)*el; u=w/norm(w,2); H=eye(n)-2*u*u1HA=H*A; ⑶test6.m clearall; clc; A=[l234; -12sqrt (2)sqrt(3); -22exp(l)pi; -sqrt(10)2-37; 0275/2]; HA=Hou_A(A) 运行结果: H= 0.2500 -0.2500 -0.5000 -0.7906 0 -0.2500 0.9167 -0.1667 -0.2635 0 -0.5000 -0.1667 0.6667 -0.5270 0 -0.7906 -0.2635 -0.5270 0.1667 0 0 0 0 0 1.0000 HA= 4.0000 -2.5811 1・4090 -6.5378 0.0000 0.4730 0.8839 -1・7805 0.0000 -1.0541 1.6576 -3.8836 0.0000 -2・8289 -4.6770 -4.1078 0 2.0000 7.0000 2.5000 7.用Jaa)b承Gauss・Sgidel送代求•解下面的方程组,输出迭代每一步的误差||%-斗||: (5ri—X2—3物=—2 一巧+2衍十仏3=1 —3^1+4x2+15^3=10 Jacobi迭代: Jaccobi・m function[x,n]=Jaccobi(Azb,xO) 务--・方程组系数阵离 $・方程组右端顶b %
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 大连理工大学 矩阵 数值 分析 上机 作业 汇编