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

大 连 理 工 大 学

矩阵与数值分析上机作业

课程名称： 矩阵与数值分析

研究生姓名：

交作业日时间：

2016 年

月 20日

12

1.1程序：

Clear all ;

n=input(' 请输入向量的长度n ：' ) for i=1:n; v(i)=1/i; end

Y1=norm(v,1) Y2=norm(v,2) Y3=norm(v,inf)

1.2结果

n=10 Y1 =2.9290 Y2 =1.2449 Y3 =1

n=100 Y1 =5.1874

Y2 =1.2787 Y3 =1

n=1000 Y1 =7.4855

Y2 =1.2822 Y3 =1

N=10000 Y1 =9.7876

Y2 =1.2825 Y3 =1

1.3 分析

一范数逐渐递增，随着n 的增加，范数的增加速度减小；二范数随着n 的增加，逐渐趋于定值，无群范数都是1.

2.1程序

clear all ; x(1)=-10^-15; dx=10^-18; L=2*10^3; for i=1:L

y1(i)=log(1+x(i))/x(i); d=1+x(i); if d == 1 y2(i)=1; else

y2(i)=log(d)/(d-1); end

x(i+1)=x(i)+dx; end

x=x(1:length(x)-1); plot(x,y1,'r' ); hold on plot(x,y2);

2.2 结果

2.3 分析

红色的曲线代表未考虑题中算法时的情况，如果考虑题中的算法则数值大小始终为1，这主要是由于大数加小数的原因。

第3题

3.1 程序

clear all ;

A=[1 -18 144 -672 2016 -4032 5376 -4608 2304 -512]; x=1.95:0.005:2.05; for i=1:length(x); y1(i)=f(A,x(i)); y2(i)=(x(i)-2)^9; end figure(3); plot(x,y1); hold on ;

plot(x,y2,'r' );

F.m 文件

function y=f(A,x) y=A(1);

for i=2:length(A); y=x*y+A(i); end ;

3.2 结果

第4题

4.1 程序

clear all ;

n=input(' 请输入向量的长度n ：' )

A=2*eye(n)-tril(ones(n,n),0); for i=1:n A(i,n)=1; end

n=length(A); U=A; e=eye(n); for i=1:n-1

[max_data,max_index]=max(abs(U(i:n,i))); e0=eye(n);

max_index=max_index+i-1; U=e0*U; e1=eye(n); for j=i+1:n

e1(j,i)=-U(j,i)/U(i,i); end U=e1*U;

P{i}=e0;%把变换矩阵存到P 中 L{i}=e1; e=e1*e0*e; end

for k=1:n-2 Ldot{k}=L{k}; for i=k+1:n-1

Ldot{k}=P{i}*Ldot{k}*P{i}; end end

Ldot{n-1}=L{n-1}; LL=eye(n); PP=eye(n); for i=1:n-1 PP=P{i}*PP; LL=Ldot{i}*LL; end

b=ones(n,2); b=e*b; %解方程 x=zeros(n,1); x(n)=b(n)/U(n,n); for i=n-1:-1:1

x(i)=(b(i)-U(i,:)*x)/U(i,i); end

X=U^-1*e^-1*eye(n);%计算逆矩阵 AN=X';

result2{n-4,1}=AN; result1{n-4,1}=x; fprintf('%d:\n',n) fprintf('%d ',AN);

4.2 结果

n=5

1.0625 0.0625 0.0625 0.0625

n=10

-0.875 1.125 0.125 0.125

-0.75 -0.75 1.25 0.25 -0.25

-0.5 -0.5 -0.5 1.5 -0.5

-0.0625 -0.0625 -0.0625 -0.0625 0.0625

-0.0625 -0.125

1.0625 -0.875 -0.75 0.0625 1.125 0.0625 0.125 0.0625 0.125

-0.75 1.25 0.25

-0.5 -0.5 -0.5 1.5 -0.5 -0.5 -0.5 -0.5 1.5

-0.0625 1.0625 -0.875 -0.75 -0.0625 0.0625 1.125 -0.0625 0.0625 0.125 -0.0625 0.0625 0.125

