大连理工大学矩阵与数值分析上机作业.docx
- 文档编号:15510700
- 上传时间:2023-07-05
- 格式:DOCX
- 页数:51
- 大小:85.27KB
大连理工大学矩阵与数值分析上机作业.docx
《大连理工大学矩阵与数值分析上机作业.docx》由会员分享,可在线阅读,更多相关《大连理工大学矩阵与数值分析上机作业.docx(51页珍藏版)》请在冰点文库上搜索。
大连理工大学矩阵与数值分析上机作业
矩阵与数值分析上机作业
学校:
大连理工大学
学院:
班级:
姓名:
学号:
授课老师:
注:
编程语言Matlab
1.琴虑计算给定働量的葩址输入向量広』(巾斑…宀产输出||工||“||工|怙㈣心请编制一牛通用程序,并用你編制的程序计算如下询量的范数:
对网1加,wm甚至更大的“计算其范数,你会发现什幺结粟?
你能否修改你的程序使得计算绪果相时赫■确呢?
程序:
Norm.m函数
functions=Norm(x,m)
%求向量x的范数
%mx1,2,inf分别表示1,2,无穷范数
n=length(x);
s=0;
switchm
case1%1-范数
fori=1:
n
s=s+abs(x(i));
end
case2%2-范数
fori=1:
n
s=s+x(if2;
end
s=sqrt(s);
caseinf%无穷-范数
s=max(abs(x));
end
计算向量x,y的范数
Test1.m
clearall;
clc;
n1=10;n2=100;n3=1000;
x1=1./[1:
n1]';x2=1./[1:
n2]';x3=1./[1:
n3]';
y1=[1:
n1]';y2=[1:
n2]';y3=[1:
n3]';
disp('n=10时');
disp('x的1-范数:
');disp(Norm(x1,1));
disp('x的2-范数:
');disp(Norm(x1,2));
disp('x的无穷-范数:
');disp(Norm(x1,inf));
disp('y的1-范数:
');disp(Norm(y1,1));
disp('y的2-范数:
');disp(Norm(y1,2));
disp('y的无穷-范数:
');disp(Norm(y1,inf));
disp('n=100时');
disp('x的1-范数:
');disp(Norm(x2,1));
disp('x的2-范数:
');disp(Norm(x2,2));
disp('x的无穷-范数:
');disp(Norm(x2,inf));
disp('y的1-范数:
');disp(Norm(y2,1));
disp('y的2-范数:
');disp(Norm(y2,2));
disp('y的无穷-范数:
');disp(Norm(y2,inf));
disp('n=1000时');
disp('x的1-范数:
');disp(Norm(x3,1));
disp('x的2-范数:
');disp(Norm(x3,2));
disp('x的无穷-范数:
');disp(Norm(x3,inf));
disp('y的1-范数:
');disp(Norm(y3,1));
disp('y的2-范数:
');disp(Norm(y3,2));
disp('y的无穷-范数:
');disp(Norm(y3,inf));
运行结果:
n=10时
范数:
1
-范数:
10
-范数:
1
x的1-范数29290;x的2-范数:
1.2449;x的无穷-
y的1-范数:
55;y的2-范数:
19.6214;y的无穷n=100时
x的1-范数:
5.1874;x的2-范数:
1.2787;x的无穷
y的1-范数:
5050;
的2-范数:
581.6786;y的无穷-范数:
100
n=1000时
x的1-范数74855;x的2-范数:
1.2822;x的无穷-范数:
1
y的1-范数:
500500;y的2-范数:
1.8271e+004;y的无穷-范数:
1000
2.耆虑砂==呼^其中定51/(0)=此时几期是连绽函戟.用此公式计算当工“―1旷巾U)-缪时的函数值*風出图像.另一方面*哮虑下面算法:
d1+j
1/(/=1tbfjj
1/=1
仙
y=liid/(d—1(
endif
用此算法计%€[-10-0io_is]时的圉数血画出图像.比校一下岌生了什么?
程序
Test2.m
clearall;
clc;
n=100;%区间
h=2*10A(-15)/n;%步长
x=-10A(-15):
h:
10A(-15);
%第一种原函数
f1=zeros(1,n+1);
fork=1:
n+1
ifx(k)~=0
f1(k)=log(1+x(k))/x(k);
else
f1(k)=1;
endend
subplot(2,1,1);
plot(x,f1,'-r');
axis([-10A(-15),10A(-15),-1,2]);
legend('原图');
%第二种算法
f2=zeros(1,n+1);
fork=1:
n+1
d=1+x(k);
if(d~=1)f2(k)=log(d)/(d-1);
else
f2(k)=1;
end
end
subplot(2,1,2);
plot(x,f2,'-r');
axis([-10A(-15),10A(-15),-1,2]);
legend('第二种算法');
运行结果:
農IQ
显然第二种算法结果不准确,是因为计算机中的舍入误差造成的,当X[1015,1015]时,
d1x,计算机进行舍入造成d恒等于1,结果函数值恒为1。
3.首先编制一个別用秦丸昶律法计算一个多项式在给定点的函數值的遇用程序,你
的程序邑括揄入多项式的系就以及给定点,嶽出函数值.利用你騙制的粒序计算
p(x)=(^-2)s=泸一18rs+IMj-7-672^+2nitkrs-lO32r4+越芮护-460SX2+230b-512在工=2邻域附近的值-画出此r)在jt€|1-航20-司上的图像・
程序:
秦九韶算法:
QinJS.m
functiony=QinJS(a,x)
%y输岀函数值
%a多项式系数,由高次到零次
%x给定点
n=length(a);
s=a
(1);
fori=2:
n
s=s*x+a(i);
end
y=s;
计算p(x):
test3.m
clearall;
clc;
x=1.6:
0.2:
2.4;%x=2的邻域
disp('x=2的邻域:
');x
a=[1-18144-6722016-40325376-46082304-512];
p=zeros(1,5);
fori=1:
5
p(i)=QinJS(a,x(i));
end
disp('相应多项式p值:
');p
xk=1.95:
0.01:
20.5;
nk=length(xk);
pk=zeros(1,nk);
k=1;
fork=1:
nk
pk(k)=QinJS(a,xk(k));
end
plot(xk,pk,'-r');
xlabel('x');ylabel('p(x)');
运行结果:
x=2的邻域:
x=
1.60001.80002.00002.20002.4000
相应多项式p值:
P=
1.0e-003*
-0.2621-0.000500.00050.2621
p(x)在X[1.95,20.5]上的图像
4,鎬制计雾给定拒阵冲的U7分解和FL17分解的通用程序、然后用你编制的程序完咸下面两个计鼻任务:
(1)考虑
II
(I
-1
自己取定上eRS并计鼻古=屉.然后用你编制的不选主元和列主无的伽u抽消去法求解谏方程组,记你计算韭的解为氏对和从5刘盹估计计界解的精度*
⑵对仕从訐到30计算其逆雑阵.
程序:
LU分解,LUDe..m
function[L,U]=LUDe.(A)
%不带列主元的LU分解
N=size(A);
n=N
(1);
L=eye(n);U=zeros(n);
fori=1:
n
U(1,i)=A(1,i);L(i,1)=A(i,1)/U(1,1);
end
fori=2:
n
forj=i:
n
z=0;
fork=1:
i-1
z=z+L(i,k)*U(k,j);
U(i,j)=A(i,j)-z;
end
forj=i+1:
n
z=0;
fork=1:
i-1z=z+L(j,k)*U(k,i);
end
L(j,i)=(A(j,i)-z)/U(i,i);
end
end
PLU分解,PLUDe..m
function[P,L,U]=PLUDe.(A)%带列主元的LU分解[m,m]=size(A);
U=A;
P=eye(m);
L=eye(m);
fori=1:
m
forj=i:
m
t(j)=U(j,i);
fork=1:
i-1
t(j)=t(j)-U(j,k)*U(k,i);
end
a=i;b=abs(t(i));
forj=i+1:
m
ifb end end ifa~=i forj=1: mc=U(i,j);U(i,j)=U(a,j);U(a,j)=c; end forj=1: mc=P(i,j);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+1: m U(j,i)=t(j)/t(i); forj=i+1: m fork=1: i-1 U(i,j)=U(i,j)-U(i,k)*U(k,j); end end end L=tril(U,-1)+eye(m); U=triu(U,0); (1) (2)程序: Test4.m clearall; clc; forn=5: 30 x=zeros(n,1); A=-ones(n); A(: n)=ones(n,1); fori=1: n A(i,i)=1; forj=(i+1): (n-1) A(i,j)=0; end x(i)=1/i; end disp('当n=');disp(n); disp('方程精确解: '); x b=A*x;%系数b disp('利用LU分解方程组的解: '); [L,U]=LUDe.(A);%LU分解 xLU=U\(L\b) disp('利用PLU分解方程组的解: '); [P,L,U]=PLUDe.(A);%PLU分解 xPLU=U\(L\(P\b)) %求解A的逆矩阵 disp('A的准确逆矩阵: '); InvA=inv(A) InvAL=zeros(n);%利用LU分解求A的逆矩阵 I=eye(n); fori=1: 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.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.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=5,6,7时A的逆矩阵的结果 当n=5 A的准确逆矩阵: InvA= 0.5000-0.2500-0.1250-0.0625-0.0625 00.5000-0.2500-0.1250-0.1250 000.5000-0.2500-0.2500 0000.5000-0.5000 0.50000.25000.12500.06250.0625 利用LU分解的A的逆矩阵: InvAL= 0.5000-0.2500-0.1250-0.0625-0.0625 00.5000-0.2500-0.1250-0.1250 000.5000-0.2500-0.2500 0000.5000-0.5000 0.50000.25000.12500.06250.0625 当n=6 A的准确逆矩阵: InvA= 0.5000-0.2500-0.1250-0.0625-0.0313-0.0313 00.5000-0.2500-0.1250-0.0625-0.0625 000.5000-0.2500-0.1250-0.1250 0000.5000-0.2500-0.2500 00000.5000-0.5000 0.50000.25000.12500.06250.03130.0313 利用LU分解的A的逆矩阵: InvAL= 0.5000-0.2500-0.1250-0.0625-0.0313-0.0313 00.5000-0.2500-0.1250-0.0625-0.0625 000.5000-0.2500-0.1250-0.1250 0000.5000-0.2500-0.2500 00000.5000-0.5000 0.50000.25000.12500.06250.03130.0313 当n=7 A的准确逆矩阵: InvA= 0.5000-0.2500-0.1250-0.0625-0.0313-0.0156-0.0156 00.5000-0.2500-0.1250-0.0625-0.0313-0.0313 000.5000-0.2500-0.1250-0.0625-0.0625 0000.5000-0.2500-0.1250-0.1250 00000.5000-0.2500-0.2500 000000.5000-0.5000 0.50000.25000.12500.06250.03130.01560.0156 利用LU分解的A的逆矩阵: InvAL= 0.5000-0.2500-0.1250-0.0625-0.0313-0.0156-0.0156 00.5000-0.2500-0.1250-0.0625-0.0313-0.0313 000.5000-0.2500-0.1250-0.0625-0.0625 0000.5000-0.2500-0.1250-0.1250 00.5000-0.2500-0.2500 00.5000-0.5000 0.50000.25000.12500.06250.03130.01560.0156 5.蜻制计算对称正定薛的OioteM洽解的遇用程序,幷用你编制的程序计算如三乩其中丸=(aij)€即".旳=b可以由你自己取忌对必从10到凶验证程序的可命出 程序: Cholesky分解: Cholesky.m functionL=Cholesky(A) N=size(A); n=N (1); L=zeros(n); L(1,1)=sqrt(A(1,1)); fori=2: n L(i,1)=A(i,1)/L(1,1); end forj=2: n s1=0; fork=1: j-1 s1=s1+L(j,k)A2; end L(j,j)=sqrt(A(j,j)-s1); fori=j+1: n 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); end end 计算Ax=b;Test5.m clearall;clc; forn=10: 20 A=zeros(n,n); b=zeros(n,1); fori=1: n forj=1: n A(i,j)=1/(i+j-1);end b(i,1)=i; end disp( 'n=');disp(n); disp( '方程组原始解');x0=A\b disp( '利用Cholesky分解的方程组的解'); L=Cholesky(A) x=L'\(L\b) end 运行结果: 只列出了n=10,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 利用Cholesky分解的方程组的解 x= 1.0e+008* -0.0000 0.0010 -0.0233 0.2330 -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 利用Cholesky分解的方程组的解 x= 1.0e+009* 0.0046 -0.0563 0.3668 -1.3972 3.2716 -4.7669 4.2094 -2.0608 0.4290 6. (1)编制强序其作用是时辅: 入的向董「输出单位向量“彳吏得(mN)』= ⑵编制换薛丹= 序并怎显式的计算出 (3)考虑距阵 汀—加十丘朗旳乘EMW肿山的程序厅扎注意、你的程 1234、 -1乂暑亞 A=-22ff, -7102-37 k0275/2; 用你编制的程序计算H使得的第一列为“I的形式,并将乩4的结果显示. 程序: (1)House.m functionu=House(x) n=length(x); e1=eye(n,1); w=x-norm(x,2)*e1; u=w/norm(w,2); functionHA=Hou_A(A) a1=A(: 1); n=length(a1); e1=eye(n,1);w=a1-norm(a1,2)*e1; u=w/norm(w,2); H=eye(n)-2*u*u' HA=H*A; (3)test6.m clearall; clc; A=[1234; -12sqrt (2)sqrt(3); -22exp (1)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 00 1.0000 HA= 0.00000.47300.8839-1.7805 0.0000-1.05411.6576-3.8836 0.0000-2.8289-4.6770-4.1078 02.00007.00002.5000 7.用代求解下面的方程组*输出迭代每一步的误差||斗—斗||: 5xj—X2—扫=—2 «—叼十2^2十4^3=1 —3g+4电4-15^3=10 程序: Jacobi迭代: Jaccobi.m function[x,n]=Jaccobi(A,b,x0)%--•方程组系数阵A %--•方程组右端顶b %--初始值x0 %--求解要求精确度eps %--迭代步数控制M %--•返回求得的解x %--•返回迭代步数n M=1000;eps=1.0e-5; D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); %求A的对角矩阵 %求A的下三角阵 %求A的上三角阵 J=D\(L+U); f=D\b; x=J*x0+fn=1;%迭代次数 err=norm(x-x0,inf) while(err>=eps) x0=x; x=J*x0+f n=n+1; err=norm(x-x0,inf) if(n>=M) ? '); disp('Warning: 迭代次数太多,可能不收敛return; end end Gauss_Seidel迭代: Gauss_Seidel.mfunction[x,n]=Gauss_Seidel(A,b,x0)%--Gauss-Seidel迭代法解线性方程组%--方程组系数阵A %--方程组右端项b %--初始值x0 %--求解要求的精确度eps %--迭代步数控制M %--返回求得的解x%--返回迭代步数n eps=1.0e-5; M=10000; D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1);
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 大连理工大学 矩阵 数值 分析 上机 作业