MATLAB实习报告-1

MATLAB --实习报告

班级：******

姓名：**

学号：***********

指导老师：***

日期：2011/11/24

目录

一．MATLAB简介 (3)

二．上机实习题目 (3)

1.第2章第1题 (3)

2.第3章第2题 (3)

3.第5章第2题 (3)

4.第4章第1题 (4)

三．题目实现过程 (4)

1.牛顿插值以及三次样条插值（第一个实验题） (4)

a.牛顿插值 (4)

b.三次样条插值 (4)

c.牛顿插值截图 (5)

d.三次样条插值截图 (6)

2.3次，4次多项式的曲线拟合 (7)

a.先输入表格中对应的数据， (7)

b.3次拟合 (8)

c. 3次拟合截图 (9)

d. 4次拟合 (9)

e. 4次拟合截图 (10)

3.高斯消去法解线性方程组 (11)

a.高斯消去法源代码 (11)

b.第一个小题 (12)

c.第二个小题 (13)

4.雅克比迭代法与SOR方法 (13)

a. 雅克比迭代法的计算公式 (13)

b. 超松弛迭代法的计算公式 (14)

c.得到希尔矩阵 (14)

d.得到b矩阵 (14)

e.雅克比迭代法实现的函数 (14)

f.SOR迭代法实现的函数 (15)

g.对于雅克比迭代法，通过执行以下代码 (16)

h. 对于SOR迭代法，执行相对应代码 (16)

四．心得与感想 (19)

一．M ATLAB简介

MATLAB是一个软件，用来进行科学计算。在实际生活或者在科学研究中，大量的计算有时候是不可避免的，在这个时候，我们可以选择很多种方式来解决我们的问题。但是选择一个好的软件对于我们来说有时候却是非常困难的，尤其是现在各种软件层出不穷，当然里面有好的，可是有大多数的软件却是鱼目混珠，在里面滥竽充数。如果我们不避开这些，那么我们可能面临的是用很大的工作量来完成一个极其简单的问题，或者说在付出了时间与精力之后却完不成问题。所以必须选择一个好的计算软件来解决问题。MATLAB 就是这样一个软件，在你能够熟悉并且流畅运用它之后，那它将是一把能够解决很多问题的钥匙。

MATLAB是Matrix Laboratory，即矩阵实验室的缩写。它是一个有Math Work公司（由Moler，Little，Bangert在1984年在加利福利亚成立）开发的软件包。用来实现数值计算或者符号操作。

另外MATLAB还拥有M-file这一功能，在这里面你可以编写自己的程序（当然这个程序不是windows程序），在这里你就可以实现很多计算，虽然说

MATLAB自带的函数在大多数情况下已经够用，但是总会有某些时候你觉得自己编写的程序运用起来更加灵活，而且在一些具体问题的时候，这个优势将会是一大突破，可以让我们很轻松解决很多问题。除了自己编写外，MATLAB 自带的很多优秀的程序，如plot（绘图）函数，cond（求矩阵条件数）函数，关于矩阵的一些相关函数，都是相当经典的。

总之，熟悉使用MATLAB可以让我们的工作变得很简单，经过这次的学习，我们以后的很多学习工作问题都可以很方便地解决了。

二．上机实习题目

1.第2章第1题

4

进行插值。用图给出{(x i,y i),x i=0.2+0.08i,i=0,1,11,10}，P4(x)及S(x)。

2.第3章第2题

拟合曲线，用图示数据曲线及相应的三种拟合曲线。

3.第5章第2题

用列主元高斯消去法解线性方程组Ax=b。

(1)[3.01 6.03 1.99

1.27 4.16?1.23

0.987?4.819.34

][

x1

x2

x3

]=[

1

1

1

];

(2) [3.00 6.03 1.99

1.27 4.16?1.23

0.990?4.819.34

][

x1

x2

x3

]=[

1

1

1

].

分别输出A，b，det (A),解向量x，(1)中A的条件数。分析比较(1),(2)。

4.第4章第1题

给出线性方程组H n x=b，其中系数矩阵H n为希尔伯特矩阵：H n=(h ij)∈

,i,j=1,2,…,n。

R n×n,h ij=1

i+j?1

假设x?=(1,1,…,1)T∈R n×n,b=H n x?。若取n=6,8,10，分别用雅克比

迭代法及SOR迭代法(ω=1,1.25,1.5)求解。比较计算结果。

三．题目实现过程

下面开始进行计算：

1.牛顿插值以及三次样条插值（第一个实验题）

此题要求利用给定点，及给定点的函数值进行牛顿插值以及三次样条插

值。

a.牛顿插值

要实现牛顿插值，要用到以下代码

% 调用格式：yi=Lagran_(x,y,xi)

% x,y 数组形式的数据表

function fi=Lagran_(x,f,xi)

