大连理工大学矩阵与数值分析上机作业

矩阵与数值分析上机作业

学校：大连理工大学

学院：

班级：

姓名：

学号：

授课老师：

注：编程语言Matlab

程序：

Norm.m函数

function s=Norm(x,m)

%求向量x的范数

%m取1,2,inf分别表示1,2，无穷范数

n=length(x);

s=0;

switch m

case 1 %1-范数

for i=1:n

s=s+abs(x(i));

end

case 2 %2-范数

for i=1:n

s=s+x(i)^2;

end

s=sqrt(s);

case inf %无穷-范数

s=max(abs(x));

end

计算向量x，y的范数

Test1.m

clear all;

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时

x的1-范数:2.9290；x的2-范数:1.2449； x的无穷-范数:1

y的1-范数:55； y的2-范数:19.6214； y的无穷-范数:10

n=100时

x的1-范数:5.1874；x的2-范数: 1.2787； x的无穷-范数:1

y的1-范数:5050； y的2-范数:581.6786； y的无穷-范数:100

n=1000时

x的1-范数:7.4855； x的2-范数:1.2822； x的无穷-范数:1

y的1-范数: 500500； y的2-范数:1.8271e+004；y的无穷-范数:1000

程序

Test2.m

clear all;

clc;

n=100;%区间

h=2*10^(-15)/n;%步长

x=-10^(-15):h:10^(-15);

%第一种原函数

f1=zeros(1,n+1);

for k=1:n+1

if x(k)~=0

f1(k)=log(1+x(k))/x(k); else

f1(k)=1;

end

end

subplot(2,1,1);

plot(x,f1,'-r');

axis([-10^(-15),10^(-15),-1,2]); legend('原图');

%第二种算法

f2=zeros(1,n+1);

for k=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([-10^(-15),10^(-15),-1,2]); legend('第二种算法');

运行结果：

显然第二种算法结果不准确，是因为计算机中的舍入误差造成的，当]10,10[1515--∈x 时，

x d +=1，计算机进行舍入造成d 恒等于1，结果函数值恒为1。

程序：

秦九韶算法:QinJS.m

function y=QinJS(a,x) %y 输出函数值

%a 多项式系数，由高次到零次 %x 给定点

n=length(a); s=a(1); for i=2:n s=s*x+a(i); end y=s;

计算p(x):test3.m

clear all ; clc;

x=1.6:0.2:2.4;%x=2的邻域

disp('x=2的邻域：');x

a=[1 -18 144 -672 2016 -4032 5376 -4608 2304 -512]; p=zeros(1,5); for i=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; for k=1:nk

pk(k)=QinJS(a,xk(k)); end

plot(xk,pk,'-r');

xlabel('x');ylabel('p(x)');

运行结果：

x=2的邻域：

x =

1.6000 1.8000

2.0000 2.2000 2.4000

相应多项式p值：

p =

1.0e-003 *

-0.2621 -0.0005 0 0.0005 0.2621 p(x)在 x[1.95,20.5]上的图像

程序：

LU分解，LUDecom.m

function [L,U]=LUDecom(A)

%不带列主元的LU分解

N = size(A);

n = N(1);

L=eye(n);U=zeros(n);

for i=1:n

U(1,i)=A(1,i);L(i,1)=A(i,1)/U(1,1); end

for i=2:n

for j=i:n

z=0;

for k=1:i-1

z=z+L(i,k)*U(k,j);

end

U(i,j)=A(i,j)-z;

end

for j=i+1:n

z=0;

for k=1:i-1

z=z+L(j,k)*U(k,i);

end

L(j,i)=(A(j,i)-z)/U(i,i);

end

end

PLU分解，PLUDecom.m

function [P,L,U] =PLUDecom(A)

%带列主元的LU分解

[m,m]=size(A);

U=A;

P=eye(m);

L=eye(m);

for i=1:m

for j=i:m

t(j)=U(j,i);

for k=1:i-1

t(j)=t(j)-U(j,k)*U(k,i);

end

end

a=i;b=abs(t(i));

for j=i+1:m

if b

b=abs(t(j));

a=j;

end

end

if a~=i

for j=1:m

c=U(i,j);

U(i,j)=U(a,j);

U(a,j)=c;

end

for j=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);

for j=i+1:m

U(j,i)=t(j)/t(i);

end

for j=i+1:m

for k=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

clear all;

clc;

for n=5:30

x=zeros(n,1);

A=-ones(n);

A(:,n)=ones(n,1);

for i=1:n

A(i,i)=1;

for j=(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]=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);

for i=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 =

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.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

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

