《数值计算方法》(第二版) 数值计算实验 matlab 有截图

数值计算方法

（第二版）

数值实验

学院信息学院专业计算机应用技术姓名张书涵学号12077503005

任课老师孙云平使用软件MA TLAB 7.0

2012 年12 月15日

一、实验目的和意义

1、通过上机学会matlab数学软件的使用，根据算法2.2的思想和步骤，实现解三对角方程组的追赶法。

2、运用所学的三种迭代法即Jacobi迭代法、Gauss-Seidel迭代法、共轭梯度法算法的思想和步骤，解各类线性方程组，编出算法程序。

3、根据插值法的思想，掌握样条插值的步骤和思想，并根据样本值用matlab软件绘制图形，熟悉了matlab的绘图函数。

4、只有通过动手编程，才能意识到自己对课本上知识的欠缺。通过本次实验，复习本学期所学知识，查漏补缺。

二、实验环境

matlab7.0

三、实验结果

第一章误差

1.利用循环语句，计算数列√5，√√5，√√√5，……的极限，要求误差小于10^-8。程序设计如下：

x =sqrt(5);

m=0;

while abs(x-m)>1e-8

m=result;

x=sqrt(x);

end

x

运行结果：

第二章 解线性方程组的直接解法

14. 考虑梯形电阻电路的设计，电路如下：

电路中的各个电流{1i ，2i ，…，8i

}须满足下列线性方程组：

R V i i =- 22 21

0 252321=-+-i i i 0 252 432=-+-i i i 0 252 543=-+-i i i 0 252 654=-+-i i i 0 252 765=-+-i i i 0 252 876=-+-i i i 052 87=+-i i

这是一个三对角方程组。设V 220=V ，Ω=27R ，运用追赶法，求各段电路的电流量。 程序设计如下：

a=[0 -2 -2 -2 -2 -2 -2 -2]; b=[2 5 5 5 5 5 5 5]; c=[-2 -2 -2 -2 -2 -2 -2]; d=[8.1481 0 0 0 0 0 0 0]; for i=2:8 a(i)=a(i)/b(i-1);

b(i)=b(i)-c(i-1)*a(i); d(i)=d(i)-a(i)*d(i-1); end;

d(8)=d(8)/b(8); for i=7:-1:1

d(i)=(d(i)-c(i)*d(i+1))/b(i); end;

x=d; x

运行结果：

第三章 解线性方程组的迭代法

14. 试分别用（1）Jacobi 迭代法；（2）Gauss-Seidel 迭代法；（3）共轭梯度法解线性方程组

??

?????

?

????????--=????????????????????????????????

--------12171427

12151

5

3

4

112323537123219143211054321x x x x x

迭代初始向量取)

0(x

=(

0，0，0，0，0T

)

。

解：（1）Jacobi 迭代法程序设计如下

clear clc i=1;

a=[10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15]; d=diag(diag(a)); l=d-tril(a); u=d-triu(a); d0=inv(d);

b=[12;-27;14;-17;12]; x0=zeros(5,1); B=d0*(l+u); f=d0*b; x=B*x0+f;

while norm(x-x0,inf)>=1e-8 x0=x; x=B*x0+f; i=i+1; end x i

运行结果：

（2）gauss-seidel迭代法程序设计如下：

clear

clc

i=1;

a=[10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15]; d=diag(diag(a));

l=d-tril(a);

u=d-triu(a);

b=[12;-27;14;-17;12];

x0=zeros(5,1);

B=inv(d-l)*u;

f=inv(d-l)*b;

x=B*x0+f;

while norm(x-x0,inf)>=1e-8

x0=x;

x=B*x0+f;

i=i+1;

end

x

i

运行结果：

（3）共轭梯度法程序设计如下：

A=[10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15];

b=[12;-27;14;-17;12];

x0=[0;0;0;0;0];

d0=b-A*x0;

r0=d0;

lamdr0=(r0*d0)/(d0*(A*d0));

x1=x0+lamdr0*d0;

k=0;Nmax=1000;

while abs(norm(x1-x0,inf))>1e-8 & k

r0=b-A*x1;

beta=-(r0*(A*d0))/(d0*(A*d0));

d0=r0+beta*d0;

lamdr0= r0*d0/(d0*(A*d0));

x0=x1;

x1=x1+lamdr0*d0;

k=k+1;

end;

x1

k

运行结果：

第五章插值法

已知直升机旋转机翼外形轮廓线上某些型值点的数据： K X k Y k K X k Y k

1 520 0

2 280 -30

3 156.6 -36

4 78 -35

5 39.62 -28.44

6 3.1 -9.4

7 0 0 8 3.1 9.4

9 39.62 28.44

10 78 35

11 156.6 36

12 280 30

13 520 0

由表中13个节点，用样条插值的方法，增加平面点为101个点，绘制比较光滑的机翼外形曲线。

x=[520,280,156.6,78,39.62,3.1,0];

x1=[0,3.1,39.62,78,156.6,280,520];

y=[0,-30,-36,-35,-28.44,-9.4,0];

y1=[0,9.4,28.44,35,36,30,0];

x2=520:-10.4:0;

y2=spline(x,y,x2);% 选用样条插值计算近似值

x3=0:10.4:520;

y3=spline(x1,y1,x3);

plot(x2,y2,'r-',x3,y3,'g-');

title('机翼外形曲线');

运行结果：

四、总结

通过本次的实验，基本掌握了这四种求解方程组的方法。并对matlab数学软件的作用有了初步的认识，数值计算这门课虽然本学期即将结束，但我相信这门课的思想和方法今后会伴随我们的研究，并且还有很多算法有待改进，有待我们去创新！