矩阵与数值分析上机作业满分.docx

矩阵与数值分析上机作业满分.docx

《矩阵与数值分析上机作业满分.docx》由会员分享，可在线阅读，更多相关《矩阵与数值分析上机作业满分.docx（36页珍藏版）》请在冰豆网上搜索。

矩阵与数值分析上机作业满分

第一题

考虑计算给定向量的范数：

输入向量

输出

。

请编制一个通用程序，并用你编制的程序计算如下向量的范数：

对n=10，100，1000甚至更大的n计算其范数，你会发现什么结果?

你能否修改你的程序使得计算结果相对精确呢？

代码：

functionProblem1（）

clc

N=input（'inputN='）;

X=zeros（1,N）;

Y=1:

N;

fori=1:

N

X（i）=1/i;

end

x1=0;

x2=0;

y1=0;

y2=0;

fori=1:

N

x1=x1+abs（X（i））;

x2=x2+X（i）*X（i）;

y1=y1+abs（Y（i））;

y2=y2+Y（i）*Y（i）;

end

x1;

x2=sqrt（x2）;

x3=max（abs（X））;

y1;

y2=sqrt（y2）;

y3=max（abs（Y））;

fprintf（'x的1范数是%.2f\n',x1）;

fprintf（'x的2范数是%.2f\n',x2）;

fprintf（'x的无穷范数是%.2f\n',x3）;

fprintf（'y的1范数是%.2f\n',y1）;

fprintf（'y的2范数是%.2f\n',y2）;

fprintf（'y的无穷范数是%.2f\n',y3）;

结果：

inputN=10

x的1范数是2.93

x的2范数是1.24

x的无穷范数是1.00

y的1范数是55.00

y的2范数是19.62

y的无穷范数是10.00

inputN=100

x的1范数是5.19

x的2范数是1.28

x的无穷范数是1.00

y的1范数是5050.00

y的2范数是581.68

y的无穷范数是100.00

inputN=1000

x的1范数是7.49

x的2范数是1.28

x的无穷范数是1.00

y的1范数是500500.00

y的2范数是18271.11

y的无穷范数是1000.00

第二题

考虑

其中定义f（0）=1，此时f（x）是连续函数。

用此公式计算当

时的函数值，画出图像。

另一方面，考虑下面算法：

d=1+x

ifd=1then

y=1

else

y=lnd/（d-1）

endif

用此算法计算

时的函数值，画出图像。

比较一下发生了什么？

代码：

functionproblem2（）

clc

N=1000;

n=2*10^（-15）/N;

%公式计算

f=zeros（1,N+1）;

t=1;

fori=-10^（-15）:

n:

10^（-15）

if（i~=0）

f（t）=log（1+i）/i;

else

f（t）=1;

end

t=t+1;

end

%算法计算

a=zeros（1,N+1）;

t=1;

fori=-10^（-15）:

n:

10^（-15）

d=1+i;

if（d~=1）

a（t）=log（d）/（d-1）;

else

a（t）=1;

end

t=t+1;

end

%画图，左边是公式算的，右边是算法算的

subplot（1,2,1）;

plot（-10^（-15）:

n:

10^（-15）,f）;

holdon

subplot（1,2,2）;

plot（-10^（-15）:

n:

10^（-15）,a）;

holdon

结果：

第三题

首先编制一个利用秦九韶算法计算一个多项式在给定点的函数值的通用程序，你的程序包括输入多项式的系数以和给定点，输出函数值。

利用你编制的程序计算

在x=2邻域附近的值。

画出

上的图像。

代码：

functionproblem3（）

%秦九韶算法

clc

n=9;

c=[1,-18,144,-672,2016,-4032,5376,-4608,2304,-512];

N=100;

m=（2.05-1.95）/N;

y=zeros（1,N+1）;

ydex=0;

forx=1.95:

m:

2.05

ydex=ydex+1;

dex=1;

temp=c（dex）;

whiledex<=n

dex=dex+1;

temp=temp*x+c（dex）;

end

y（ydex）=temp;

end

plot（1.95:

m:

2.05,y）;

结果：

第四题

编制计算给定矩阵A的LU分解和PLU分解的通用程序，然后用你编制的程序完成下面两个计算任务：

（1）考虑

自己取定

，并计算b=Ax。

然后用你编制的不选主元和列主元的Gauss消去法求解该方程组，记你计算出的解为

。

对n从5到30估计计算解的精度。

（2）对n从5到30计算其逆矩阵。

代码：

%LU分解

functionproblem4（）

