MATLAB实习报告-1

1、MATLAB实习报告班级:*姓名:*学号:*指导老师:*日期：2011/11/24目录.MATLAB简介上机实习题目第2章第第3章第第5章第第4章第.四.1题2题2题1题题目实现过程牛顿插值以及三次样条插值（第一个实验题）a. 牛顿插值b. 三次样条插值c. 牛顿插值截图d. 三次样条插值截图 3次，4次多项式的曲线拟合a. 先输入表格中对应的数据， b. 3次拟合c. 3次拟合截图d. 4次拟合e. 4次拟合截图高斯消去法解线性方程组 高斯消去法源代码第一个小题第二个小题雅克比迭代法与 SOF方法a. 雅克比迭代法的计算公式 b. 超松弛迭代法的计算公式 c. 得

2、到希尔矩阵d. 得到b矩阵e. 雅克比迭代法实现的函数 f. SOR迭代法实现的函数g. 对于雅克比迭代法，通过执行以下代码 .h. 对于SOR迭代法，执行相对应代码 .心得与感想a.b.c.99101111121313131414141415161619.MATLABS介MATLA是一个软件，用来进行科学计算。在实际生活或者在科学研究中， 大量的计算有时候是不可避免的，在这个时候，我们可以选择很多种方式来 解决我们的问题。但是选择一个好的软件对于我们来说有时候却是非常困难 的，尤其是现在各种软件层出不穷，当然里面有好的，可是有大多数的软件 却是鱼目混珠，在里面滥竽充数。如果我们不避开这些，那

3、么我们可能面临 的是用很大的工作量来完成一个极其简单的问题， 或者说在付出了时间与精 力之后却完不成问题。所以必须选择一个好的计算软件来解决问题。MATLAB就是这样一个软件，在你能够熟悉并且流畅运用它之后，那它将是一把能够解决很多问题的钥匙。MATLA是Matrix Laboratory ，即矩阵实验室的缩写。它是一个有 Math Work公司(由Moler，Little ，Bangert在1984年在加利福利亚成立)开发的 软件包。用来实现数值计算或者符号操作。另外MATLA还拥有M-file这一功能，在这里面你可以编写自己的程序(当 然这个程序不是windows程序)，在这里你就可以实现

4、很多计算，虽然说 MATLA自带的函数在大多数情况下已经够用，但是总会有某些时候你觉得自 己编写的程序运用起来更加灵活，而且在一些具体问题的时候，这个优势将 会是一大突破，可以让我们很轻松解决很多问题。 除了自己编写外，MATLAB 自带的很多优秀的程序，如plot (绘图)函数，co nd (求矩阵条件数)函数， 关于矩阵的一些相关函数，都是相当经典的。总之，熟悉使用MATLA可以让我们的工作变得很简单，经过这次的学习， 我们以后的很多学习工作问题都可以很方便地解决了。.上机实习题目1.第2章第1题已知函数在下列各点的值为Xi0.81.0f(Xi)0.980.920.810

5、.640.38试用4次牛顿插值多项式P4(x)及三次样条函数S(x)(自然边界条件)对数据 进行插值。用图给出(Xi,yi),Xi = 0.2+ 0.0 8i, i= 0,1,11,10，F4(x)及 S(x)。2.第3章第2题由实验给出数据表x0.0y1.00.410.500.610.912.022.46试用3次，4次多项式的曲线拟合，再根据数据曲线形状，求另外一个函数的 拟合曲线，用图示数据曲线及相应的三种拟合曲线。3.第5章第2题用列主元高斯消去法解线性方程组Ax = b。13.016.031.99X1(1) 1.274.16-1.23 X2=1;0.