-0.75 1.25 0.25

-0.5 -0.5 -0.5 1.5 -0.5 -0.5 -0.5 -0.5 1.5 -0.5

-0.0625 -0.0625 -0.0625 -0.0625 0.0625 -0.0625 -0.0625 -0.0625 -0.0625 0.0625

-0.0625 -0.125 -0.25 1.0625 -0.875 -0.75 0.0625 1.125 0.0625 0.125 0.0625 0.125

-0.75 1.25 0.25

0.0625 -0.0625 -0.125 -0.25 -0.0625 1.0625 -0.875 -0.75 -0.0625 0.0625 1.125 -0.0625 0.0625 0.125 -0.0625 0.0625 0.125

-0.75 1.25 0.25

-0.0625 -0.125 -0.25 -0.5 0.0625 -0.0625 -0.125 -0.25

同样的方法可以算出n=20,n=30时的结果，这里就不罗列了。

第5题

5.1 程序

clear all ;

n=input(' 请输入向量的长度n ：10至20' )

for i=1:n for j=1:n

a(i,j)=1/(i+j-1); end end for j=1:n sum=0; for k=1:j-1

sum = sum + l(j,k)^2; end

l(j,j)=sqrt(a(j,j)-sum); for i=j+1:n sum=0; for k=1:j-1

sum =sum + l(i,k)*l(j,k); end

l(i,j)=(a(i,j)-sum)/l(j,j); end end

b=ones(n,1); y=zeros(n,1); y(n)=b(n)/l(n,n); for i=n-1:-1:1

y(i)=(b(i)-l(i,:)*y)/l(i,i); end l=l'; x=zeros(n,1); x(n)=y(n)/l(n,n); for i=n-1:-1:1

x(i)=(y(i)-l(i,:)*x)/l(i,i); end

fprintf('%d\t',x); fprintf('\n');

5.2 结果

n=10

n=11

n=12

n=13

n=14

n=15

n=16 445748685

688

23

.85

.6

.5

.9

.7

.8

.6

69820595. -[**************]90-77484610-80547115-82914903-72327488-62111481-60010758-1.76915E -1.68795E 08

9.9

5

111

920

907

690

037

946

+11

+11

n=17

n=18

n=19

n=20

-746517.83111493.4-[***********][***********][***********][**************]30

-[**************]40-549934613.06265E+3.18327E+3.23484E+2.86169E+2.48062E+2.55037E+7.40163E+6.93865E+97

7

592

12

12

12

12

12

12

12

12

152439746-1.17236E 7.93546E+-5.20373E -5.40791E -5.42806E -4.8713E+-4.26882E -4.65382E -1.3123E+-1.21156E 43

+11

11

+13

+13

+13

13

+13

+13

14

+14

-761846207.35352E+-6.11155E 4.7524E+14.93812E+4.89563E+4.46792E+3.97243E+4.57473E+

1.22E+15 1.11E+15

048

11

+12

4

14

14

14

14

14

2.18036E+-2.70378E 2.80302E+

-2.62E+15 -2.72E+15 -2.66E+15 -2.48E+15 -2.25E+15 -2.72E+15 -6.64E+15 -5.98E+15

11

+12

13

-8.11E+13 9.24E+15 9.59E+15 9.27E+15 8.90E+15 8.30E+15 1.04E+16 2.17E+16 1.94E+16

+11

12

-2.17E+16 -2.25E+16 -2.14E+16 -2.14E+16 -2.07E+16 -2.60E+16 -4.18E+16 -3.74E+16

11

+12

14

3.40E+16 3.53E+16 3.34E+16 3.52E+16 3.56E+16 4.28E+16 4.06E+16 3.77E+16

+11

12

+14

-3.55E+16 -3.68E+16 -3.47E+16 -3.89E+16 -4.12E+16 -4.38E+16 -4.89E+15 -1.01E+16

30

11

13