clc

%生成矩阵A

forn=5:

30

A=-ones（n）;

A（:

n）=ones（n,1）;

foria=1:

n

A（ia,ia）=1;

forja=（ia+1）:

（n-1）

A（ia,ja）=0;

end

end

x=rand（n,1）;

b=A*x;

%LU分解

[n,n]=size（A）;

L=eye（n）;

M=eye（n）;

fori=1:

n-1

ifA（i,i）==0

fprintf（'Error:

A（%d,%d）=0!

\n',i,i）;return;

end

forj=i+1:

n

L（j,i）=A（j,i）/A（i,i）;

end

fork=i+1:

n

forl=i:

n

A（k,l）=A（k,l）-L（k,i）*A（i,l）;

end

end

end

U=A;

M=inv（U）*inv（L）*M

y1=inv（L）*b;

x1=inv（U）*y1;

error=norm（x-x1）

end

return;

%PLU分解

functionproblem4（）

clc

%生成矩阵A

n=5;

A=-ones（n）;

A（:

n）=ones（n,1）;

foria=1:

n

A（ia,ia）=1;

forja=（ia+1）:

（n-1）

A（ia,ja）=0;

end

end

x=rand（n,1）;

b=A*x;

%PLU分解

[n,n]=size（A）;

Q=eye（n）;

R=eye（n）;

S=eye（n）;

L=zeros（n,n,n-1）;

foril=1:

（n-1）

L（:

:

il）=eye（n）;

end

P=zeros（n,n,n-1）;

foril=1:

（n-1）

P（:

:

il）=eye（n）;

end

fori=1:

n-1

m=i;

a=abs（A（i,i））;

fork=i+1:

n

ifabs（A（k,i））>a

a=abs（A（k,i））;

m=k;

end

end

A（[i,m],:

）=A（[m,i],:

）;

P（[i,m],:

i）=P（[m,i],:

i）;

forj=i+1:

n

L（j,i,i）=-A（j,i）/A（i,i）;

forl=i:

n

A（j,l）=A（j,l）+L（j,i,i）*A（i,l）;

end

end

end

U=A;

U

forx=1:

n-1

Q=L（:

:

x）;

fory=x+1:

n-1

Q=P（:

:

y）*Q*P（:

:

y）;

end

S=S*inv（Q）;

end

%L矩阵

S

forz=1:

（n-1）

R=P（:

:

z）*R;

end

%P矩阵

R

b1=R*b;

b2=inv（S）*b1;

x2=inv（U）*b2;

M=inv（U）*inv（S）*R

error=norm（x-x2）

结果：

由于篇幅有限仅列出n=5时的精度和逆矩阵

LU分解

M=

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

error=

8.8818e-16

PLU分解

M=

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

error=

7.8516

第五题

编制计算对称正定阵的Cholesky分解的通用程序，并用你编制的程序计算Ax=b，其中

。

b可以由你自己取定，对n从10到20验证程序的可靠性。

代码：

functionproblem5（）

%cholesky分解

clc

%生成矩阵A

N=10;

A=zeros（N）;

foria=1:

N

forja=1:

N

A（ia,ja）=1/（ia+ja-1）;

end

end

L=chol（A）;

L

结果：

由于篇幅有限仅列出n=10时的L矩阵

L=

1.00000.50000.33330.25000.20000.16670.14290.12500.11110.1000

00.28870.28870.25980.23090.20620.18560.16840.15400.1417

000.07450.11180.12780.13310.13310.13040.12650.1220

0000.01890.03780.05250.06300.07020.07480.0777

00000.00480.01190.01950.02650.03260.0378

000000.00120.00360.00680.01030.0139

0000000.00030.00110.00220.0038

00000000.00010.00030.0007

000000000.00000.0001

0000000000.0000

第六题

（1）编制程序House（x），其作用是对输入的向量x，输出单位向量u使得

。

（2）编制Householder变换阵

乘以

的程序HA，注意，你的程序并不显式的计算出H。

（3）考虑矩阵

用你编制的程序计算H使得HA的第一列为

的形式，并将HA的结果显示。

代码：

functionproblem6（）

clc

x=[1;-1;-2;-10^（1/2）;0];

ifx

（1）>0

sg=1;

else

sg=-1;

end

Nx=length（x）;

e1=zeros（Nx,1）;

e1

（1）=1;

sum=0;

fori=1:

Nx

sum=sum+x（i）*x（i）;

end

x2=sqrt（sum）;

clearsum;

w=x-sg*x2*e1;

