大连理工大学矩阵上机Word文档格式.docx
- 文档编号:16587680
- 上传时间:2022-11-24
- 格式:DOCX
- 页数:47
- 大小:153.99KB
大连理工大学矩阵上机Word文档格式.docx
《大连理工大学矩阵上机Word文档格式.docx》由会员分享,可在线阅读,更多相关《大连理工大学矩阵上机Word文档格式.docx(47页珍藏版)》请在冰豆网上搜索。
0.7499990000005
结果分析
计算机做加减法时,运算次序会影响结果,计算和时应先安排绝对值小的数参与运算,这样能取得较高的精度。
2.秦九韶算法。
已知n次多项式
,用秦九韶算法编写通用的程序计算函数在
点的值,并计算
在23点的值。
(提示:
编写程序时,输入系数向量和点
,输出结果,多项式的次数可以通过向量的长度来判断).
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);
%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);
x=zeros(n,1);
x(n)=b(n)/a(n,n);
fori=n-1:
1
x(i)=(b(i)-a(i,[i+1:
n])*x([i+1:
n]))/a(i,i);
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];
max=k;
fori=k+1:
n
ifa(i,k)>
a(max,k)
max=i;
temp=a(k,k:
n+1);
a(k,k:
n+1)=a(max,k:
a(max,k:
n+1)=temp;
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:
x(n,1)=a(n,n+1)/a(n,n);
fori=n-1:
1
sum=0;
forj=i+1:
sum=sum+x(j,1)*a(i,j);
x(i,1)=(a(i,n+1)-sum)/a(i,i);
[x]=gauss(a,b)
-0.289233816015754
0.345435715779115
-0.712811731086879
-0.220608510570529
-0.430400432704022
0.154308739838311
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;
U(1,1)=a(1,1)/L(1,1);
ifL(1,1)*U(1,1)==0
Factorizationimpossible'
)
else
forj=2:
U(1,j)=a(1,j);
L(j,1)=a(j,1)/U(1,1);
fori=2:
sum=0;
fork=1:
i-1
sum=sum+L(i,k)*U(k,i);
U(i,i)=(a(i,i)-sum)/L(i,i);
ifL(i,i)*U(i,i)==0
)
else
forj=i+1:
h=0;
s=0;
h=h+L(i,k)*U(k,j);
s=s+L(j,k)*U(k,i);
U(i,j)=1/L(i,i)*(a(i,j)-h);
L(j,i)=1/U(i,i)*(a(j,i)-s);
sum=sum+L(n,k)*U(k,n);
U(n,n)=(a(n,n)-sum)/L(n,n);
[L,U]=mylu(a)
L=
Columns1through3
100
-0.41935483870967710
0-0.3045851528384281
00-0.353872899362565
000
Columns4through6
-0.39755492585681310
0-0.1569436359494821
00-0.653976718762685
0-0.112102597106773-0.0175453757582425
-0.119266477757044-0.0833908821051642-0.0142268147798013
Columns7through9
-0.024618525643364810
-0.0199621375629645-0.09231647025909681
U=
31-130
029.5483870967742-9
0028.2587336244541
00-10
0-11-4.19354838709677
-10-3.35043668122271-1.27729257641921
75.4612710063743-31.185********5-0.451999227351748
044.6019996774714-7.17969451931716
0045.8731926371318
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)
6.1817e+13
inv(U)*inv(L)
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
ifdet(a)==0
Ü
â
u=a;
p=eye(m);
l=eye(m);
m
forj=i:
t(j)=u(j,i);
t(j)=t(j)-u(j,k)*u(k,i);
a=i;
b=abs(t(i));
ifb<
abs(t(j))
b=abs(t(j));
a=j;
ifa~=i
forj=1:
c=u(i,j);
u(i,j)=u(a,j);
u(a,j)=c;
c=p(i,j);
p(i,j)=p(a,j);
p(a,j)=c;
c=t(a);
t(a)=t(i);
t(i)=c;
u(i,i)=t(i);
m;
u(j,i)=t(j)/t(i);
u(i,j)=u(i,j)-u(i,k)*u(k,j);
end
l=tril(u,-1)+eye(m);
u=triu(u,0);
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,其中
,
.
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:
dim
ifi>
=(j+1)
L(i,j)=(A(i,j)-sum(L(i,1:
(j-1)).*L(j,1:
(j-1))))/L(j,j);
L(j,j)=sqrt(A(j,j)-sum(L(j,1:
j-1).*L(j,1:
j-1)));
L_T=L.'
;
fori=1:
ifi==1
y
(1)=B
(1)/L(1,1);
y(i)=(B(i)-sum(L(i,1:
(i-1)).*y(1:
(i-1))))/L(i,i);
fori=dim:
1
ifi==dim
x(i)=y(i)/L(i,i);
x(i)=(y(i)-sum(L((i+1):
dim,i).*x((i+1):
dim)))/L(i,i);
%disp('
(((((((('
%disp(x(i))
disp('
x'
disp(x)
x
52.2500
-38.7500
30.7500
-52.7500
6.用追赶法编制程序求解方程组Ax=b,其中
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];
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);
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):
x(i)=y(i)-u(i)*x(i+1);
7.已知
编程求解矩阵A的QR分解.
function[Q,R]=qrhouseholder(A)
dim=size(A,1);
R=A;
Q=eye(dim);
(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;
function[H,rho]=householder(x,y)
x=x(:
y=y(:
iflength(x)~=length(y)
TheColumnVectorsXandYMustHaveTheSameLength!
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;
y=zeros(n,1);
succ=1;
while1
fori=1:
y(i)=b(i);
ifj~=i
y(i)=y(i)-A(i,j)*x(j);
ifk==it_max
succ=0;
y(i)=y(i)/A(i,i);
ifnorm(y-x,2)<
eps
break;
x=y;
k=k+1;
A=[4,-1,1;
4,8,1;
-2,1,5];
b=[7,-21,15];
[x,k,succ]=Jacobi(A,b,1e-6,1000)
0.0606
-3.1111
3.6465
k=
21
succ=
G-
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 大连理工大学 矩阵 上机
![提示](https://static.bdocx.com/images/bang_tan.gif)