MATLAB实习报告-1.doc

MATLAB-实习报告班级：*姓名：*学号：*指导老师：*日期：2011/11/24目录一MATLAB简介3二上机实习题目31.第2章第1题32.第3章第2题33.第5章第2题34.第4章第1题4三题目实现过程41.牛顿插值以及三次样条插值（第一个实验题）4a.牛顿插值4b.三次样条插值4c.牛顿插值截图5d.三次样条插值截图62.3次，4次多项式的曲线拟合7a.先输入表格中对应的数据， 7b.3次拟合7c. 3次拟合截图8d. 4次拟合9e. 4次拟合截图93.高斯消去法解线性方程组10a.高斯消去法源代码10b.第一个小题11c.第二个小题124.雅克比迭代法与SOR方法13a. 雅克比迭代法的计算公式13b. 超松弛迭代法的计算公式13c.得到希尔矩阵13d.得到b矩阵13e.雅克比迭代法实现的函数13f.SOR迭代法实现的函数14g.对于雅克比迭代法，通过执行以下代码15h. 对于SOR迭代法，执行相对应代码15四心得与感想18一 MATLAB简介MATLAB是一个软件，用来进行科学计算。在实际生活或者在科学研究中，大量的计算有时候是不可避免的，在这个时候，我们可以选择很多种方式来解决我们的问题。但是选择一个好的软件对于我们来说有时候却是非常困难的，尤其是现在各种软件层出不穷，当然里面有好的，可是有大多数的软件却是鱼目混珠，在里面滥竽充数。如果我们不避开这些，那么我们可能面临的是用很大的工作量来完成一个极其简单的问题，或者说在付出了时间与精力之后却完不成问题。所以必须选择一个好的计算软件来解决问题。MATLAB就是这样一个软件，在你能够熟悉并且流畅运用它之后，那它将是一把能够解决很多问题的钥匙。MATLAB是Matrix Laboratory，即矩阵实验室的缩写。它是一个有Math Work公司（由Moler，Little，Bangert在1984年在加利福利亚成立）开发的软件包。用来实现数值计算或者符号操作。另外MATLAB还拥有M-file这一功能，在这里面你可以编写自己的程序（当然这个程序不是windows程序），在这里你就可以实现很多计算，虽然说MATLAB自带的函数在大多数情况下已经够用，但是总会有某些时候你觉得自己编写的程序运用起来更加灵活，而且在一些具体问题的时候，这个优势将会是一大突破，可以让我们很轻松解决很多问题。除了自己编写外，MATLAB自带的很多优秀的程序，如plot（绘图）函数，cond（求矩阵条件数）函数，关于矩阵的一些相关函数，都是相当经典的。总之，熟悉使用MATLAB可以让我们的工作变得很简单，经过这次的学习，我们以后的很多学习工作问题都可以很方便地解决了。二 上机实习题目1. 第2章第1题已知函数在下列各点的值为xi0.20.40.60.81.0f(xi)0.980.920.810.640.38试用4次牛顿插值多项式P4(x)及三次样条函数S(x)（自然边界条件）对数据进行插值。用图给出xi, yi, xi=0.2+0.08i, i=0, 1, 11, 10，P4(x)及S(x)。2. 第3章第2题由实验给出数据表x0.00.10.20.30.50.81.0y1.00.410.500.610.912.022.46试用3次，4次多项式的曲线拟合，再根据数据曲线形状，求另外一个函数的拟合曲线，用图示数据曲线及相应的三种拟合曲线。3. 第5章第2题用列主元高斯消去法解线性方程组Ax=b。(1)3.016.031.991.274.16-1.230.987-4.819.34x1x2x3=111;(2) 3.006.031.991.274.16-1.230.990-4.819.34x1x2x3=111.分别输出A，b，det(A),解向量x，(1)中A的条件数。分析比较(1),(2)。 4. 第4章第1题给出线性方程组Hnx=b，其中系数矩阵Hn为希尔伯特矩阵：Hn=hijRnn, hij=1i+j-1, i,j=1, 2, n。假设x*=(1, 1, , 1)TRnn, b=Hnx*。若取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);endb.三次样条插值要实现三次样条插值，要用到以下代码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);endM(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次拟合通过以下代码 xx =0 0.1000 0.2000 0.3000 0.5000 0.8000 1.0000 yy =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次拟合通过下面代码 xx =0 0.1000 0.2000 0.3000 0.5000 0.8000 1.0000 yy =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;Abdisp(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);endend这个程序要求输入两个参数，一个是系数矩阵，另外一个是矩阵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.0063b = 1.0000 0.5781 3.1004det(A) is . -0.0305cond(A) is . 4.0724e+003ans = 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.0844b = 1.0000 0.5767 3.1097det(A) is . -0.4070cond(A) is . 305.2933ans = 119.5273 -47.1426 -36.8403 到这里这一大题已经完全解决，下面进入下一题。 4. 雅克比迭代法与SOR方法这一题考查的是迭代法的使用，雅克比迭代法与SOR方法，在这道题中用。a. 雅克比迭代法的计算公式x(0)=(x10, x20, , xn(0)T,xik+1=(bi-j=1jinaijxik)aii,i=1,2,n;k=0,1,表示迭代次数b. 超松弛迭代法的计算公式x(0)=(x10, x20, , xn(0)T,xik+1=xik+(bi-j=inaijxik-j=1i-1aijxik+1)aii,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); endendendd.得到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); endendende. 雅克比迭代法实现的函数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); endendX = X_new;endf. 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) 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);endX = X_new;Endg. 对于雅克比迭代法，通过执行以下代码 p211_1_JJ(6) p211_1_JJ(8) p211_1_JJ(10)可以分别得到：Now Jacobi method!Start with the vector that (0, 0, 0, .)Tans = 2.4500 1.1036 0.6265 0.4060 0.2831 0.2071Now Jacobi method!Start with the vector that (0, 0, 0, .)Tans =2.7179 1.4101 0.8524 0.5809 0.4221 0.3198 0.2497 0.1995Now Jacobi method!Start with the vector that (0, 0, 0, .)Tans = 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.1951h. 对于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, .)Tans = 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, .)Tans = 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, .)Tans = 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, .)Tans =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, .)Tans =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, .)Tans =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, .)Tans = 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