数值分析上机实验报告.docx

数值分析上机报告学 院： 专 业： 学 号： 姓 名： 指导老师： 联系电话： 目录1 序 言31.1 MATLAB语言特点31.2 选用原因42 第一题52.1 题目52.2 计算结果分析63 第二题93.1 题目93.2 结果分析104 第三题114.1 题目114.2 计算结果分析12附录一 多项式拟合程序代码15附录二 多项式插值程序代码16附录三 jacobi迭代程序代码17附录四 Gausssidel 迭代程序代码18附录五 四阶Runge-Kutta算法程序代码19总结与体会201 序 言 1.1 MATLAB语言特点近年来MATLAB之所以能在各个行业得到广泛的应用，主要得益于MATLAB语言与其他编程语言相比不同的特点：高效的矩阵运算机制；多样化的操作途径；功能强大的工具箱；良好的扩充能力；完善的帮助系统。1高效的矩阵运算机制MATLAB软件是基于矩阵计算开发的，在其他编程语言中需要使用多个for语句才能完成的操作，在MATLAB中直接使用矩阵即可完成计算，因而MATLAB在数据计算分析，特别是对海量数据的处理方面表现出相比其他的编程语言更大的优势。2多样化的操作途径MATLAB语言为用户提供了多种操作方式选择。用户可以编写代码实现各种功能，代码可重复利用，同时，不擅长编程的用户也可以通过MATLAB图形界面操作，完成MATLAB的相应功能。3功能强大的工具箱对于算法的开发，MATLAB提供了大量现成的函数，用户可以直接调用。MATLAB软件对目前主流的算法都提供了现成的函数，并封装成一个个适用于不同领域的算法工具箱。常用的工具箱包括数学类、经济类、信号处理类工具箱，同时随着MATLAB软件版本的更新，不断有功能更强大的工具箱添加到MATLAB中。MATLAB语言的这一特点，将大大节省算法开发的时间，用户无须详细了解算法的细节，只要掌握MATLAB中算法相应函数的调用即可。4良好的扩展能力利用MATLAB语言编写的程序具有良好的扩展能力，可以方便地与各种编程语言链接。用户可以方便地在MATLAB中调用其他语言已编写好的程序，同时在其他语言中也可以方便地调用MATLAB的程序。MATLAB语言具有良好的接口编程技术。5完善的帮助系统完善的帮助系统是MATLAB的又一突出特点，MATLAB向用户提供了多种帮助途径。通过MATLAB的帮助系统，用户可以获取MATLAB常用函数的使用方法及应用实例，而且这种帮助可以是实时的、在线的。同时，为了便于用户更好地使用MATLAB软件，在MATLAB中的主要算法都是可以直接看到源代码的。1.2 选用原因本次数值分析上机实习采用MATLAB语言进行编程，主要是基于以下几方面的考虑：1编程效率高MATLAB语言是一种面向科学与工程计算的高级语言，允许用数学形式的语言编写程序，且比BASIC、FORTRAN和C等语言更加接近我们书写计算公式的思维方式，用MATLAB编写程序犹如在演算纸上排列公式与求解问题。因此，也可通俗地称MATLAB语言为演算纸式科学算法语言。由于它编程简单，所以编程效率高，易学易懂。2高效方便的矩阵和数组运算MATLAB语言像BASIC、FORTRAN和C语言一样规定了矩阵的算术运算符、关系运算符、逻辑运算符、条件运算符及赋值运算符，而且这些运算符大部分可以毫无改变的照搬到数组建的运算中。有些如算术运算符只要增加“.”就可以用于数组建间的运算。另外，它不需定义数组的维数，并给出矩阵函数、特殊矩阵专门的库函数，使之在求解数字图像处理问题时显得大为简洁、高效、方便，这是其他高级语言所不能相比的。3用户使用方便MATLAB语言是一种解释执行的语言(在没有被专门的工具编译之前)，它灵活、方便，起跳时程序手段丰富，调试速度快，需要学习时间少。人们用任何一种语言编程和调试一般都要经过四个步骤：编辑、编译、连接，以及执行和调试。各个步骤之间是顺序关系，编程的过程就是在他们之间作瀑布型的循环。MATLAB语言与其它语言相比，较好的解决了上述问题，把编辑、编译、连接和执行融为一体。它能在同一画面上进行灵活操作，快速排除输入程序的书写错误、语法错误甚至语义错误，从而加快了用户编写、修改和调试程序的速度，可以说在编程和调试过程中它是一种比VB还要简单的语言，更不用谈C和C+。4扩充性强，交互性好MATLAB语言有丰富的库函数，在进行复杂的数学运算时可以直接调用，且这些库函数同用户文件在形式上一样，所以用户文件也可以作为MATLAB的库函数来调用。婴儿，用户可以根据自己的需要方便的建立和扩充新的库函数，提高MATLAB使用效率和扩充它的功能。另外，为了充分利用FORTRAN、C等语言的资源，包括用户自己编好的FORTRAN、C语言程序，通过建立M文件的形式，混合编程，方便地调用有关的FORTRAN、C语言子程序，还可以在C语言和FORTRAN语言中方便的使用MATLAB的数值计算功能。良好的交互性增加了代码的重用性。综上所述，MATLAB最适用的范围是科学计算，对于本实习需要处理矩阵运算来实现图像处理算法来说，MATLAB再合适不过了。2 第一题2.1 题目某过程测涉及两变量x 和y, 拟分别用插值多项式和多项式拟合给出其对应规律的近似多项式，已知xi与yi之间的对应数据如下，xi=1,2,10yi = 34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392（1）请用次数分别为3，4，5，6的多项式拟合并给出最好近似结果f(x)。（2）请用插值多项式给出最好近似结果下列数据为另外的对照记录，它们可以作为近似函数的评价参考数据。xi = Columns 1 through 7 1.5000 1.9000 2.3000 2.7000 3.1000 3.5000 3.9000 Columns 8 through 14 4.3000 4.7000 5.1000 5.5000 5.9000 6.3000 6.7000 Columns 15 through 17 7.1000 7.5000 7.9000yi = Columns 1 through 7 42.1498 41.4620 35.1182 24.3852 11.2732 -1.7813 -12.3006 Columns 8 through 14 -18.1566 -17.9069 -11.0226 2.0284 19.8549 40.3626 61.0840 Columns 15 through 17 79.5688 93.7700 102.36772.2 计算结果分析（1）根据实验数据，分别使用3,4,5,6次多项式拟合，可得拟合的多项式如下，并分别作图：n=3；p=-1.0326 19.334-94.479 131.79 n=4；p=-0.381857.368-42.143 73.5330.74498n=5；p=0.098075-3.0789 34.502-163.51 304.73-139.5n=6;p=0.019359 -0.540795.1137-16.897 -0.8669666.375 -18.699 上四幅图是多项式拟合曲线与原散点连线曲线的比较。明显可以看到，随着多项式拟合次数的增加，曲线的拟合度也越来越好，基本与连点曲线一致。最好近似结果是n=6时多项式拟合，表达式为：fx=0.019359*x6-0.54079*x5+5.1137*x4-16.897*x3-0.86696*x2+66.375*x1-18.699（2）根据原散点构建拉格朗日插值多项式，对给出的xi的数值进行插值得到yi，在下图可以看到利用拉格朗日插值多项式插值的图像与曲线拟合度很好。yi = Columns 1 through 13 42.3840 41.4947 35.0742 24.3601 11.2792 -1.7683 -12.2977 -18.1626 -17.9118 -11.0210 2.0333 19.8565 40.3584 Columns 14 through 17 61.0794 79.5709 93.7788 102.37133 第二题3.1 题目用雅格比法与高斯赛德尔迭代法解下列方程组Ax=b，研究其收敛性，上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。(1)A行分别为A1=6,2,-1,A2=1,4,-2,A3=-3,1,4； b1=-3,2,4T, b2=100,-200,345T,(2) A行分别为A1=1,0,8,0.8,A2=0.8,1,0.8,A3=0.8,0.8,1；b1=3,2,1 T, b2=5,0,-10T,(3)A行分别为A1=1,3,A2=-7,1；b=4,6T,3.2 结果分析在以下的程序中，默认初始值为0，最大迭代次数N为50，精度为0.00001（1）根据题意可得方程：对此方程分别采用两种方法，通过运行jacobi.m与gausssidel.m程序进行收敛性分析：Jacobi：满足误差限的迭代近似解为：-0.7273 0.8081 0.2525需要迭代18次；Gausssidel：满足误差限的迭代近似解为：-0.7273 0.8081 0.2525需要迭代11次。 Jacobi：满足误差限的迭代近似解为：36.3636 -2.0707 114.0404需要迭代24次；Gausssidel：满足误差限的迭代近似解为： 36.3636 -2.0707 114.0404需要迭代16次。 通过以上计算分析，此方程的解收敛。根据迭代次数的比较，可以得到Gausssidel迭代法要比Jacobi收敛性更快，更好。同时，上述两个方程系数阵都相同，但右端不同，用同样方法不同右端比较，得到结论，右端对迭代收敛有影响。（2）用jacobi到最大迭代次数50次时仍不能收敛。用Gausssidel满足误差限的迭代近似解为：5.7692 0.7692 -4.2308需要迭代16次。用jacobi到最大迭代次数50次时仍不能收敛。用Gausssidel满足误差限的迭代近似解为：32.6923 7.6923 -42.3077需要迭代43次。通过上述计算分析，用jacobi迭代时，两个方程都不能收敛，用Gaussssidel迭代时，可以得到收敛的近似解，既，Gausssidel方法具有良好的收敛性。（3）对此方程组，无论jacobi还是Gausssidel迭代法都不能得到收敛的近似解。可以推断是给出的行列式无法满足这两种方法的收敛条件。系数阵既不是行对角占优矩阵，也不是列占优阵，因此用这两种方法，得不到收敛的近似解。4 第三题4.1 题目用Runge-Kutta 4阶算法对初值问题y/=-20*y，y(0)=1按不同步长求解，用于观察稳定区间的作用，推荐两种步长h=0.1,0.2。注：此方程的精确解为：y=e-20x4.2 计算结果分析取不同步长得到求解的微分方程的解同精确解的误差分别如表1、表2、表3、表4、表5。从上表中可以得出，在h=0.0001,0.001,0.01,0.1 时，绝对误差在逐步变小，因此在这种情况下认为此方法是稳定的。当h=0.2 时，绝对误差逐渐扩大，因此在h=0.2 时，算法不稳定。在不同h下面的精度问题，取不同的步长求y(0.1)的值。从表6中可以看到步长越小，标准四阶Runge-Kutta 算法得到的值精度越高，计算所需的步骤也越长。附录一 多项式拟合程序代码%多项式拟合xx= 1 2 3 4 5 6 7 8 9 10;yy= 34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392;n=input(n=);%判断x和y向量大小是否一致if isequal(size(xx),size(yy) errorendx = xx(:);y = yy(:);%构建范德蒙德矩阵V(:,n+1) = ones(length(x),1,class(x);for j = n:-1:1 V(:,j) = x.*V(:,j+1)end%求解正规方程组的唯一解%所得的解称为最小二乘多项式%对矩阵V进行有选择的qr分解.%当矩阵V为mn并且mn,那么只会产生具有前n列的正交矩阵QQ,R = qr(V,0); p = R(Q*y); % p的保存 p = Vy;r = y - V*p;p = p.; %p的转置 一般采用行向量表示多项式系数%S结构包括cholesky分解的范德蒙德矩阵R，自由度df，标准残差normrS.R = R;S.df = length(y) - (n+1);S.normr = norm(r);%画图X=xx;Y=yy;Y1=polyval(p,xx);plot(xx,yy,:o,xx,Y1,-*);legend(yuanshi,nihe);附录二 多项式插值程序代码xx=1 2 3 4 5 6 7 8 9 10;yy=34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392;xi=1.5000 1.9000 2.3000 2.7000 3.1000 3.5000 3.9000 4.3000 4.7000 5.1000 5.5000 5.9000 6.3000 6.7000 7.1000 7.5000 7.9000;yi=1.5000 1.9000 2.3000 2.7000 3.1000 3.5000 3.9000 4.3000 4.7000 5.1000 5.5000 5.9000 6.3000 6.7000 7.1000 7.5000 7.9000%xi为标量或向量，被估计函数的自变量；%yi为xi处的函数估计值。x=xx(:);y=yy(:);xi=xi(:);n=length(x);m=length(y);%插值点与它的函数值应有相同个数if n=merror(The lengths of X and Y must be equal!);return;endyi=zeros(size(xi);for k=1:nw=ones (size(xi);for j=1:k-1 k+1:n;%输入的插值节点必须互异if abs(x(k)-x(j)eps;error(the data is error);return;endw=(xi-x(j)/(x(k)-x(j).*w; %计算Lagrange基函数endyi=yi+w*y(k); %计算Lagrange插值函数endyi%画图X=xi;Y1=load(f:matlab7shuzhifenxiY1.txt);Y2=yi(:);plot(xi,Y1,:o,xi,Y2,-*);legend(yuanshi,chazhi);附录三 jacobi迭代程序代码a=input(请输入系数矩阵a：);b=input(请输入矩阵b:);N=input(请输入最大迭代次数N：);esp=input(请输入近似解的误差限：);if any(diag(a)=0 error(系数矩阵错误，迭代终止！)endD=diag(diag(a);X0=zeros(size(b);x1=0;x2=0;x3=0;X1=x1;x2;x3;h=inv(D)*b;B=inv(D)*(D-a);B1=triu(B);B2=tril(B);k=1;fprintf(雅可比迭代法 );t=0;Y0=zeros(size(b);while t=N Y1=h+B*Y0; if norm(Y1-Y0,inf)esp fprintf(满足误差限 n) break ; end Y0=Y1; fprintf(第%2d次迭代得：,t) disp(Y1); t=t+1; endfprintf(满足误差限的雅可比迭代近似解为：)disp(Y1); 附录四 Gausssidel 迭代程序代码a=input(请输入系数矩阵a：);b=input(请输入矩阵b:);N=input(请输入最大迭代次数N：);esp=input(请输入近似解的误差限：);if any(diag(a)=0 error(系数矩阵错误，迭代终止！)endD=diag(diag(a);X0=zeros(size(b);x1=0;x2=0;x3=0;X1=x1;x2;x3;h=inv(D)*b;B=inv(D)*(D-a);B1=triu(B);B2=tril(B);k=1;fprintf(高斯-赛德尔迭代法 n);fprintf(第0次迭代得：)disp(X1);while k=N x1=h(1,1)+B1(1,:)*X0; X1=x1;x2;x3; x2=h(2,1)+B1(2,:)*X0+B2(2,:)*X1; X1=x1;x2;x3; x3=h(3,1)+B2(3,:)*X1; X1=x1;x2;x3; if norm(X1-X0,inf)esp fprintf(已满足误差限。 ) break ; end X0=X1; fprintf(第%2d次迭代得：,k) disp(X1); k=k+1; endfprintf(满足误差限的高斯-赛德尔迭代近似解为：)disp(X1);附录五 四阶Runge-Kutta算法程序代码h=0.0001,0.001,0.01,0.1k=1for xk=0.1:0.1:1 %取不同的x进行计算for i=1:1:4 %取不同的步长计算每个xhj=h(i)yk=1for xkp=hj:hj:xk %四阶标准Runge-Kutta迭代算法k1=hj*(-20*yk);format longk1;k2=hj*(-20*(yk+k1*0.5