6、987-4.819.34X313.006.031.99X111.274.16-1.23X2=1.0.990-4.819.34X31分别输出A, b, det (A),解向量x,中A勺条件数。分析比较(1),(2)。4.第4章第1题给出线性方程组Hnx= b，其中系数矩阵Hn为希尔伯特矩阵：Hn = (hij) 1Rnn,hij = T,i,j = 1,2,n。假设x? = (1, 1,1)T Rnn,b = HnX?。若取n = 6,8,10 ,分别用雅克比 迭代法及SO迭代法(3 = 1,1.25,1.5)求解。比较计算结果。.题目实现过程下面开始进行计算：1. 牛顿插值以及三次样条插值(第

7、一个实验题) 此题要求利用给定点，及给定点的函数值进行牛顿插值以及三次样条插 值。a. 牛顿插值要实现牛顿插值，要用到以下代码% 调用格式：yi=Lagran_(x,y,xi)% x,y数组形式的数据表function fi=Lagra n_(x,f,xi) fi=zeros(size(xi);np 1=le ngth(f);for i=1: np1z=on es(size(xi);for j=1: np1if i=j,z=z.*(xi-x(j)/(x(i)-x(j);e ndendfi=fi+z*f(i);endb. 三次样条插值要实现三次样条插值，要用到以下代码fun ctio n S=c

8、sfit(X,Y,dxO,dx n)N=le ngth(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-1temp=A(k-1)/B(k-1);B(k)=B(k)-tem p*C(k-1);U(k)=U(k)-tem p*U(k-1);endM(N)=U(N-1)/B(N-1);for k=N-2:-1:1

9、M(k+1)=(U(k)-C(k)*M(k+2)/B(k);endM(1)=3*(D(1)-dxO)/H(1)-M( 2)/2;M(N+1)=3*(dx n-D(N)/H(N)-M(N)/2;for k=0:N-1S(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);endc. 牛顿插值截图在输入对应的两个数组后，首先进行4次牛顿插值，下面是运行截图片i rr屋印是&!吁 nvni旳 mt H*LWWWfH時LgMitin

10、a-対打ZB"4筑1门11百日IhHU明日 di>ud-皤 dmlJ-i acEpXJI J 4f 4£144 4mbl4 ug * dtLI-i kpb* + dZO* bCT7“ dju41xle lit厂HIIEPVj：壬1 Ikii hi Cap J, SC)然后根据结果绘图d.三次样条插值截图 Uy 1 D! Lt；_；卞一"銮辆沁 s rr 斤：f H an-iiiMd wild*J 4kjurte cAm m事H逝旳mMiairiJJi4IIIrLJriLinL1庙Clue1*11dMtl4 dnMf 伽hl»z“ir ZZ 才vin

11、eP4 tjtHlh4lI. IICl -fl. HW 爼r*M血心De -a珀 I r 吨-TffTP山顽D于 EDU-rSUOdwH-iac EB7 I IL.*Cr3T fc E町 祚E ire石 II!*-polyvalO 函数首先得到的是三次样条插值后的多项式系数矩阵，然后利用计算在对应点的值,L-d44DlF 柯些 Eg 血a ft屋总A 'Airta .4" m .L 呼 "FH、hi7 CgiBiiid*trtiWinLwlbtal31'J1XITir-n-sr-S-I. IIbL 'fl.肛H 羽，t专H-2.0D1-irwLJUD

12、LL 爭®DDlJIOD>：kh = iHrKniii?. » xLJ山即ma-is>. rl jnhmLuMii. 和 iiF E h ?.aiic.:) il2 - lUM；- palTTwLSrX； .1 言-hJzhOhiiJi.I芥pld-别Ml.-叭-»l巧丄一4 冊 llDtfiLk rLp Hl占 a illj f '. J J | H4 ir*并且通过绘图函数将结果表示在图像上面,可以发现，利用三次样条插通过与上面的牛顿插值所得到的图像进行比较, 值得到的结果更加的精确。2. 3次，4次多项式的曲线拟合a.先输入表格中对应的数

13、据，然后首先通过调用plotO函数，绘图Lvlri-4 blh*.Urj L4IB rt S! £1 a；，左 M讥±cli m；14 F* :Jag卜 A U> V =(1 k1 C'. 3 Q.«1 1.41?.az ?.4C； plM毎- >Si_(W7 初尸MtKimi耐? Wlit311M1-53 -t.WWl icSVTLAETiwrTt * :| 也)i fc - 一 3 "tPfcilLMj J -bwb-AiUFile Edit 'Tievy Insert Tools Desktop Wind ow Help

14、fe您Q的®辰261.54-O.fi0.4O'. 50.60.7 OS 0日b.3次拟合通过以下代码>> x0 0.10000.20000.30000.50000.80001.0000>> y1.00002.46000.41000.50000.61000.91002.0200>> s1 = p olyfit(x, y, 3);>> temp = 0:0.001:1;>> y_temp = zeros(size(te mp);>> y_temp = pol yval(s1, tem p);

15、 >> hold on;>> pl ot(te mp, y_temp)c. 3次拟合截图可以得到下面的截图，这是看到3次拟合后的多项式与原来的图形比较File Edit View Ir&ert Tools Desktop Windtsvw Help怎殴氏膜 E B Bd. 4次拟合通过下面代码>> x0.100001.0000>> yy =1.00002.4600>> s2 = p olyfit(x, y, 4);>> y_temp = pol yval(s2, tem p);0.41000.20000.50000

16、.30000.61000.50000.91000.80002.0200>> hold on;>> plot(te mp, y_temp, 'r')e. 4次拟合截图通过进行4次拟合，我们会感觉不同，从图像上可以看出（不同次数的拟合 在图像上用不同的颜色标示出来以便观察， 其中3次拟合是蓝色，4次的是红色）以上都是用多项式进行拟合，下面代码用一个线性函数来进行拟合>> xmea n = mean( x);>> ymea n = mean( y);>> sumx2 = (x - xmea n) * (x - xmea n)

17、'>> sumxy = (y - ymea n) * (y - ymea n)'>> A = sumxy / sumx2;>> B = ymea n - A * xmea n;>> temp = 0:0.01:1;>> y_temp = A * temp + B;>> plot(te mp, y_temp, 'gree n')的到下图其中绿色的线是这次拟合的结果，很明显，误差很大BB Figure 1口回 ssFile Edit View Insert Tools Desktop Windo

18、vyr Help DQHSlfe題记目=酉在进行解决这个题目时候用的拟合函数是MATLAB自带的拟合函数polyfito ,这里是多项式拟合。polyfitO函数得到的是一个矩阵，然后运用另一个 MATLAB自带的函数polyvalO来进行计算给定点的函数值。到这里这一题完成。3. 高斯消去法解线性方程组对于这道题目，要求用高斯消去法解两个线性方程组，所以利用MATLAB的M-file编写一个文件便可以解决两个小题。而高斯消去法的核心是逐渐利用消去法是原来的系数矩阵变成一个上三 角矩阵。然后一个迭代便可以解出方程的每一个解。a.高斯消去法源代码fun ction X = p178_2_GAUS

19、S(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): nif (A(k, j) = 0)j) / A(i, j)；temp= A(k, :)+ A(i, :) * (-A(k,temp_b = b(k) + b(i) * (-A(k, j) / A(i, j); A(k, :) = temp;b(k) = temp_b;endendendi = i + 1;en d;Abdis p('det

20、(A) is .');x = det(A);dis p(x);dis p('con d(A) is .');x = con d(A);dis p(x);X(n) = b( n) / A( n, n);for i = (n - 1):-1:1temp_b = 0;for j = (i + 1): ntemp_b = temp_b + A(i, j) * X(j);endX(i) = (b(i) - temp_b) / A(i, i);endendb。从这个程序要求输入两个参数，一个是系数矩阵，另外一个是矩阵程序中我们可以看出程序首先初始化一个解向量，全部以0为元素。然后通

21、过循环来进行系数矩阵的上三角化，在这个过程中，需要判断每一个系数矩阵的行向量的第一元素是否为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;>> P 178_2_GAUSS(A, b)可以解出方程：A =3.0100 6.03001.990001.6158 -2.069600 -0.0063b =1.0000

22、0.57813.1004 det(A) is .-0.0305 con d(A) is .4.0724e+003ans =1.0e+003 *1.5926-0.6319-0.4936c. 第二个小题通过以下过程：>> A = 3.00 6.03 1.99;1.27 4.16 -1.23; 0.990 -4.81 9.34;>> b = 1; 1; 1;>> p 178_2_GAUSS(A, b)可以解出方程：A =3.0000 6.03001.990001.6073 -2.072400 -0.0844b =1.00000.57673.1097det(A) i

23、s .-0.4070con d(A) is .305.2933ans =119.5273-47.1426-36.8403到这里这一大题已经完全解决，下面进入下一题。4. 雅克比迭代法与SO方法这一题考查的是迭代法的使用，雅克比迭代法与SO方法，在这道题中用。a.雅克比迭代法的计算公式x(k+1 )X=(x(0),x20), £，ni = 1,2,n;k =b.超松弛迭代法的计算公式(0)Z(0)(0)X 丿=(X1 ,X2，3 (b - Zj=i aijX(k)(k+1Xij=1j Mi0,1,表示迭代次数,xn0)T,=x(k) + 3 (bi -耳i aijX(k)- Zj=11

24、 aijX(k+1)?aii ,i = 1,2,n; k = 0,1,表示迭代次数3为松弛因子首先要通过程序获得希尔矩阵与b矩阵，然后才能进行下面的工作，而要 实现这个功能需要以下的函数：=(bi-刀 aijX(k)? aii ,c.得到希尔矩阵fun ctio n Hn = GET_Hn(n) Hn = zeros( n);for i = 1:nfor j = 1:nHn(i, j) = 1 / (i + j - 1); endendendd.得到b矩阵fun ctio n b = GET_b( n) Hn_temp = GET_ Hn(n);b = zeros( n, 1);for i =

25、 1:nfor j = 1:n b(i) = b(i) + Hn_te mp(i, j);endendende. 雅克比迭代法实现的函数fun ction X = p211_1_JJ( n)Hn = GET_H n(n);b = GET_b( n);temp = 0;X0 = zeros(1, n);X_old = zeros(1, n);X_new = zeros(1, n);dis p('Now Jacobi method!');dis p('Start with the vector that (0, 0, 0, .FT');for i = 1:nfor

26、k = 1:nX old = X new;temp = 0;for j = 1:n if(j = i)temp = temp + Hn(i, j) * X_old(j);endendX_new(i) = (b(i) - tem p) / Hn (i, i);endendX = X_new;endf. SOI迭代法实现的函数fun ction X = p211_1_SOR( n, w)Hn = GET_H n(n);b = GET_b( n);tempOl = 0;temp02 = 0;XO = zeros(1, n);X_old = zeros(1, n);X_new = zeros(1, n

27、);dis p('Now Successive Over Relaxtio n method!');dis p('Start with the vector that (0, 0, 0, .)5);for i = 1:nfor k = 1:nX_old = X_new;temp01 = 0;temp02 = 0;for j = 1:nif(j < i)temp01 = temp01 + Hn (i, j) * X_new(j);endif(j > i)temp02 = temp02 + Hn (i, j) * X_old(j);endendendX_new(

28、i) = w * (b(i)- temp01- temp02)/ Hn(i, i) + X_old(i);endX = X_new;End对于雅克比迭代法，通过执行以下代码P211_1_JJ (6)p211_1_JJ(8) p211_1_JJ(10)g.>>>>>>可以分别得到：Now Jacobi method!Start with the vector that (0, 0, 0, .FTans =2.45001.10360.62650.40600.28310.2071Now Jacobi method!Start with the vector tha

29、t (0, 0, 0, .FT ans =2.71791.41010.85240.58090.42210.31980.24970.1995Now Jacobi method!Start with the vector that (0, 0, 0, .FTans =Colu mns 1 through 92.92901.66621.05170.74230.55540.43150.34450.28070.2325Colu mn 100.1951h.对于SO迭代法，执行相对应代码n=6,沪1,1.25,1.5 的时候>> P 211_1_SOR(6, 1)Now Successive O

30、ver Relaxti on method!Start with the vector that (0, 0, 0, .FT ans =2.45001.10360.62650.40600.28310.2071>> p 211_1_SOR(6, 1.25)Now Successive Over Relaxti on method!Start with the vector that (0, 0, 0, .FTans =3.06250.23100.87040.33890.31410.2097>> p 211_1_SOR(6, 1.5)Now Successive Over

31、Relaxti on method!Start with the vector that (0, 0, 0, .FTans =3.6750 -1.1009 2.0106 -0.39940.7670 -0.0384与 n=8, 0=1,1.25,1.5 的时候>> p 211_1_SOR(8, 1)Now Successive Over Relaxti on method!Start with the vector that (0, 0, 0, .FTans =2.71791.41010.85240.58090.24970.19950.42210.3198>> p 211

32、_1_SOR(8, 1.25)Now Successive Over Relaxti on method!Start with the vector that (0, 0, 0, .FTans =3.39730.48871.08980.50620.25730.20420.45010.3203>> p 211_1_SOR(8, 1.5)Now Successive Over Relaxti on method!Start with the vector that (0, 0, 0, .FTans =4.0768 -0.9424 2.2923 -0.27530.40710.1275与n=10,沪1,1.25,1.5 的时候>> p 211_1_SOR(10, 1)0.92520.0578Now Successive Over Relaxti on method!