2014数值分析MATLAB上机实验

1、数值分析实习报告姓名：gestepoA学号：201*班级：* 班序言随着计算机技术的迅速发展， 数值分析在工程技术领域中的应用越来越广泛， 并且成为数学与计算机之间的桥梁。 要解决工程问题， 往往需要处理很多数学模 型，不仅要研究各种数学问题的数值解法， 同时也要分析所用的数值解法在理论 上的合理性，如解法所产生的误差能否满足精度要求： 解法是否稳定、 是否收敛 及熟练的速度等。而且还能减少大量的人工计算。由于工程实际中所遇到的数学模型求解过程迭代次数很多，计算量很大，所 以需要借助如MATLAB C+, VB, JAVA的辅助软件来解决，得到一个满足误差限 的解。本文所计算题目，均采用 MA

2、TLABS行编程，MATLA被称为第四代计算机 语言，利用其丰富的函数资源，使编程人员从繁琐的程序代码中解放出来MATLAB最突出的特点就是简洁，它用更直观的、符合人们思维习惯的代码。 它具有以下优点：1友好的工作平台和编程环境。MATLAB面精致，人机交互性强，操作简单。2简单易用的程序语言。MATLAB!个高级的矩阵/阵列语言，包含控制语 言、函数、数据结构，具有输入、输出和面向对象编程特点。用户可以在命令窗 口中将输入语句与执行命令同步，也可以先编好一个较大的复杂的应用程序（ M 文件）后再一起运行。3 强大的科学计算机数据处理能力。包含大量计算算法的集合，拥有 600 多 个工程中要用

3、到的数学运算函数。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,拟分别用插值多项式和多项式拟合给出

4、其对应规律的 近似多项式，已知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试验原理和方法:实验原理：拉格朗日插值法中先构造插值基础函数:?= n?=o?初?加（?

5、= 0,12? , ?，然后构 j *k?造出拉格朗日多项式：?（?=空=。（n?=0? ?。j r ? ?最佳平方逼近中，设逼近函数？?？ ?=? + ?+ ? + ? + ?,逼近函数和真实函? 1?数之差 r = ? - ? 1 ? = ?蔦-?,即:?= ? ?根据最小二?1?乘准则令En?=o? = ?,?可以得到？?=仔的?-1 ?实验方法：逼近法采用最佳平方逼近，依据最小二乘原则：马?=0? = ?,?油已知条件采用离散型。插值法采用拉格朗日插值法。在逼近法中，由于是离散型的，所以法方程系数阵设计成求和。分别求出3、4、5、6次的多项式，逼近结果和真实值有一定差距， 最小二乘正是

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.

7、033*xA3 + 19.33*xA2 - 94.48*x + 131.8拟合结果:12345678910x12345678910y55.61711.896-5.561-2.95213.52537.67263.29184.18494.15387.0000000000000四次多项式：-0.3818*xA4 + 7.368*xA3 - 42.14*乂人2 + 73.53*x+ 0.745 拟合结果：12345678910x12345678910y39.12132.08010.085-5.563-2.73021.56061.117100.588115.45772.045222802222017 I

8、10D83B-l102?五次多项式：0.09807*L5 - 3.079P4 + 34.5*tA3 -163.5*2 + 304.7*t -139.5 拟合结果：12345678910x12345678910y33.21945.7749.032-16.500-8.90626.90870.983100.07399.41674.5001203335840六次多项式：0.01936*tA6 - 0.5408*tA5 + 5.114十4 - 16.93 - 0.8672 + 66.38*t - 18.7 拟合结果：12345678910x12345678910y34.50541.14913.270-1

9、3.94-12.2223.71173.950105.1696.34578.4006408650409460对比可知，六次多项式拟合结果最好。1.3.2拉格朗日插值法插值多项式 5.353*10A( -5)*xA9 - 0.003088*xA8 + 0.07229*xA7 - 0.8792*xA6 + 5.932*xA5-22.41*xA4 + 50.1计3 - 86.47*乂人2 + 113.5*x - 25.2注：此多项式为拉格朗日多项式的近似式，当x=10的时候偏差可以达到 23以上。对比数据：123456789x1.50001.90002.30002.70003.10003.50003.