当n= 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

程序：

Cholesky分解：Cholesky.m

function L=Cholesky(A)

N = size(A);

n = N(1);

L=zeros(n);

L(1,1)=sqrt(A(1,1));

for i=2:n

L(i,1)=A(i,1)/L(1,1);

end

for j=2:n

s1=0;

for k=1:j-1

s1=s1+L(j,k)^2;

end

L(j,j)=sqrt(A(j,j)-s1);

for i=j+1:n

s2=0;

for k=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

clear all;clc;

for n=10:20

A=zeros(n,n);

b=zeros(n,1);

for i=1:n

for j=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

-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

利用Cholesky分解的方程组的解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

程序：

（1）House.m

function u=House(x)

n=length(x);

e1=eye(n,1);

w=x-norm(x,2)*e1;

u=w/norm(w,2);

(2)Hou_A.m

function HA=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

clear all;

clc;

A=[1 2 3 4;

-1 2 sqrt(2) sqrt(3);

-2 2 exp(1) pi;

-sqrt(10) 2 -3 7;

0 2 7 5/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

程序：

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)); %求A的对角矩阵

L=-tril(A,-1); %求A的下三角阵

U=-triu(A,1); %求A的上三角阵

J=D\(L+U);

f=D\b;

x=J*x0+f

n=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.m

function [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); %求A的上三角阵

G=(D-L)\U;

f=(D-L)\b;

x=G*x0+f

n=1; %迭代次数

err=norm(x-x0,inf)

while(err>=eps)

x0=x;

x=G*x0+f

n=n+1;

err=norm(x-x0,inf)

if(n>=M)

disp('Warning: 迭代次数太多，可能不收敛！'); return;

end

end

解方程组，test7.m

clear all;

clc;

A=[5 -1 -3;

-1 2 4;

-3 4 15];

b=[-2;1;10];

disp('精确解');

x=A\b

disp('迭代初始值');

x0=[0;0;0]

disp('Jacobi迭代过程：');

[xj,nj]=Jaccobi(A,b,x0);

disp('Jacobi最终迭代结果：');

xj

disp('迭代次数');

nj

disp('Gauss-Seidel迭代过程：'); [xg,ng]=Gauss_Seidel(A,b,x0);

disp('Gauss-Seidel最终迭代结果：'); xg

disp('迭代次数');

ng

运行结果：

精确解

x =

-0.0820

-1.8033

1.1311

迭代初始值

x0 =

Jacobi迭代过程：

x =

-0.4000

0.5000

0.6667

err =

0.6667

x =

0.1000

-1.0333

0.4533

err =

1.5333

...

x =

-0.0820

-1.8033

1.1311

err =

9.6603e-006

Jacobi最终迭代结果：

xj =

-0.0820

-1.8033

1.1311

迭代次数

nj =

281

Gauss-Seidel迭代过程：

x =

-0.4000

0.3000

0.5067

err =

0.5067

x =

-0.0360

-0.5313

0.8012

err =

0.8313

x =

-0.0256

-1.1151

0.9589

err =

0.5838

...

x =

-0.0820

-1.8033

1.1311

err =

9.4021e-006

Gauss-Seidel最终迭代结果：xg =

-0.0820

-1.8033

1.1311

迭代次数

ng =

20

程序：

Newton迭代法：Newtoniter.m

function [x,iter,fvalue]=Newtoniter(f,df,x0,eps,maxiter) %牛顿法 x得到的近似解

%iter迭代次数

%fvalue函数在x处的值

%f，df被求的非线性方程及导函数

%x0初始值

%eps 允许误差限

%maxiter 最大迭代次数

fvalue=subs(f,x0);dfvalue=subs(df,x0);

for iter=1:maxiter

x=x0-fvalue/dfvalue

err=abs(x-x0)

x0=x;

fvalue=subs(f,x0)

dfvalue=subs(df,x0);

if(err

end

弦截法：secant.m

function [x,iter,fvalue]=secant(f,x0,x1,eps,maxiter)

%弦截法 x得到的近似解

%iter迭代次数

%fvalue函数在x处的值

%f被求的非线性方程

%x0，x1初始值

%eps 允许误差限

%maxiter 最大迭代次数

fvalue0=subs(f,x0);fvalue=subs(f,x1);

for iter=1:maxiter

x=x1-fvalue*(x1-x0)/(fvalue-fvalue0)

err=abs(x-x1)

x0=x1;x1=x;

fvalue0=subs(f,x0);fvalue=subs(f,x1)

if(err

end

求方程的实根：test8.m

clear all；

clc;

syms x