大连理工大学矩阵上机.docx
- 文档编号:3682842
- 上传时间:2022-11-24
- 格式:DOCX
- 页数:47
- 大小:153.99KB
大连理工大学矩阵上机.docx
《大连理工大学矩阵上机.docx》由会员分享,可在线阅读,更多相关《大连理工大学矩阵上机.docx(47页珍藏版)》请在冰豆网上搜索。
大连理工大学矩阵上机
大连理工大学
DALIANUNIVERSITYOFTECHNOLOGY
研究生作业
课程名称:
矩阵与数值分析
研究生姓名:
学号:
作业成绩:
任课教师签名
能源与动力学院
上交作业时间:
2015年12月25日
1.设
,其精确值为
.
(1)编制按从大到小顺序
,计算
的通用程序;
(2)编制按从小到大顺序
,计算
的通用程序;
(3)按两种顺序分别计算
,并指出有效位数(编制时的单精度);
(4)通过本次上机题,你明白了什么?
MATLAB程序
function[s1,s2]=s(N)
formatlongg
s1=0;s2=0;
forj=2:
N
s1=1.0/(j*j-1)+s1;
end
forj=N:
-1:
2
s2=1.0/(j*j-1)+s2;
end
end
计算结果
>>[s1,s2]=s(1.0e2)
s1=
0.740049504950495
s2=
0.740049504950495
>>[s1,s2]=s(1.0e4)
s1=
0.749900004999506
s2=
0.7499000049995
>>[s1,s2]=s(1.0e6)
s1=
0.749999000000522
s2=
0.7499990000005
结果分析
计算机做加减法时,运算次序会影响结果,计算和时应先安排绝对值小的数参与运算,这样能取得较高的精度。
2.秦九韶算法。
已知n次多项式
,用秦九韶算法编写通用的程序计算函数在
点的值,并计算
在23点的值。
(提示:
编写程序时,输入系数向量和点
,输出结果,多项式的次数可以通过向量的长度来判断).
MATLAB程序
A=input('ÇëÊäÈëϵÊý,ÓɸߴÎÃÝ¿ªÊ¼');
n=input('ÇëÊäÈë¼ÆËã±äÁ¿µÄÖµ');
len=length(A);
val=zeros(len);
val
(1)=A
(1);
%%printf('len=%c',len)
fori=2:
1:
len
%disp(val(i-1))
%disp(n)
val(i)=val(i-1)*n+val(i);
end
%printf('?
?
?
?
?
¨¢?
?
?
?
%f',val(len))
val(len)
计算结果
请输入系数,由高次幂开始[73-511]
请输入计算变量的值23
ans=
85169
3.分别用Gauss消元法和列主元消去法编程求解方程组Ax=b,其中
。
Gauss消元法
functionx=gauss(a,b)
n=length(b);
fork=1:
n-1
ifa(k,k)==0
fprintf('Error:
thepivotelementequaltozero!
\n',k);
return;
end
index=[k+1:
n];
m=-a(index,k)/a(k,k);
a(index,index)=a(index,index)+m*a(k,index);
b(index)=b(index)+m*b(k);
end
x=zeros(n,1);
x(n)=b(n)/a(n,n);
fori=n-1:
-1:
1
x(i)=(b(i)-a(i,[i+1:
n])*x([i+1:
n]))/a(i,i);
end
计算结果
>>a=[31,-13,0,0,0,-10,0,0,0;-13,35,-9,0,-11,0,0,0,0;0,-9,31,-10,0,0,0,0,0;0,0,-10,79,-30,0,0,0,-9;0,0,0,-30,57,-7,0,-5,0;0,0,0,0,-7,47,-30,0,0;0,0,0,0,0,-30,41,0,0;0,0,0,0,-5,0,0,27,-2;0,0,0,-9,0,0,0,-2,29]
>>b=[-15;27;-23;0;-20;12;-7;7;10]
x=gau(a,b)
x=
-0.289233816015754
0.345435715779115
-0.712811731086879
-0.220608510570529
-0.430400432704022
0.154308739838311
.0578********
0.201053894823681
0.290228661879745
列主元消去法
function[x]=gauss(a,b)
n=length(a);
x=zeros(n,1);
a=[ab];
fork=1:
n-1
max=k;
fori=k+1:
n
ifa(i,k)>a(max,k)
max=i;
end
end
temp=a(k,k:
n+1);
a(k,k:
n+1)=a(max,k:
n+1);
a(max,k:
n+1)=temp;
fori=k+1:
n
a(i,k)=-a(i,k)/a(k,k);
a(i,k+1:
n+1)=a(i,k+1:
n+1)+a(i,k)*a(k,k+1:
n+1);
end
end
x(n,1)=a(n,n+1)/a(n,n);
fori=n-1:
-1:
1
sum=0;
forj=i+1:
n
sum=sum+x(j,1)*a(i,j);
end
x(i,1)=(a(i,n+1)-sum)/a(i,i);
end
计算结果
>>a=[31,-13,0,0,0,-10,0,0,0;-13,35,-9,0,-11,0,0,0,0;0,-9,31,-10,0,0,0,0,0;0,0,-10,79,-30,0,0,0,-9;0,0,0,-30,57,-7,0,-5,0;0,0,0,0,-7,47,-30,0,0;0,0,0,0,0,-30,41,0,0;0,0,0,0,-5,0,0,27,-2;0,0,0,-9,0,0,0,-2,29]
>>b=[-15;27;-23;0;-20;12;-7;7;10]
>>[x]=gauss(a,b)
x=
-0.289233816015754
0.345435715779115
-0.712811731086879
-0.220608510570529
-0.430400432704022
0.154308739838311
.0578********
0.201053894823681
0.290228661879745
4.编程求解题3中矩阵A的LU分解及列主元的LU分解(求出L,U和P),并用LU分解的方法求A的逆矩阵及A的行列式.
LU分解
function[L,U]=mylu(a)
[n,n]=size(a);
fori=1:
n
L(i,i)=1;
end
U(1,1)=a(1,1)/L(1,1);
ifL(1,1)*U(1,1)==0
fprintf('Factorizationimpossible')
else
forj=2:
n
U(1,j)=a(1,j);
L(j,1)=a(j,1)/U(1,1);
end
end
fori=2:
n-1
sum=0;
fork=1:
i-1
sum=sum+L(i,k)*U(k,i);
end
U(i,i)=(a(i,i)-sum)/L(i,i);
ifL(i,i)*U(i,i)==0
fprintf('Factorizationimpossible')
else
forj=i+1:
n
h=0;
s=0;
fork=1:
i-1
h=h+L(i,k)*U(k,j);
s=s+L(j,k)*U(k,i);
end
U(i,j)=1/L(i,i)*(a(i,j)-h);
L(j,i)=1/U(i,i)*(a(j,i)-s);
end
end
end
sum=0;
fork=1:
n-1
sum=sum+L(n,k)*U(k,n);
U(n,n)=(a(n,n)-sum)/L(n,n);
end
end
计算结果
>>[L,U]=mylu(a)
L=
Columns1through3
100
-0.41935483870967710
0-0.3045851528384281
00-0.353872899362565
000
000
000
000
000
Columns4through6
000
000
000
100
-0.39755492585681310
0-0.1569436359494821
00-0.653976718762685
0-0.112102597106773-0.0175453757582425
-0.119266477757044-0.0833908821051642-0.0142268147798013
Columns7through9
000
000
000
000
000
000
100
-0.024618525643364810
-0.0199621375629645-0.09231647025909681
U=
Columns1through3
31-130
029.5483870967742-9
0028.2587336244541
000
000
000
000
000
000
Columns4through6
00-10
0-11-4.19354838709677
-10-3.35043668122271-1.27729257641921
75.4612710063743-31.185********5-0.451999227351748
044.6019996774714-7.17969451931716
0045.8731926371318
000
000
000
Columns7through9
000
000
000
00-9
0-5-3.57799433271131
-30-0.784718179747412-0.561543439982356
21.3806984371195-0.513187420344639-0.367236336322372
026.4130849214705-2.41999576495225
0027.3895043327769
LU分解计算a的行列式和逆矩阵
>>det(U)
ans=
6.1817e+13
>>inv(U)*inv(L)
ans=
0.03900.01600.00580.00360.00730.01760.01290.00140.0012
0.01590.03780.01310.00650.01210.00970.00710.00240.0022
0.00490.01160.03810.00770.00690.00390.00280.00150.0025
0.00080.00200.00640.01810.01050.00330.00240.00240.0058
0.00050.00110.00360.01010.02430.00700.00510.00480.0035
0.00010.00030.00100.00280.00680.04190.03060.00130.0010
0.00010.00020.00070.00210.00500.03060.04680.00100.0007
0.00010.00020.00080.00230.00480.00140.00100.03820.0033
0.00030.00060.00200.00580.00360.00110.00080.00340.0365
列主元LU分解
function[l,u,p]=lielu(a)
[m,n]=size(a);
ifm~=n
error('¾ØÕó²»ÊÇ·½Õó')
return
end
ifdet(a)==0
error('¾ØÕó²»Äܱ»Èý½Ç·Ö½â')
end
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
end
a=i;b=abs(t(i));
forj=i+1:
m
ifb b=abs(t(j)); a=j; end end ifa~=i forj=1: m c=u(i,j); u(i,j)=u(a,j); u(a,j)=c; end forj=1: m c=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); end 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); end 计算结果 >>a=[31,-13,0,0,0,-10,0,0,0;-13,35,-9,0,-11,0,0,0,0;0,-9,31,-10,0,0,0,0,0;0,0,-10,79,-30,0,0,0,-9;0,0,0,-30,57,-7,0,-5,0;0,0,0,0,-7,47,-30,0,0;0,0,0,0,0,-30,41,0,0;0,0,0,0,-5,0,0,27,-2;0,0,0,-9,0,0,0,-2,29]; >>[l,u,p]=lielu(a) l= 1.000000000000 -0.41941.00000000000 0-0.30461.0000000000 00-0.35391.000000000 000-0.39761.00000000 0000-0.15691.0000000 00000-0.65401.000000 0000-0.1121-0.0175-0.02461.00000 000-0.1193-0.0834-0.0142-0.0200-0.09231.0000 u= 31.0000-13.0000000-10.0000000 029.5484-9.00000-11.0000-4.1935000 0028.2587-10.0000-3.3504-1.2773000 00075.4613-31.1856-0.452000-9.0000 000044.6020-7.17970-5.0000-3.5780 0000045.8732-30.0000-0.7847-0.5615 00000021.3807-0.5132-0.3672 000000026.4131-2.4200 0000000027.3895 p= 100000000 010000000 001000000 000100000 000010000 000001000 000000100 000000010 000000001 5.编制程序求解矩阵A的cholesky分解,并用程序求解方程组Ax=b,其中 , . MATLAB程序 A=[21-51 1-507 021-1 16-1-4 ]; B=[13,-9,6,0]; %disp(A) dim=size(A,1); L=zeros(dim,dim); y=[] x=zeros(dim,1); fori=1: 1: dim forj=1: 1: dim ifi>=(j+1) L(i,j)=(A(i,j)-sum(L(i,1: (j-1)).*L(j,1: (j-1))))/L(j,j); else L(j,j)=sqrt(A(j,j)-sum(L(j,1: j-1).*L(j,1: j-1))); end end end L_T=L.'; fori=1: 1: dim ifi==1 y (1)=B (1)/L(1,1); else y(i)=(B(i)-sum(L(i,1: (i-1)).*y(1: (i-1))))/L(i,i); end end fori=dim: -1: 1 ifi==dim x(i)=y(i)/L(i,i); else x(i)=(y(i)-sum(L((i+1): dim,i).*x((i+1): dim)))/L(i,i); end %disp('((((((((') %disp(x(i)) end disp('x') disp(x) 计算结果 x 52.2500 -38.7500 30.7500 -52.7500 6.用追赶法编制程序求解方程组Ax=b,其中 , . MATLAB程序 A=[4,2,0,0;3,-2,1,0;0,2,5,3;0,0,-1,6]; a=[032-1];c=[213];b=[4-256];d=[62105]; n=length(b); u0=0;y0=0;a (1)=0; L (1)=b (1)-a (1)*u0; y (1)=(d (1)-y0*a (1))/L (1); u (1)=c (1)/L (1); fori=2: (n-1) L(i)=b(i)-a(i)*u(i-1); y(i)=(d(i)-y(i-1)*a(i))/L(i); u(i)=c(i)/L(i); end L(n)=b(n)-a(n)*u(n-1); y(n)=(d(n)-y(n-1)*a(n))/L(n); x(n)=y(n); fori=(n-1): -1: 1 x(i)=y(i)-u(i)*x(i+1); end x' 计算结果 ans= 1 1 1 1 7.已知 编程求解矩阵A的QR分解. MATLAB程序 function[Q,R]=qrhouseholder(A) dim=size(A,1); R=A; Q=eye(dim); fori=1: (dim-1) x=R(i: dim,i); y=[1;zeros(dim-i,1)]; Ht=householder(x,y); H=blkdiag(eye(i-1),Ht); Q=Q*H; R=H*R; end end function[H,rho]=householder(x,y) x=x(: ); y=y(: ); iflength(x)~=length(y) error('TheColumnVectorsXandYMustHaveTheSameLength! '); end rho=-sign(x (1))*norm(x)/norm(y); y=rho*y; v=(x-y)/norm(x-y); I=eye(length(x)); H=I-2*v*v'; End 计算结果 C= 1.00001.000000 -1.00003.0000-0.50000.5000 -2.00002.00001.50000.5000 -2.00002.0000-0.50002.5000 >>[Q,R]=qrhouseholder(C) Q= -0.3162-0.7071-0.25820.5774 0.3162-0.70710.2582-0.5774 0.6325-0.0000-0.7746-0.0000 0.6325-0.00000.51640.5774 R= -3.16233.16230.47432.0555 0.0000-2.82840.3536-0.3536 0.00000.0000-1.54921.0328 -0.0000-0.0000-0.00001.1547 8.分别应用Jacobi迭代法和Gauss-Seidel迭代法求解如下方程组 Jacobi法主程序 function[x,k,succ]=Jacobi(A,b,eps,it_max) n=length(A); k=0; x=zeros(n,1); y=zeros(n,1); succ=1; while1 fori=1: n y(i)=b(i); forj=1: n ifj~=i y(i)=y(i)-A(i,j)*x(j); end end ifk==it_max succ=0;return; end y(i)=y(i)/A(i,i); end ifnorm(y-x,2) break; end x=y;k=k+1; end 计算结果 >>A=[4,-1,1;4,8,1;-2,1,5]; >>b=[7,-21,15]; >>[x,k,succ]=Jacobi(A,b,1e-6,1000) x= 0.0606 -3.1111 3.6465 k= 21 succ= 1 G-
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 大连理工大学 矩阵 上机