数值分析MATLAB上机实验

数值分析实习报告姓名： gestepoA学号：201*班级：序言随着计算机技术的迅速发展， 数值分析在工程技术领域中的应用越来越广泛，并且成为数学与计算机之间的桥梁。 要解决工程问题， 往往需要处理很多数学模型， 不仅要研究各种数学问题的数值解法， 同时也要分析所用的数值解法在理论上的合理性， 如解法所产生的误差能否满足精度要求： 解法是否稳定、 是否收敛及熟练的速度等。而且还能减少大量的人工计算。由于工程实际中所遇到的数学模型求解过程迭代次数很多，计算量很大，所以需要借助如 MATLAB ， C+ ， VB ， JAVA 的辅助软件来解决，得到一个满足误差限的解。本文所计算题目，均采用 MATLAB 进行编程， MATLAB 被称为第四代计算机语言， 利用其丰富的函数资源， 使编程人员从繁琐的程序代码中解放出来MATLAB 最突出的特点就是简洁， 它用更直观的、 符合人们思维习惯的代码。它具有以下优点：1 友好的工作平台和编程环境。 MATLAB 界面精致， 人机交互性强， 操作简单。2 简单易用的程序语言。 MATLAB 是一个高级的矩阵/阵列语言，包含控制语言、函数、数据结构，具有输入、输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编好一个较大的复杂的应用程序（ M 文件）后再一起运行。3 强大的科学计算机数据处理能力。 包含大量计算算法的集合， 拥有 600 多个工程中要用到的数学运算函数。4 出色的图像处理功能， 可以方便地输出二维图像， 便于我们绘制函数图像。目录1 第一题 41.1 实验目的 41.2 实验原理和方法 41.3 实验结果 51.3.1 最佳平方逼近法51.3.2 拉格朗日插值法 71.3.3 对比 82 第二题 92.1 实验目的 92.2 实验原理和方法102.3 实验结果 102.3.1 第一问 102.3.2 第二问 112.3.3 第三问 113 第三题 123.1 实验目的 123.2 实验原理和方法123.3 实验结果 124 MATLAB 程序 141 第一题某过程涉及两变量x 和 y， 拟分别用插值多项式和多项式拟合给出其对应规律的近似多项式，已知xi 与 yi 之间的对应数据如下：12345678910x12345678910y34.640.314.6-14.2-13.324.875.2103.597.478.2588719448721570234795743847392请用次数分别为3, 4, 5, 6的多项式拟合并给出最好近似结果f(x)请用插值多项式给出最好近似结果。1.1 实验目的：学习逼近和插值的原理和编程方法，由给出的已知点构造多项式， 在某个范围内近似代替已知点所代表的函数，以便于简化对未知函数的各种计算。1.2 试验原理和方法：实验原理：拉格朗日插值法中先构造插值基础函数：-，然后构造出拉格朗日多项式：。最佳平方逼近中，设逼近函数，逼近函数和真实函数之差，即：，根据最小二乘准则令，可以得到。实验方法：逼近法采用最佳平方逼近，依据最小二乘原则：，由已知条件采用离散型。插值法采用拉格朗日插值法。在逼近法中，由于是离散型的，所以法方程系数阵设计成求和。分别求出3、4、5、6次的多项式，逼近结果和真实值有一定差距,理论上多最小二乘正是让这些差距达到最小，项式次数越高结果和真实值差距越小。拉格朗日插值法中 la=la*(p-x(j)/(x(k)-x(j) ”语句实现的是我们通常书写的连乘形式拉格朗日插值多项式，但是表示不方便，而如果用“ s=collect(s) ”函数将其展开成降哥排列多项式以后，由于余项问题结果会和原本的多项式有偏差，这种偏差随着x的增大而增大。求出多项式后和题目中给出的参考点进行比较。最后，选择六次最佳平方逼近多项式和拉格朗日插值多项式(九次)进行比较，选取xi=a+ih=1+0.2*i(i=0,1,45),分别绘制两者的图像进行比较。1.3 试验结果1.3.1 最佳平方逼近法三次多项式：-1.033*人3 + 19.33*xA2 - 94.48*x + 131.8拟合结果：12345678910x12345678910y55.6111.89-5.56-2.9513.5237.6763.2984.1894.1587.0070601020502010403000120100* *2011111111112E4E57B910四次多项式：-0.3818*xA4 + 7.368*xA3 - 42.14*xA2 + 73.53*x+ 0.745拟合结果:12345678910x12345678910y39.1232.0810.08-5.56-2.7321.5661.11100.58115.4572.041202523800027282725017 I100SO6 口402。0五次多项式：0.09807*tA5 - 3.079*tM + 34.5*tA3 - 163.5*tA2 + 304.7*t - 139.5拟合结果：12345678910x12345678910y33.2145.779.03-16.50-8.9026.9070.98100.0799.4174.5091422003638335386400六次多项式：0.01936*tA6 - 0.5408*5 + 5.114*tA4 - 16.9*tA3 - 0.867*tA2 + 66.38*t -18.7拟合结果:12345678910x12345678910y34.5041.1413.27-13.94-12.2223.7173.95105.1696.3478.4056940086501400945600对比可知，六次多项式拟合结果最好。1.3.2 拉格朗日插值法插值多项式 5.353*10人(-5)*人9- 0.003088*xA8+ 0.07229*人7- 0.8792*人6+5.932*xA5 - 22.41*xA4 + 50.133 - 86.47*人2 + 113.5*x - 25.2注：此多项式为拉格朗日多项式的近似式，当x=10的时候偏差可以达到 23以上。对比数据:123456789x1.50001.90002.30002.70003.10003.50003.90004.30004.7000y42.14941.46235.11824.38511.273-1.781-12.300-18.156-17.9068022236691011121314151617x5.10005.50005.90006.30006.70007.10007.50007.9000y-11.0222.028419.85440.36261.08479.56893.770102.3676960807插值结果:123456789x1.50001.90002.30002.70003.10003.50003.90004.30004.7000y42.38441.49435.07424.36011.279-1.768-12.297-18.162-17.9110721237681011121314151617x5.10005.50005.90006.30006.70007.10007.50007.9000y-11.0212.033319.85640.35861.07979.57093.778102.371o544983其中红点表示参考点。1.3.3比较选取 xi=a+ih=1+0.2*i(i=0,1,45),分别绘制六次多项式拟合和拉格朗日插值结果图:120其中绿线表示拉格朗日插值多项式图像，蓝线表示六次多项式拟合图像。两者效果近似但后者比前者低三次。2第二题用雅格比法与高斯赛德尔迭代法解下列方程组 Ax=b1或Ax=b2,研究其收敛 性。上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛 有无影响。A 行分别为 A1=6,2,- 1,A2=1,4,- 2,A3=- 3,1,4; b1=-3,2,4为b2=100,- 200,345 To(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,- 10TOA 行分别为 A1=1,3,A2=-7,1 ; b1=4,6T2.1 试验目的学习 jacobi 迭代法和 GuassSeidel 迭代法的原理和编程方法，研究方程组系数阵和右边项对方程的解及其收敛性的影响，判断迭代法的收敛条件。2.2 实验原理和方法实验原理：将方程组系数阵A 分解为 ,其中 D 为对角阵， L 为减去 D 的下三角阵， U为减去 D 的上三角阵。Jacobi 迭代法中构造如下迭代公式：而 Gauss-Seidel 迭代法的迭代公式为：初始值直接选取为 0 。在判断其收敛性时，分别求解其迭代矩阵的谱半径， 为迭代矩阵的特征值。实验方法：分别编写 jacobi 迭代及其收敛判别函数和Seidel 迭代及其收敛判别函数。如果在初试迭代步数之内还未收敛就进行收敛判别，收敛判别的依据是迭代矩阵的谱半径是否小于1 。比较同一方程组的 jacobi 迭代法和 Seidel 迭代法的结果是否相同，在达到精度要求后比较两种方法的迭代次数，比较哪一个的效率更高。比较方程组系数阵和等号右边的变化会对方程的解和收敛速度造成什么影响。如果迭代不收敛，那么考虑为什么不收敛，如果把方程组系数阵进行强对角占优处理，是否会收敛。2.3 实验结果规定误差界： 1e-42.3.1 第一问,由 jacobi 迭代法求得9.4022e-005 。由 seidel 迭代法求得9.0769e-005，由 jacobi 迭代法求得6.9948e-005 。由 Seidel 迭代法求得20 次，实际迭代 16 次，精度为20 次，实际迭代10 次，精度为：40 次，实际迭代 23 次，精度为20 次，实际迭代15 次，精度为8.6384e-005通过对比可知： 1 、 Seidel 迭代的收敛速度明显高于jacobi 迭代。 2 、 b 矩阵对收敛速度和误差精度有影响， b 中元素较大时会放慢收敛速度并加大误差。2.3.2 第二问 ，由jacobi迭代法求解，100次迭代尚且不能达到精度。此时调用 jacobi迭代法的收敛 判别函数，求得特征值为：，迭代不收敛。、由于Seidel迭代法求解，迭代次数31 ,精度为8.7826e-005 ，Jacobi迭代不收敛。Seldel迭代法求得，迭代次数38,精度为8.4552e-005 。比较得知 A矩阵元素如果相差很小，迭代次数会大幅增加，综合比较可知b矩阵元素如果相差很大会增加迭代次数。2.3.3 第三问试验结果和讨论此时，jacobi迭代法和Seidel迭代法都不收敛。如果交换A中行的顺序，得到，用jacobi迭代计算，迭代8次，解得用Seidel迭代法计算，只需迭代 5次，得到，精度为2.6326e-005 。此时,从此可以看出收敛速度的快慢。3第三题1给止函数f(x)=2， -5 MxM5 ,及下点为=-5+i, i = 0,1川10 ,求其二次样1 x条插值多项式(可取I型或II型边界条件)，并画图及与f(x)的图形进行比较分 析。注：涉及到线性方程组求解问题需采用适当的求解算法。3.1 实验目的学习三弯矩法的原理和编程方法，对比原函数和三次样条插值的结果。3.2 实验原理和方法实验原理：记，给出插值两端处的二阶导数。组成递推式：由于系数阵按行对角占优，方程组存在唯一确定解，可以使用高斯列主元消去法来解方程。最后将各个参数带入样条函数即可求得样条函数。实验方法由于在本题中x(i+1)-x(i)=1,所以h(i)=1。在编程中直接将 h设置成常数，简化了运算。首先求解法、g,然后列出方阵求解 M(i),在求解方程组的过程中采用列主元素高斯消去法，分为消元和回代两个过程，编写这两个函数，解出除了两端的M,而两端点的 M值等于两端点的函数二阶导数值。编写函数求出样条函数的系数，然后求出方程，对于三弯矩法三次样条函数，如果有n个点，则有n-1个样条函数，除了两端需要求解n-2个M值，即解n-2阶方程组。在表达样条函数的时候采用if语句，对不同的区间进行划分，然后细分-5,5这个区间，间隔0.1将其分为100份，这样可以体现出连续性， 此时绘图对比三次样条函数和原函数。本题中原函数的二阶导不是计算机解出的没有编写相关程序。本题中求解的样条函数，MATLAB系统自动的将公因式提取，并且合并同类项，所以表达出的函数并不整齐规律，为了更好地体现三次样条函数的结构和性质，我专门手写了规整的样条函数。3.3 实验结果三弯矩法求解代入 x=-5-4 -3 -2 -10 1 2 3 4 5y=0.03850.05880.10000.20000.50001.00000.50000.20000.10000.0588 0.0385M=0.01410.05990.09930.7431-1.87150.74310.09930.05990.0141求得样条函数系数阵为：S=0.00140.00240.03710.05650.00240.01000.05650.09000.01000.01650.09000.18350.01650.12380.18350.37620.1238-0.31190.37621.3119-0.31190.12381.31190.37620.12380.01650.37620.18350.01650.01000.18350.09000.01000.00240.00240.00140.09000.05650.05650.0371样条函数：(97*X)/5000 - (7*(X + 4)人3)/5000 + (3*(X + 543)/1250 + 1341/10000 (67*X)/2000 - (3*(X + 3)A3)/1250 + (X + 4)A3/100 + 381/2000(187*X)/2000 - (X + 2)A3/100 + (33*(X + 3)A3)/2000 + 741/2000(1927*X)/10000 - (33*(X + 1)A3)/2000 + (619*(X + 2)A3)/5000 + 5689/10000 (9357*X)/10000 - (3119*(X + 1)A3)/10000 - (619*XA3)/5000 + 13119/10000(3199*(X - 1)A3)/10000 - (9357*X)/10000 + (619*XA3)/5000 + 13119/10000 (33*(X - 1)A3)/2000 - (1927*X)/10000 - (619*(X - 2)A3)/5000 + 5689/10000 (X - 2)A3/100 - (187*X)/2000 - (33*(X - 3)A3)/2000 + 741/2000 (3*(X - 3)A3)/1250 - (67*X)/2000 - (X - 4)A3/100 + 381/2000(7*(X - 4)A3)/5000 - (97*X)/5000 - (3*(X - 5)A3)/1250 + 1341/10000写出标准形式的样条函数：S1(x)= 0.0014*(-4-x)A3+0.0024*(x+5)A3+0.0371*(-4-x)+0.0565*(x+5)x-5,-4S2(x)= 0.0024*(-3-x)A3+0.0100*(x+4)A3+0.0565*(-3-x)+0.0900*(x+4)x-4,-3S3(x)= 0.0100*(-2-x)A3+0.0165*(x+3)A3+0.0900*(-2-x)+0.1835*(x+3)x-3,-2S4(x)= 0.0165*(-1-x)A3+0.1238*(x+2)A3+0.1835*(-1-x)+0.3762*(x+2)x-2,-1S5(x)= 0.1238*(-x)A3-0.3119*(x+1)A3+0.3762*(-x)+1.3119(x+1)x-1,0S6(x)= -0.3119*(1-x)A3+0.1238*xA3+1.3119*(1-x)+0.3762*xx0,1S7(x)= 0.1238*(2-x)A3+0.0165*(x-1)A3+0.3762*(2-x)+0.1835*(x-1)S8(x)= 0.0165*(3-x)A3+0.0100*(x-2)A3+0.1835*(3-x)+0.0900*(x-2)S9(x)= 0.0100*(4-x)A3+0.0024*(x-3)A3+0.0900*(4-x)+0.0565*(x-3)S10(x)= 0.0024*(5-x)A3+0.0014*(x-4)A3+0.0565*(5-x)+0.0371*(x-4)输入：x1,2x2,3x3,4x4,5a=-5:0.1:5;for i=1:length(a)b(i)=f(a(i)；end4 MATLAB 程序第一题：离散型最佳平方逼近函数function f=squar_approx_ls(x,y,n) syms t;N=length(x);P=zeros(n+1);Q=zeros(n+1,1);a=0;for i=0:nfor j=0:nfor k=1:Na=a+x(kF(i+j);endP(i+1,j+1)=a;a=0;endendb=0;for i=0:nfor k=1:Nb=b+y(k)*x(kiendQ(i+1,1)=b;b=0;end s=inv(P)*Q;f=0;for i=1:n+1f=f+s(i)*tA(i-1);simplify(f);end f=collect(f);f=vpa(f,4);拉格朗日插值函数function s=Lagrange(x,y,x0)syms p;n=length(x);s=0;for(k=1:n)la=y(k);for(j=1:k-1)la=la*(p-x(j)/(x(k)-x(j);end;for(j=k+1:n)la=la*(p-x(j)/(x(k)-x(j);end;s=s+la;simplify(s);endif(nargin=2)s=collect(s);s=vpa(s,4);elsem=length(x0);for i=1:mtemp(i)=subs(s,p,x0(i);ends=temp;end第二题Jacobi 迭代法函数function x=jacobi(A,b,P,delta,n)N=length(b);for k=1:nfor j=1:Nx(j)=(b(j)-A(j,1:j-1,j+1:N)*P(1:j-1,j+1:N)/A(j,j);enderr=abs(norm(x-P);P=x;if(errdelta)break;endPendx=x;k,errjacobi 迭代收敛判别函数：function s=liansanxing(A)N,N=size(A);D=zeros(N);LU=zeros(N);for i=1:ND(i,i)=A(i,i);endLU=A-D;p=-inv(D)*LU;V,s=eig(p);Seidel 迭代法函数function x=Seidel(A,b,P,delta,n)N=length(b);for k=1:nfor j=1:Nif j=1x(1)=(b(1)-A(1,2:N)*P(2:N)/A(1,1);elseif j=Nx(N)=(b(N)-A(N,1:N-1)*(x(1:N-1)/A(N,N);elsex(j)=(b(j)-A(j,1:j-1)*x(1:j-1)-A(j,j+1:N)*P(j+1:N)/A(j,j);endenderr=abs(norm(x-P);P=x;if(err=-5&x=-4&x=-3&x=-2&x=-1&x=0&x=1&x=2&x=3&x=4&x=5y=0.0024*(5-x)A3+0.0014*(x-4)A3+0.0565*(5-x)+0.0371*(x-4);endWelcome ToDownload !欢迎您的下载，资料仅供参考！