u=w/sqrt（w'*w）;

u

%生成矩阵A

A=[1,2,3,4;-1,3,2^（1/2）,3^（1/2）;-2,-2,exp

（1）,pi;-10^（1/2）,2,-3,7;0,2,7,5/2];

H=eye（Nx）-2*u*u';

H*A

结果：

u=

-0.6124

-0.2041

-0.4082

-0.6455

0

ans=

4.0000-0.83111.4090-6.5378

0.00002.05630.8839-1.7805

0.0000-3.88741.6576-3.8836

0.0000-0.9843-4.6770-4.1078

02.00007.00002.5000

第七题

用Jacobi和Gauss-Seidel迭代求解下面的方程组，输出迭代每一步的误差

:

代码：

functionproblem7（）

clc

inter=10;%迭代次数

rx=[69*7/61-8,69*3/（61*2）-7/2,69/61]';

%根据题意生成A,b

A=[5-1-3;-124;-3415];

b=[-2110]';

N=length（b）;

x=ones（N,1）;

%Jacobi迭代法

L=zeros（N,N）;

D=L;

U=L;

foria=1:

N

D（ia,ia）=A（ia,ia）;

ifia

foriu=（ia+1）:

N

U（ia,iu）=A（ia,iu）;

end

end

ifia>1

foril=1:

（ia-1）

L（ia,il）=-A（ia,il）;

end

end

end

forintt=1:

inter

x=inv（D）*（L+U）*x+inv（D）*b;

wuca1=0;

%求误差

forix=1:

N

wuca1=wuca1+（x（ix）-rx（ix））^2;

end

wuca1

end

%G-S迭代法

forintt=1:

inter

x=inv（D-L）*U*x+inv（D-L）*b;

wuca2=0;

%求误差

forix=1:

N

wuca2=wuca2+（x（ix）-rx（ix））^2;

end

wuca2

end

结果:

wuca1=

24.6036

wuca1=

12.3259

wuca1=

1.9276

wuca1=

5.4520

wuca1=

15.9787

wuca1=

10.0748

wuca1=

3.6731

wuca1=

6.2770

wuca1=

11.9670

wuca1=

8.9591

wuca2=

5.5096

wuca2=

9.1425

wuca2=

6.9372

wuca2=

8.0757

wuca2=

7.4675

wuca2=

7.7807

wuca2=

7.6170

wuca2=

7.7018

wuca2=

7.6577

wuca2=

7.6806

第八题

取不同的初值用Newton迭代以和弦截法求方程

的实根，列表或者画图说明收敛速度。

代码：

functionproblem8（）

%inNew记录Newton迭代法的收敛速度

%inS记录弦截法的收敛次数

clc

chuzhi=4;%初值

f=[];

%Newton迭代

x=chuzhi;

forinNew=1:

inter

xx=x-（x^3+2*x^2+10*x-100）/（3*x^2+4*x+10）;

f（inNew）=x;

ifabs（xx-x）<10^（-5）

break;

end

x=xx;

end

subplot（1,2,1）;

plot（1:

inNew,f）;

holdon

clearf;

inNew

xx

%弦截法

xk=chuzhi;

xkr=chuzhi-1;

forinS=1:

inter

xka=xk-（xk^3+2*xk^2+10*xk-100）*（xk-xkr）/（（xk^3+2*xk^2+10*xk-100）-（xkr^3+2*xkr^2+10*xkr-100））;

f（inS）=xk;

ifabs（xka-xk）<10^（-5）%精度10^（-5）

break;

end

xkr=xk;

xk=xka;

end

subplot（1,2,2）;

plot（1:

inS,f）;

holdon

inS

xka

结果：

inNew=

7

xx=

3.4606

inS=

10

xka=

3.4606

结果：

第九题

用二分法求方程

在区间

上的所有根。

代码：

functionproblem9（）

%二分法求解方程

clc

%查看根的个数

N=50;

yk=zeros（1,N+1）;

h=4*pi/N;

e=[];

c=[];

i=1;

forti=0:

h:

4*pi

e（i）=-2/（exp（ti））;

c（i）=cos（ti）;

i=i+1;

end

plot（0:

h:

4*pi,e）;

holdon

plot（0:

h:

4*pi,c）;

%根的取值范围

min=[14710];

max=[25811];

accr=zeros（1,4）;

acc=0.0001;%精度

forr=1:

4

while（max（r）-min（r）>acc）

ymax=exp（max（r））*cos（max（r））+2;

temp=（max（r）+min（r））/2;