fi=zeros(size(xi));

np1=length(f);

for i=1:np1

z=ones(size(xi));

for j=1:np1

if i~=j,z=z.*(xi-x(j))/(x(i)-x(j));end

end

fi=fi+z*f(i);

end

b.三次样条插值

要实现三次样条插值，要用到以下代码

function S=csfit(X,Y,dx0,dxn)

N=length(X)-1;

H=diff(X);

D=diff(Y)./H;

A=H(2:N-1);

B=2*(H(1:N-1)+H(2:N));

C=H(2:N);

U=6*diff(D);

B(1)=B(1)-H(1)/2;

U(1)=U(1)-3*(D(1));

B(N-1)=B(N-1)-H(N)/2;

U(N-1)=U(N-1)-3*(-D(N));

for k=2:N-1

temp=A(k-1)/B(k-1);

B(k)=B(k)-temp*C(k-1);

U(k)=U(k)-temp*U(k-1);

end

M(N)=U(N-1)/B(N-1);

for k=N-2:-1:1

M(k+1)=(U(k)-C(k)*M(k+2))/B(k);

end

M(1)=3*(D(1)-dx0)/H(1)-M(2)/2;

M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;

for k=0:N-1

S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));

S(k+1,2)=M(k+1)/2;

S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;

S(k+1,4)=Y(k+1);

end

c.牛顿插值截图

在输入对应的两个数组后，首先进行4次牛顿插值，下面是运行截图

然后根据结果绘图

d.三次样条插值截图

首先得到的是三次样条插值后的多项式系数矩阵，然后利用polyval()函数计算在对应点的值，

并且通过绘图函数将结果表示在图像上面，

通过与上面的牛顿插值所得到的图像进行比较，可以发现，利用三次样条插值得到的结果更加的精确。

2.3次，4次多项式的曲线拟合

a.先输入表格中对应的数据，

然后首先通过调用plot()函数，绘图

b.3次拟合

通过以下代码

>> x

x =

0 0.1000 0.2000 0.3000 0.5000 0.8000 1.0000

>> y

y =

1.0000 0.4100 0.5000 0.6100 0.9100

2.0200

2.4600

>> s1 = polyfit(x, y, 3);

>> temp = 0:0.001:1;

>> y_temp = zeros(size(temp));

>> y_temp = polyval(s1, temp);

>> hold on;

>> plot(temp, y_temp)

c. 3次拟合截图

可以得到下面的截图，这是看到3次拟合后的多项式与原来的图形比较

d. 4次拟合

通过下面代码

>> x

x =

0 0.1000 0.2000 0.3000 0.5000 0.8000

1.0000

>> y

y =

1.0000 0.4100 0.5000 0.6100 0.9100

2.0200

2.4600

>> s2 = polyfit(x, y, 4);

>> y_temp = polyval(s2, temp);

>> hold on;

>> plot(temp, y_temp, 'r')

e. 4次拟合截图

通过进行4次拟合，我们会感觉不同，从图像上可以看出（不同次数的拟合在图像上用不同的颜色标示出来以便观察，其中3次拟合是蓝色，4次的是红色）

以上都是用多项式进行拟合，下面代码用一个线性函数来进行拟合

>> xmean = mean(x);

>> ymean = mean(y);

>> sumx2 = (x - xmean) * (x - xmean)';

>> sumxy = (y - ymean) * (y - ymean)';

>> A = sumxy / sumx2;

>> B = ymean - A * xmean;

>> temp = 0:0.01:1;

>> y_temp = A * temp + B;

>> plot(temp, y_temp, 'green')

的到下图其中绿色的线是这次拟合的结果，很明显，误差很大

在进行解决这个题目时候用的拟合函数是MATLAB自带的拟合函数polyfit()，这里是多项式拟合。polyfit()函数得到的是一个矩阵，然后运用另一个MATLAB自带的函数polyval()来进行计算给定点的函数值。

到这里这一题完成。

3.高斯消去法解线性方程组

对于这道题目，要求用高斯消去法解两个线性方程组，所以利用MATLAB 的M-file编写一个文件便可以解决两个小题。

而高斯消去法的核心是逐渐利用消去法是原来的系数矩阵变成一个上三角矩阵。然后一个迭代便可以解出方程的每一个解。

a.高斯消去法源代码

function X = p178_2_GAUSS(A, b)

[n, m] = size(A);

X = zeros(n, 1);

temp = zeros(1, m);

temp_b = 0;

i = 1;

for j = 1: (m - 1)

if (A(i, j) ~= 0)

for k = (i + 1):n

if (A(k, j) ~= 0)

temp = A(k, :) + A(i, :) * (-A(k, j) / A(i, j));

temp_b = b(k) + b(i) * (-A(k, j) / A(i, j));

A(k, :) = temp;

b(k) = temp_b;