10、90004.30004.7000y42.149841.46235.11824.38511.273-1.7813-12.300-18.1566-17.9060222691011121314151617x5.10005.50005.90006.30006.70007.10007.50007.9000y-11.0222.028419.85440.36261.08479.56893.7700102.367696087插值结果:123456789x1.50001.90002.30002.70003.10003.50003.90004.30004.7000y42.384041.49435.07424.36

11、011.279-1.7683-12.297-18.1626-17.9117212781011121314151617x5.10005.50005.90006.30006.70007.10007.50007.9000y-11.0212.033319.85640.35861.07979.57093.7788102.371054493其中红点表示参考点。1.3.3比较选取xi=a+ih=1+0.2*i(i=0,1, ? ,45)，分别绘制六次多项式拟合和拉格朗日插值结果图:120其中绿线表示拉格朗日插值多项式图像，蓝线表示六次多项式拟合图像。两者效果近似但后者比前者低三次。2第二题用雅格比法与高斯-

12、赛德尔迭代法解下列方程组 Ax=b1或Ax=b2，研究其收敛性。 上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无 影响。(1) A行分别为 A仁6,2, - 1, A2=1,4, - 2, A3=- 3,1,4 ;b1=- 3,2,4 T;b2=100,-200,345 T。(2) A行分别为 A1=1,0,8,0.8,A2=0.8,1,0.8,A3=0.8,0.8,1; b1=3,2,1T;b2=5,0, - 10T。 A行分别为 A仁1,3, A2=-7,1 ; b1=4,6 T。2.1试验目的学习jacobi迭代法和GuassSeidel迭代法的原理和编程方法，研

13、究方程组系数阵和右边 项对方程的解及其收敛性的影响，判断迭代法的收敛条件。2.2 实验原理和方法实验原理：将方程组系数阵 A分解为A= D+ L+ U,其中D为对角阵，L为减去D的下三角阵，U 为减去 D 的上三角阵。Jacobi 迭代法中构造如下迭代公式：?( ?+1) = -? -1 (?+ ?) ?( ?) + ?-1 ? 而 Gauss-Seidel 迭代法的迭代公式为：?(?+1) = - (?+ ?)-1 ?(?) + (?+ ?)-1 ?/ f初始值直接选取为0。在判断其收敛性时，分别求解其迭代矩阵的谱半径p?, P? =/?为迭代矩阵的特征值。实验方法：分别编写 jacobi

14、迭代及其收敛判别函数和 Seidel 迭代及其收敛判别函数。 如果在初试迭 代步数之内还未收敛就进行收敛判别，收敛判别的依据是迭代矩阵的谱半径是否小于1。比较同一方程组的 jacobi 迭代法和 Seidel 迭代法的结果是否相同， 在达到精度要求后比 较两种方法的迭代次数，比较哪一个的效率更高。比较方程组系数阵和等号右边的变化会对方程的解和收敛速度造成什么影响。 如果迭代不收敛，那么考虑为什么不收敛，如果把方程组系数阵进行强对角占优处理， 是否会收敛。2.3 实验结果规定误差界： 1e-42.3.1 第一问62-1-3 ?= 14-2 ,?= 2 .-3144-0.7273由 jacobi

15、迭代法求得?= 0.8081，设定迭代 20 次，实际迭代 16 次，精度为 9.4022e-005。0.2525-0.7272由 seidel 迭代法求得 ?= 0.8081 ，设定迭代 20 次，实际迭代 10 次，精度为：0.25259.0769e-00562-1100 ?= 14-2 ，?= -200 -31434536.3636由 jacobi 迭代法求得 ?= -2.0707 ，设定迭代 40 次，实际迭代 23 次，精度为114.04046.9948e-005。36.3637由 Seidel 迭代法求得 ?= -2.0707 ，设定迭代 20 次，实际迭代 15 次，精度为114