1.50E+15 1.55E+15 3.06E+15 -2.94E+15 -6.68E+15 1.72E+15 1.63E+16 1.04E+16

2.17168E+-7.94676E

12

14

4.80384E+-3.1296E+

14

14

2.304E+14 -9.01E+15 1.40E+17 1.19E+17

1.99E+15 -8.07E+16 -6.35E+16

1.70E+16 1.09E+16

7.5453E+1

4

1.48E+16 -8.62E+16 -7.98E+16

+14 1.58892E+

-2.50E+15 -2.19E+15 -8.57E+15 -2.64E+16 -2.01E+16 5.08E+15 7.29E+15 -1.45E+15 3.84E+16 3.85E+16

+13 1.08203E+

-8.98E+15 -9.29E+15 -1.03E+16 -9.10E+15 -7.78E+15 -4.50E+15 -1.23E+16 -1.58E+16

+12

4

2.35E+16 2.44E+16 2.37E+16 2.73E+16 2.94E+16 2.40E+16 -1.64E+16 -5.28E+15

7.13565E+-5.87483E

449059794-3.52275E 1.3792E+1-1.99261E 7.40507E+-1.83339E 3.69292E+-8.64269E 1.51789E+-3.70513E 6.12295E+

第6题

6.1 程序

clear all ;

A=[ 1 2 3 4; -1 3 sqrt(2) sqrt(3); -2 2 exp(1) pi; -sqrt(10) 2 -3 7; 0 2 7 5/2;]; U=f61(A(:,2)); HA=f62(U,A);

f.m 文件

function U=f61(x) e1=eye(length(x),1);

U=x-sign(x(1))*sqrt(dot(x,x))*e1; U=U./sqrt(dot(U,U)); function HA=f62(U,A) HA=A-2*U*U'*A;

6.2 结果

第7题

7.1 程序

clear all ; max=1000; x(1,:)=[1 2 3]; for i=1:max

x(i+1,1)=(-2+x(i,2)+3*x(i,3))/5; x(i+1,2)=(1+x(i,1)-4*x(i,3))/2; x(i+1,3)=(10+3*x(i,1)-4*x(i,2))/15;