end

end

end

i = i + 1;

end;

A

b

disp('det(A) is ...');

x = det(A);

disp(x);

disp('cond(A) is ...');

x = cond(A);

disp(x);

X(n) = b(n) / A(n, n);

for i = (n - 1):-1:1

temp_b = 0;

for j = (i + 1):n

temp_b = temp_b + A(i, j) * X(j);

end

X(i) = (b(i) - temp_b) / A(i, i);

end

end

这个程序要求输入两个参数，一个是系数矩阵，另外一个是矩阵b。从程序中我们可以看出程序首先初始化一个解向量，全部以0为元素。然后通过循环来进行系数矩阵的上三角化，在这个过程中，需要判断每一个系数矩阵的行向量的第一元素是否为0，若是0的话，则需要进行换行操作，或者跳过此次循环，直接进入下一次循环（因为有可能是这一行已经达到要求，不需要再进行操作。）

b.第一个小题

通过以下过程：

>> A = [3.01 6.03 1.99;1.27 4.16 -1.23; 0.987 -4.81 9.34];

>> b = [1; 1; 1];

>> p178_2_GAUSS(A, b)

可以解出方程：

A =

3.0100 6.0300 1.9900

0 1.6158 -2.0696

0 0 -0.0063

b =

1.0000

0.5781

3.1004

det(A) is ...

-0.0305

cond(A) is ...

4.0724e+003

ans =

1.0e+003 *

1.5926

-0.6319

-0.4936

c.第二个小题

通过以下过程：

>> A = [3.00 6.03 1.99;1.27 4.16 -1.23; 0.990 -4.81 9.34];

>> b = [1; 1; 1];

>> p178_2_GAUSS(A, b)

可以解出方程：

A =

3.0000 6.0300 1.9900

0 1.6073 -2.0724

0 0 -0.0844

b =

1.0000

0.5767

3.1097

det(A) is ...

-0.4070

cond(A) is ...

305.2933

ans =

119.5273

-47.1426

-36.8403

到这里这一大题已经完全解决，下面进入下一题。

4.雅克比迭代法与SOR方法

这一题考查的是迭代法的使用，雅克比迭代法与SOR方法，在这道题中用。

a. 雅克比迭代法的计算公式