16、.04048.6384e-005通过对比可知： 1、 Seidel 迭代的收敛速度明显高于 jacobi 迭代。 2、 b 矩阵对收敛速 度和误差精度有影响， b 中元素较大时会放慢收敛速度并加大误差。2.3.2 第二问1 ?= 0.80.80.8 0.831 0.8，?= 20.811由 jacobi 迭代法求解， 100 次迭代尚且不能达到精度。此时调用 jacobi 迭代法的收敛判别函数，求得特征值为：1 = -1.6 ,2、3 = 0.8 , P (?= 1.6 1，迭代不收敛。5.7691由于 Seidel 迭代法求解 ?= 0.7693 ，迭代次数 31，精度为 8.7826e-0

17、05-4.230710.80.85 ?= 0.810.8，?= 0 0.80.81-10Jacobi 迭代不收敛。32.6922Seldel 迭代法求得 ?= 7.6922 ，迭代次数 38，精度为 8.4552e-005。-42.3076比较得知 A矩阵元素如果相差很小，迭代次数会大幅增加，综合比较可知b矩阵元素如果相差很大会增加迭代次数。2.3.3 第三问试验结果和讨论?= 1-731，?= 4?= 6此时 p （第=4.5826 1 ,p(?3= 21 1 , jacobi迭代法和Seidel迭代法都不收敛。如果交换 A 中行的顺序， 得到 -71；，用jacobi迭代计算，迭代8次，解

18、得x = I；：64 。用Seidel迭代法计算，只需迭代5次，得到x = -0.6364 ，精度为2.6326e-005。此时p (?= 1.54550.2182 p(?) = 0.0476 ，从此可以看出收敛速度的快慢。3第三题1给定函数f(x)r，5x5,及节点Xi 5 i, i 0,1,L 10,求其三次样1 x2条插值多项式(可取I型或II型边界条件)，并画图及与f(x)的图形进行比较分析。注：涉及到线性方程组求解问题需采用适当的求解算法。3.1实验目的学习三弯矩法的原理和编程方法，对比原函数和三次样条插值的结果。3.2实验原理和方法实验原理:记 h?= ? ?_1， ?尸?+1CC

19、6?+1? ?_1，人丄十 /+?= 1 - ? ?=厂（一 -?），给出插值?+? ?+1?+1re? (?) = ?) ?(? = ?(? = ?。组成递推式：?_1 + 2?+ ?+1= ?由于系数阵按行对角占优，方程组存在唯一确定解，可以使用高斯列主元消去法来解 方程。最后将各个参数带入样条函数(?_ ?3(? ?_1)3?_1 2?= ?_1 二 + ? + (?菊1_厂？?)6? ? 6? ? 6即可求得样条函数。两端处的二阶导数?(?)=? ? 2-7+ （?毎-?）?6?- ?_1?实验方法由于在本题中x(i+1)_x(i)=1,所以h(i)=1。在编程中直接将 h设置成常数，

20、简化了运算。首先求解、g，然后列出方阵求解 M(i)，在求解方程组的过程中采用列主元素高斯消去法，分为消元和回代两个过程，编写这两个函数，解出除了两端的M，而两端点的 M值等于两端点的函数二阶导数值。编写函数求出样条函数的系数，然后求出方程，对于三弯矩法三次样条函数，如果有n个点，则有n_1个样条函数，除了两端需要求解n_2个M值，即解n_2阶方程组。在表达样条函数的时候采用if语句，对不同的区间进行划分， 然后细分_5,5这个区间，间隔0.1将其分为100份，这样可以体现出连续性，此时绘图对比三次样条函数和原函数。本题中原函数的二阶导不是计算机解出的没有编写相关程序。本题中求解的样条函数，M

21、ATLAB系统自动的将公因式提取，并且合并同类项，所以表达出的函数并不整齐规律，为了更好地体现三次样条函数的结构和性质，我专门手写了规整的样条函数。3.3实验结果三弯矩法求解代入 x=-5-4-3-2-10123451 y = f(?=21 + ?0.0588? ( - 5 ) = ? (5)0.00842y=0.03850.05880.10000.20000.50001.00000.50000.20000.10000.0385S=M=0.01410.05990.09930.7431-1.87150.74310.09930.05990.01410.00140.00240.03710.05650

22、.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.09000.05650.00240.00140.05650.0371求得样条函数系数阵为：样条函数：(97*X)/5000 - (7*(X + 4)A3)/5000 + (3*(X + 5)A3)/1250 + 1341/10000(

23、67*X)/2000 - (3*(X + 3)A3)/1250 + (X + 4)人3/100 + 381/2000 (187*X)/2000 - (X + 2)A3/100 + (33*(X + 3)人3)/2000 + 741/2000 (1927*X)/10000 - (33*(X + 1)人3)/2000 + (619*(X + 2)人3)/5000 + 5689/10000 (9357*X)/10000 - (3119*(X + 1)人3)/10000 - (619*乂人3)/5000 + 13119/10000 (3199*(X - 1)A3)/10000 - (9357*X)/1

24、0000 + (619*乂人3)/5000 + 13119/10000 (33*(X - 1)A3)/2000 - (1927*X)/10000 - (619*(X - 2)人3)/5000 + 5689/10000 (X - 2)A3/100 - (187*X)/2000 - (33*(X - 3)人3)/2000 + 741/2000 (3*(X - 3)A3)/1250 - (67*X)/2000 - (X - 4)人3/100 + 381/2000 (7*(X - 4)A3)/5000 - (97*X)/5000 - (3*(X - 5)人3)/1250 + 1341/10000 写出

25、标准形式的样条函数：-4-x)+0.0565*(x+5)-3-x)+0.0900*(x+4)-2-x)+0.1835*(x+3)-1 -x)+0.3762*(x+2)x-5,-4x-4,-3x-3,-2x-2,-1 x-1,0 x0,1S1(x)= 0.0014*( -4-x)A3+0.0024*(x+5)A3+0.0371*(S2(x)= 0.0024*( -3-x)A3+0.0100*(x+4)A3+0.0565*(S3(x)= 0.0100*( -2-x)A3+0.0165*(x+3)A3+0.0900*(S4(x)= 0.0165*( -1 -x)A3+0.1238*(x+2)A3+0

26、.1835*(S5(x)= 0.1238*( -x)A3-0.3119*(x+1)A3+0.3762*( -x)+1.3119(x+1)S6(x)= -0.3119*(1 -x)A3+0.1238*xA3+1.3119*(1 -x)+0.3762*xS7(x)= 0.1238*(2-x)A3+0.0165*(x-1)人3+0.3762*(2-x)+0.1835*(x-1)x1,2S8(x)= 0.0165*(3-x)A3+0.0100*(x-2)人3+0.1835*(3-x)+0.0900*(x-2)x2,3S9(x)= 0.0100*(4-x)A3+0.0024*(x-3)人3+0.0900

27、*(4-x)+0.0565*(x-3)x3,4S10(x)= 0.0024*(5 -x)A3+0.0014*(x -4)人3+0.0565*(5 -x)+0.0371*(x -4)x4,5输入：a=-5:0.1:5;for i=1:le ngth(a) b(i)=f(a(i);end求得原函数和样条函数的对比图:4 MATLAB 程序第一题：离散型最佳平方逼近函数fun cti onf=squar_approx_ls(x, y,n)syms t ;N=le ngth(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(k)AiendQ(i+1,1)=b;b=0; end s=inv(P)*Q; f=0; for i=1:n+1f=f+s(i)*tq-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: