完整版山建大数值分析实验报告.docx

完整版山建大数值分析实验报告.docx

《完整版山建大数值分析实验报告.docx》由会员分享，可在线阅读，更多相关《完整版山建大数值分析实验报告.docx（29页珍藏版）》请在冰豆网上搜索。

完整版山建大数值分析实验报告

山东建筑大学

《数值计算A》实验报告

二级学院：

理学院

专业：

信息与计算科学

指导教师：

XXX

班级学号：

信计XXXX

姓名：

XXXX

实验一高斯消去法

（1）题目：

用选主元素法和高斯消去法求解下列方程组：

已知方程增广矩阵：

a=

2-1006

-1-3-201

-13-200

00-351

（2）原理：

（1）高斯消去法：

相对于约当消去法而言，总的来说就是将增广矩阵化为下三角阵。

（2）顺序选主元素法:

相对于高斯消去法的唯一不同是要先按当前要排列元素所在列大小进行排列。

（3）设计思想：

（1）高斯消去法：

首先让每一行的元素除以该行的主对角线元素。

然后利用此行使位于下一行主对角线以前的元素变为0，依次类推。

（2）选主元素法：

在高斯消去法的基础上，每次进行化上三角阵之前，重新排列各方程的位置。

（4）对应程序：

（a）高斯消去法

functiony=gauss1（a,b）

%高斯顺序消去法

[m,n]=size（a）;

ifm~=n

disp（'输入错误，系数矩证阵只能是方阵'）

end

ifn~=length（b）

disp（'输入错误,常数项的个数应与方程的个数相同'）

end

fork=1:

n-1

fori=k+1:

n

a（i,k）=a（i,k）/a（k,k）;

b（i）=b（i）-a（i,k）*b（k）;

forj=k+1:

n

ifa（k,k）==0

disp（'主元素为零，消去法无法继续'）;

break;

else

a（i,j）=a（i,j）-a（i,k）*a（k,j）;

end

end

end

end

b（n）=b（n）/a（n,n）;

fori=（n-1）:

-1:

1

w=0;

forj=（i+1）:

n

w=w+a（i,j）*b（j）;

end

b（i）=（b（i）-w）/a（i,i）;

end

y=b;

（b）高斯列主元消去法

functionz=gauss2（a,b,ep）

%高斯列主元素消元法

ifnargin==2

ep=0.000001

end

[m,n]=size（a）;

ifm~=n

disp（'输入错误，系数矩证阵只能是方阵'）

end

ifn~=length（b）

disp（'输入错误,常数项的个数应与方程的个数相同'）

end

fork=1:

n-1

p=a（k,k）;I=k;

fori=k:

n

ifabs（a（i,k））>abs（p）

p=a（i,k）;I=i;

end

end

ifp<=ep

z=0;

end

ifI~=k

forj=k:

n

w=a（k,j）;a（k,j）=a（I,j）;a（I,j）=w;

end

u=b（k）;b（k）=b（I）;b（I）=u;

end

fori=k+1:

n

a（i,k）=a（i,k）/a（k,k）;

b（i）=b（i）-a（i,k）*b（k）;

forj=k+1:

n

a（i,j）=a（i,j）-a（i,k）*a（k,j）;

end

end

end

b（n）=b（n）/a（n,n）;

fori=（n-1）:

-1:

1

w=0;

forj=（i+1）:

n

w=w+a（i,j）*b（j）;

end

b（i）=（b（i）-w）/a（i,i）;

end

%disp（'方程组的解为'）;

z=b;

（5）实验结果：

>>a=

2-100

-1-3-20

-13-20

00-35

>>b=

6

1

0

1

>>gauss1（a,b）

方程组的解为

ans=

35/12

-1/6

-41/24

-33/40

（6）实验体会：

主元消去法和高斯消去法的确是两个非常锻炼人编程的方法，在编写程序时，需要使用的大量的循环和分支结构，但无论是高斯消去法还是高斯列主元法，它们的原理还算不难理解，通过变成能够较好的理解它们。

实验二解线性方程组的迭代法

2.1实验目的

掌握解线性方程组的雅可比迭代和高斯-塞德尔迭代算法；

培养编程与上机调试能力.

2.2实验要求：

（1）选择一种计算机语言（Matlab）设计出雅可比（Jacobi）Gauss-Seidel、SOR迭代法，迭代法的算法程序，并且选择不同的迭代次数，观察输出结果；

（2）利用Matlab求方程组

2.3实验内容

计算书上的习题

（1）题目：

P61例2.5.1

a=[2023;181;2-315];b=[24;12;30];x0=[0;0;0];  

（2）对应程序：

Jacobi迭代法：

functionX=jacobi（a,b,X0,ep）

%Jacobi迭代法求解方程组

ifnargin==3

ep=1.0e-6;

elseifnargin<3

error

return

end

D=diag（diag（a））;

D=inv（D）;

L=tril（a,-1）;

U=triu（a,1）;

B=-D*（L+U）;

f=D*b;

X=B*X0+f;

k=1;

while（norm（X-X0,inf）>=ep）&（k<=1000）

X0=X;

X=B*X0+f;

k=k+1;

end

disp（'迭代次数为'）

k

return

高斯-塞德尔（Gauss-Seidel）：

functionz=gauss2（a,b,ep）

%高斯列主元素消元法

ifnargin==2

ep=0.000001

end

[m,n]=size（a）;

ifm~=n

disp（'输入错误，系数矩证阵只能是方阵'）

end

ifn~=length（b）

disp（'输入错误,常数项的个数应与方程的个数相同'）

end

fork=1:

n-1

p=a（k,k）;I=k;

fori=k:

n

ifabs（a（i,k））>abs（p）

p=a（i,k）;I=i;

end

end

ifp<=ep

z=0;

end

ifI~=k

forj=k:

n

w=a（k,j）;a（k,j）=a（I,j）;a（I,j）=w;

end

u=b（k）;b（k）=b（I）;b（I）=u;

end

fori=k+1:

n

a（i,k）=a（i,k）/a（k,k）;

b（i）=b（i）-a（i,k）*b（k）;

forj=k+1:

n

a（i,j）=a（i,j）-a（i,k）*a（k,j）;

end

end

end

b（n）=b（n）/a（n,n）;

fori=（n-1）:

-1:

1

w=0;

forj=（i+1）:

n

w=w+a（i,j）*b（j）;

end

b（i）=（b（i）-w）/a（i,i）;

end

%disp（'方程组的解为'）;

z=b;

gauss2（a,b,ep）

（3）实验结果：

  Jacobi迭代法：

迭代次数为

k=

6

ans=

0.33

1.33

2.67

高斯消元法：

ans=

0.15

1.94

2.37 

（4）实验体会：

Jacobi迭代法和高斯消元法都能很好的解决方程组的求解问题，在上机程中遇到了也遇到了不少问题，但最后在老师的悉心辅导下都得到了很好的解答，这两个程序使我明白了要变出好的程序就需要我们积极思考问题，勇于发现和解决问题。

实验三矩阵特征值问题计算

3.1实验目的

掌握求矩阵的特征值和主特征向量的幂法；

培养编程与上机调试能力.

3.2实验要求

（1）选择一种计算机语言设计出幂法求主特征值和相应特征向量的程序，并且选择不同的初值，观察所需的迭代次数和迭代结果.

（2）利用Matlab求特征值和特征向量

调用格式1：

eig（A）%得到特征值列向量

调用格式2：

，其中为由特征列向量构成的方阵，为由特征值构成的对角阵.%得到特征值和所对应的特征向量

3.3实验内容

计算书上的习题

（1）题目：

  P81例3.1.1

A=[246;3915;41636]x0=[1;1;1]

（2）对应程序：

乘幂法

  functiony=chm（a,x0,k）

%乘幂法求主特征值及特征向量

ifnargin<3

k=10

end

fori=1:

k

[mj]=max（abs（x0））;

y1=（1/x0（j））*x0;%标准化

x0=a*y1;

end

[mj]=max（abs（x0））

lanmda=x0（j）

u1=y1

A=[246;3915;41636]

x0=[1;1;1]

chm（A,x0,6）

经典Jacobi法

functionyy=Jeig（a,eps）

%经典Jacobi法-求全部特征值及特征向量

ifa~=a'

disp（'输入错误,应输入对称矩阵'）

return

end

n=size（a）;

ifnargin==1

eps=1.0e-6;

end

I=eye（n）;

SS=0;

fori=2:

n

forj=1:

i-1

SS=SS+a（i,j）^2;

end

end

k=1;

G=a;

R1=I;

p=0;

q=0;

while（SS>=eps）&（k<10000）

a=G;

w=0;

R=R1;

fori=2:

n

forj=1:

i-1

ifabs（a（i,j））>w

w=abs（a（i,j））;

p=i;

q=j;

end

end

end

ifa（p,p）==a（q,q）

theta=sign（a（p,q））*pi/4;

G（p,p）=a（p,p）*（cos（theta））^2+a（q,q）*（sin（theta））^2+a（p,q）*sin（2*theta）;

G（q,q）=a（p,p）*（sin（theta））^2+a（q,q）*（cos（theta））^2-a（p,q）*sin（2*theta）;

G（p,q）=（a（q,q）-a（p,p））*（cos（theta））*（sin（theta））+a（p,q）*（cos（2*theta））;

G（q,p）=G（p,q）;

fori=1:

n

R1（i,p）=R（i,p）*cos（theta）+R（i,q）*sin（theta）;

R1（i,q）=-R（i,p）*sin（theta）+R（i,q）*cos（theta）;

if（i~=p）&（i~=q）

G（p,i）=a（i,p）*（cos（theta））+a（i,q）*（sin（theta））;

G（i,p）=G（p,i）;

G（q,i）=-a（i,p）*（sin（theta））+a（i,q）*（cos（theta））;

G（i,q）=G（q,i）;

end

end

else

gasi=（a（