{

x (0)=(x 1(0),x 2(0),…,x n (0)

)T ,

x i (k+1)=(b i ?∑a ij x i (k )

n

j=1j≠i

)a ii ?,

i =1,2,…n;k =0,1,…表示迭代次数

b. 超松弛迭代法的计算公式

{

x (0)=(x 1(0),x 2(0),…,x n (0)

)T ,

x i (k+1)=x i

(k )+ω((b i ?∑a ij x i (k )n j=i ?∑a ij x i (k+1)i?1j=1))a ii ?,

i =1,2,…n;k =0,1,…表示迭代次数

ω为松弛因子

首先要通过程序获得希尔矩阵与b 矩阵，然后才能进行下面的工作，而要实现这个功能需要以下的函数： c.得到希尔矩阵

function Hn = GET_Hn(n) Hn = zeros(n); for i = 1:n for j = 1:n

Hn(i, j) = 1 / (i + j - 1); end end end d.得到b 矩阵

function b = GET_b(n) Hn_temp = GET_Hn(n); b = zeros(n, 1); for i = 1:n for j = 1:n

b(i) = b(i) + Hn_temp(i, j); end end end

e. 雅克比迭代法实现的函数

function X = p211_1_JJ(n) Hn = GET_Hn(n); b = GET_b(n); temp = 0;

X0 = zeros(1, n);

X_old = zeros(1, n);

X_new = zeros(1, n);

disp('Now Jacobi method!');

disp('Start with the vector that (0, 0, 0, ...)^T');

for i = 1:n

for k = 1:n

X_old = X_new;

temp = 0;

for j = 1:n

if(j ~= i)

temp = temp + Hn(i, j) * X_old(j);

end

end

X_new(i) = (b(i) - temp) / Hn(i, i);

end

end

X = X_new;

end

f.SOR迭代法实现的函数

function X = p211_1_SOR(n, w)

Hn = GET_Hn(n);

b = GET_b(n);

temp01 = 0;

temp02 = 0;

X0 = zeros(1, n);

X_old = zeros(1, n);

X_new = zeros(1, n);

disp('Now Successive Over Relaxtion method!');

disp('Start with the vector that (0, 0, 0, ...)^T');

for i = 1:n

for k = 1:n

X_old = X_new;

temp01 = 0;

temp02 = 0;

for j = 1:n

if(j < i)

temp01 = temp01 + Hn(i, j) * X_new(j);

end

if(j > i)

temp02 = temp02 + Hn(i, j) * X_old(j);

end

end

end

X_new(i) = w * (b(i) - temp01 - temp02) / Hn(i, i) + X_old(i);

end

X = X_new;

End

g.对于雅克比迭代法，通过执行以下代码

>> p211_1_JJ(6)

>> p211_1_JJ(8)

>> p211_1_JJ(10)

可以分别得到：

Now Jacobi method!

Start with the vector that (0, 0, 0, ...)^T

ans =

2.4500 1.1036 0.6265 0.4060 0.2831 0.2071 Now Jacobi method!

Start with the vector that (0, 0, 0, ...)^T

ans =

2.7179 1.4101 0.8524 0.5809 0.4221 0.3198

0.2497 0.1995

Now Jacobi method!

Start with the vector that (0, 0, 0, ...)^T

ans =

Columns 1 through 9

2.9290 1.6662 1.0517 0.7423 0.5554 0.4315 0.3445 0.2807 0.2325

Column 10

0.1951

h. 对于SOR迭代法，执行相对应代码

n=6, ω=1,1.25,1.5的时候

>> p211_1_SOR(6, 1)

Now Successive Over Relaxtion method!

Start with the vector that (0, 0, 0, ...)^T

ans =

2.4500 1.1036 0.6265 0.4060 0.2831 0.2071

>> p211_1_SOR(6, 1.25)

Now Successive Over Relaxtion method!

Start with the vector that (0, 0, 0, ...)^T

ans =

3.0625 0.2310 0.8704 0.3389 0.3141 0.2097

>> p211_1_SOR(6, 1.5)

Now Successive Over Relaxtion method!

Start with the vector that (0, 0, 0, ...)^T

ans =

3.6750 -1.1009 2.0106 -0.3994 0.7670 -0.0384

与n=8, ω=1,1.25,1.5的时候

>> p211_1_SOR(8, 1)

Now Successive Over Relaxtion method!

Start with the vector that (0, 0, 0, ...)^T

ans =

2.7179 1.4101 0.8524 0.5809 0.4221 0.3198 0.2497 0.1995

>> p211_1_SOR(8, 1.25)

Now Successive Over Relaxtion method!

Start with the vector that (0, 0, 0, ...)^T

ans =

3.3973 0.4887 1.0898 0.5062 0.4501 0.3203 0.2573 0.2042

>> p211_1_SOR(8, 1.5)

Now Successive Over Relaxtion method!

Start with the vector that (0, 0, 0, ...)^T

ans =

4.0768 -0.9424 2.2923 -0.2753 0.9252 0.0578 0.4071 0.1275

与n=10, ω=1,1.25,1.5的时候

>> p211_1_SOR(10, 1)

Now Successive Over Relaxtion method!

Start with the vector that (0, 0, 0, ...)^T

ans =

Columns 1 through 9

2.9290 1.6662 1.0517 0.7423 0.5554 0.4315 0.3445 0.2807 0.2325

Column 10

0.1951

>> p211_1_SOR(10, 1.25)

Now Successive Over Relaxtion method!

Start with the vector that (0, 0, 0, ...)^T

ans =

Columns 1 through 9

3.6612 0.7098 1.2835 0.6617 0.5807 0.4299 0.3506 0.2844 0.2363

Column 10

0.1984

>> p211_1_SOR(10, 1.5)

Now Successive Over Relaxtion method!

Start with the vector that (0, 0, 0, ...)^T

ans =

Columns 1 through 9

4.3935 -0.7958 2.5326 -0.1523 1.0720 0.1565 0.5050 0.2041 0.2819

Column 10

0.1766

到此这道题目结束了。

四．心得与感想

在此次实验的四个题目中，第一个题使用的是老师的程序代码，第二个题是用的是MATLAB自带的函数，第三个题和第四个题中的5个函数（包括p178_2_GAUSS(A, b),GET_Hn(n)，GET_b(n)，

p211_1_JJ(A,b),p211_1_SOR(A,b)）是自己写的程序代码。

在学习数值分析的同时，进行学习MATLAB，有助于将两者都更加的熟悉并加以充分利用。比如在学习数值分析过程中对某个问题不是很明白，也许我们就会找出各种借口来将它跳过去，以至于过了一段时间之后连自己有什么问题都忘了，就给以后的生活留下了一个炸弹，有可能在某天我们急需要解决一个问题，而要解决那个问题就需要在数值分析中有个一个知识点，但是因为在学的时候我们没有搞懂，所以在那时候就不能够及时做出来，有可能造成其他问题。但是由于要进行MATLAB实验，我们就必须要把那个问题解决，要不然就不能够让MATLAB运行。同理对于MATLAB的学习是一样的。两者促进另外一者的学习。

MATLAB是一个软件，更可以是一门语言，还有好多的功能我们未曾用过，那些都是值得去搞懂的，但是由于时间原因，不可能停下学习，所以现在有这么多，先满足吧。