ytemp=exp（temp）*cos（temp）+2;

ifytemp*ymax==0

break;

else

ifytemp*ymax>0

max（r）=temp;

else

min（r）=temp;

end

end

end

accr（r）=temp;

end

accr

结果：

accr=

1.88074.69407.854810.9955

第十题

考虑函数

。

用等距节点作

的Newton插值，画出插值多项式以和

的图像，观察收敛性。

代码：

functionf=Newton（x,y）

symst;

n=length（x）;

c（1:

n）=0;

f=y

（1）;

y1=0;

p=1;

fori=1:

n-1

forj=i+1:

n

y1（j）=（y（j）-y（i））/（x（j）-x（i））;

end

c（i）=y1（i+1）;

p=p*（t-x（i））;

f=f+c（i）*p;

y=y1;

end

f=simplify（f）;

end

所以考虑函数f（x）=sin（πx），x∈[0,1]，取二次、三次、四次插值函数。

二次：

>>x=linspace（0,1,3）

>>y=sin（pi*x）;

>>f2=Newton（x,y）

f2=-4*t*（t-1）

三次：

>>x=linspace（0,1,4）;

>>y=sin（pi*x）;

>>f3=Newton（x,y）

f3=-（9*3^（1/2）*t*（t-1））/4

四次：

>>x=linspace（0,1,5）;

>>y=sin（pi*x）;

>>f4=Newton（x,y）

f4=

（t*（32968342914557536*t^3-65936685829115088*t^2+5181631624232978*t+481984*2^（1/2）+23187））/90992

结果：

第十一题

对函数

，取不同的节点数n，用等距节点lagrange插值，观察Runge现象。

代码：

functionp=Lagrange（n,t）

x=linspace（-5,5,n）;

y=1./（1+x.^2）;

m=length（t）;

fori=1:

m

s=0;

fork=1:

n

l（k）=1;

forj=1:

n

ifj~=k

l（k）=l（k）*（t（i）-x（j））/（x（k）-x（j））;

end

end

s=l（k）*y（k）+s;

end

p（i）=s;

end

end

>>t=-5:

0.1:

5;

>>p=1./（1+x.^2）;

>>plot（x,p）

>>holdon

>>p5=Lagrange（5,x）;

>>plot（x,p5,'r'）

>>p10=Lagrange（10,x）;

>>plot（x,p10,'g'）

>>p15=Lagrange（15,x）;

>>plot（x,p15,'m'）

>>p20=Lagrange（20,x）;

>>plot（x,p20,'y'）

结果：

第十二题

令

，考虑积分

。

区间分为50，100,200,500，1000等，分别用复合梯形以和复合Simpson积分公式计算积分值，将数值积分的结果与精确值比较，列表说明误差的收敛性。

代码：

functionproblen12（）

clc

N=input（'请输入N='）;%改变区间份数

h=2*pi/N;

%复化梯形公式

x=0:

h:

2*pi;

sum=0;

forin=2:

N

sum=sum+exp（3*x（in））*cos（pi*x（in））;

end

Tn=2*pi*（1+2*sum+exp（6*pi）*cos（pi*2*pi））/（2*N）

%Simapson

sum1=0;

sum2=sum;

forsi=1:

N

xx=（h/2）+（si-1）*h;

sum1=sum1+exp（3*xx）*cos（pi*xx）;

end

Sn=2*pi*（1+4*sum1+2*sum2+exp（6*pi）*cos（pi*2*pi））/（6*N）

结果：

请输入N=50

Tn=

3.5125e+07

Sn=

3.5231e+07

请输入N=100

Tn=

3.5205e+07

Sn=

3.5232e+07

请输入N=200

Tn=

3.5226e+07

Sn=

3.5232e+07

请输入N=500

Tn=

3.5231e+07

Sn=

3.5232e+07

请输入N=1000

Tn=

3.5232e+07

Sn=

3.5232e+07

第十三题

分别用2点，3点以和5点的Gauss型积分公式计算如下定积分：

（1）

（2）

代码：

functiony=Chebyshev（f,n）

%使用Chebyshev多项式进行Gauss型数值积分

%输入：

f为表达式，n为Gauss点数

%输出：

y为积分值

x=zeros（1,n）;

A=zeros（1,n）;

fori=1:

n

x（i）=cos（（2*n-1）*pi/（2*n））;

A（i）=pi/n;

end

y=sum（A（1:

n）.*subs（f,x（1:

n）））;

y=double（y）;

end

functiony=Lengder（f,a,b,x,w）

%使用Lengder