err(i)=sqrt(dot(x(i+1,:)-x(i,:),x(i+1,:)-x(i,:))); if (err(i)

x(i+1,1)=(-2+x(i,2)+3*x(i,3))/5; x(i+1,2)=(1+x(i+1,1)-4*x(i,3))/2; x(i+1,3)=(10+3*x(i+1,1)-4*x(i+1,2))/15;

err(i)=sqrt(dot(x(i+1,:)-x(i,:),x(i+1,:)-x(i,:))); if (err(i)

7.2 结果

误差越来越小。

第8题

8.1 程序

clear all ; max=100; x(1)=1; for i=1:max

x(i+1)=x(i)-(x(i)^3+2*x(i)^2+10*x(i)-100)/(3*x(i)^2+4*x(i)+10); if (abs(x(i+1)-x(i))

x(i+1)=x(i)-(x(i)^3+2*x(i)^2+10*x(i)-100)/(x(i)^3+2*x(i)^2+10*x(i)-(x(i-1)^3+2*x(i-1)^2+10*x(i-1)))*(x(i)-x(i-1)); if (abs(x(i+1)-x(i))

break ; end end hold on plot(x,'r' );

8.2结果

8.3 分析

由计算结果可知，弦截法的收敛速度比牛顿法的收敛速度快。

第9题

9.1 程序

clear all ; f(0,4*pi);

f.m 文件

function []= f( l,r ) if (r-l

fprintf('%g,',(r-l)/2+l); return end

if f9(l)*f9(r)

if f9((l+r)/2+l)*f9(l)

f(l,(l+r)/2+l); end

if f9((l+r)/2+l)*f9(r)

if f9(l)*f9(r) > 0 f(l,(l+r)/2+l); f((l+r)/2+l,r); end

9.2 结果

X=4.71239、8.24668、17.6715、14.1372

第10题

10.1 程序

clear all ; n=3;% 节点个数 Xj=0:1/n:1; y=sin(pi*Xj); for i=1:n+1 f(i,1)=y(i); end

for j=2:n+1 for i=1:n-j+2; dx=((j-1)/n);

f(i,j)=(f(i+1,j-1)-f(i,j-1))/dx; end end

for i=1:n+1 a(i)=f(1,i); end

x=0:0.001:1; for i=1:1/0.001+1; y1(i)=sin(pi*x(i)); y2(i)=f10(a,Xj,x(i)); end

figure(10); plot(x,y1,'r' );

hold on ; plot(x,y2,'' );

10.2结果

10.3分析

有图像可知插值函数的值已经很接近原函数的值了。

第11题

11.1 程序

clear all ;

n=input(' 请输入n ：' ) % n代表节点

Xj=-5:1/n:5; Yj=1./(1.+Xj.^2);

x=-5:0.01:5; for i=1:10/0.01+1; y1(i)=1/(1+x(i)^2); y2(i)=f(Yj,Xj,x(i)); end

figure(11); plot(x,y1,'r' ); hold on ;

plot(x,y2,,'' );

f.m 文件

function y=f(Yj,Xj,x) y=0;

for i=1:length(Yj) l=1;

for j=1:length(Xj) if i==j continue ; end

l=l*(x-Xj(j))/(Xj(i)-Xj(j)); end

y=y+Yj(i)*l; end

11.2 结果

从左往右n 依次取1、2、3、4

11.3 分析

随着n 的不断增加，插值越来越接近真实值。

第12题

12.1 程序

clear all ;

n=input(' 请输入n ：' ) %n=50 100 200 500 1000

x=0:2*pi/n:2*pi; for i=1:n+1 X=x(i);

Y(i)=exp(3*X)*cos(pi*X); end Y(1)=1; for i=1:n

X=(x(i)+x(i+1))/2; Y2(i)=exp(3*X)*cos(pi*X); end

T=(x(n+1)-x(1))/(2*n)*(Y(1)+2*sum(Y(2:n))+Y(n+1));

S=(x(n+1)-x(1))/(6*n)*(Y(1)+4*sum(Y2)+2*sum(Y(2:n))+Y(n+1));

第13题

13.1 程序

clear all ;

G2=f(1/sqrt(3))+f(-1/sqrt(3));

G3=0.555555556*f(-0.7745966692)+0.555555556*f(+0.7745966692)+0.[1**********]*f(0);

G5=0.2369268851*f(-0.9061798459)+0.2369268851*f(+0.9061798459)+... 0.4786286705*f(-0.5384693101)+0.4786286705*f(+0.5384693101)+... 0.5688888889*f(0);

f.m 文件1

function y=f(x) y=x^2/sqrt(1-x^2);

f.m 文件2

function y=f(x) x=pi/2*(x+1)/2; y=sin(x)/x;

13.2

第14题

14.1 程序

clear all ; U0=2;%初值

step=0.2;%step表示步长 U1(1)=U0; Len=1;

for i=1:floor(Len/step)

U1(i+1)=U1(i)+step*f14(U1(i),(i-1)*step); end

% 改进的Euler 法 U2(1)=U0;

for i=1:floor(Len/step)

U2(i+1)=U2(i)+step/2*(f14(U2(i),(i-1)*step)+f14(U1(i+1),(i+1-1)*step)); end

% Runge_KUtta

U3(1)=U0;

for i=1:floor(Len/step) t=(i-1)*step; k1=f14(U3(i),t);

k2=f14(U3(i)+step/2*k1,t+step/2); k3=f14(U3(i)+step/2*k2,t+step/2); k4=f14(U3(i)+step*k3,t+step);

U3(i+1)=U3(i)+step/6*(k1+2*k2+2*k3+k4); end

for i=1:floor(Len/step)+1 t(i)=(i-1)*step; end figure(14) plot(t,U1,'' ); plot(t,U2,'g' ); plot(t,U3,'r' );

f.m 文件

function y=f(x,t) y=1/((t+1)^2)*(t*x-x^2);

14.2 结果