《数值分析大作业(牛顿下山法-拉格朗日法-切比雪夫法)及Matlab程序.doc》由会员分享，可在线阅读，更多相关《数值分析大作业(牛顿下山法-拉格朗日法-切比雪夫法)及Matlab程序.doc（73页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、Four short words sum up what has lifted most successful individuals above the crowd: a little bit more.-author-date数值分析大作业(牛顿下山法-拉格朗日法-切比雪夫法)及Matlab程序数值分析大作业(牛顿下山法-拉格朗日法-切比雪夫法)及Matlab程序课 程 设 计课程名称： 数值分析 设计题目： 学 号： 姓 名： 完成时间： 2014.11.18 题目一： 解线性方程组的直接法设方程组，其中，矩阵中，由相应的矩阵元素计算，使解向量。(1) 不变，对的元素加一个扰动，求解方程
2、组；(2) 不变，对的元素和分别加一个扰动，求解方程组；(3) 对上述两种扰动方程组的解做误差分析。一.数学原理： 本计算采用直接法中的列主元高斯消元法，高斯列主元消元法原理如下：1、设有n元线性方程组如下：2、第一步：如果a11!=0, 令li1= ai1/a11, I= 2,3,n用（-li1）乘第一个方程加到第i个方程上,得同解方程组：a(1)11 a(1)12 . . . a(1)1n x1 b(1)1a(1)21 a(1)22 . . . a(1)2n x2 b(1)2 . . . . . . . = .a(1)n-11 a(1)n-12 . . a(1)n-1n xn-1 b(1)
3、n-1a(1)n1 a(1)n2 . . . a(1)nn xn b(1)n简记为： A(2) x = b(2)其中 a(2)ij = a(1)ij li1 * a(1)1j , I ,j = 2,3,.,n b(2)I = b(1)I li1 * b(1)1 , I = 2,3,.,n第二步：如果a(2)22 != 0,令li2= a（2）i2/a（2）22, I= 3,n依据同样的原理，对矩阵进行化间（省略），依次下去，直到完成！最后，得到上三角方程组： a(1)11 a(1)12 . . . a(1)1n x1 b(1)10 a(1)22 . . . a(1)2n x2 b(1)2 .
4、. . . . . . = .0 0 . . a(n-1)n-1n xn-1 b(n-1)n-10 0 . . . a(n)nn xn b(n)n简记为： A(n) x = b(n)最后从方程组的最后一个方程进行回代求解为： Xn = b(n) / a(n)nn Xi = ( b(k)k - a(k)kjxj ) / a(k)kk二解题过程：1.由题中所给条件可求出b。B = 6.0000 7.7156 9.9299 12.7560 16.3238 20.7813（1） 不变，对的元素加一个扰动，求解方程组。B =6.0000 7.7156 9.9299 12.7560 16.3238 20.
5、7813+0.0001解得x =0.5997 2.6920 -1.8500 3.3917 0.0000 1.1667（2）不变，对的元素和分别加一个的扰动，求解方程组。 A = 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.100001000000000 1.210000000000000 1.331000000000000 1.464100000000001 1.61051000000
6、0001 1.000000000000000 1.200000000000000 1.440000000000000 1.728000000000000 2.073600000000000 2.488319999999999 1.000000000000000 1.300000000000000 1.690000000000000 2.197000000000000 2.856100000000001 3.712930000000001 1.000000000000000 1.400000000000000 1.960000000000000 2.743999999999999 3.84159
7、9999999999 5.378239999999998 1.000000000000000 1.500000000000000 2.250000000000000 3.375000000000000 5.062500000000000 7.593751000000000x =0.825832305593523 1.742109746543727 -0.259537880024995 2.064585114452380 0.551832152872303 1.075178560563062三、误差分析：从上面计算结果可以看出，当系数矩阵或右端向量发生极小的扰动，方程组的解也会产生很大的误差，产
8、生的原因是范德蒙阵为变态阵。由数值计算知识可知 其中 为条件数，从上式看到，当A的条件数很大时，解的相对误差也很大，此时的对应的线性方程为病态线性方程组。计算条件数时，取矩阵的无穷范数，经计算得矩阵、受扰动后的矩阵和的条件数为；可以看到三个矩阵的条件数非常大，即使系数矩阵或右端向量发生很小的变化，也会导致解产生很大的误差。四、收获与体会：运用matlab编程解决数学问题很方便，病态阵的条件数非常大，给系数矩阵或者右端向量一个很微小的扰动，方程组的解也会产生很大的变动。通过做这个题目，我对让课上抽象病态现象有了直观的认识。题目二：多项式插值在区间上对龙格函数做插值，分别用等距节点(节点步长)和切
9、比雪夫多项式的零点做插值节点，画出原函数和两个插值函数的图像进行比较，并利用对两种插值方法做误差分析。一.数学原理：1.拉格朗日等距节点插值拉格朗日插值多项式为 其中 余项为 与有关。2.拉格朗日以切比雪夫零点为节点插值切比雪夫多项式 在区间 上有 个零点，为 由上式得到切比雪夫多项式的个零点，并作为插值节点做拉格朗日插值，得到多项式。一.计算过程：原函数为y=1./(1+x.2)，通过matlab绘图得到原函数图像为：1. 利用拉格朗日差分(等距节点)：x =(-5 -4 -3 -2 -1 0 1 2 3 4 5)y=1./(1+x.2)= (0.0385 0.0588 0.1000 0.2
10、000 0.5000 1.0000 0.5000 0.2000 0.1000 0.0588 0.0385)拉格朗日插值多项式为：=(t*(t/8 + 5/8)*(t - 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/4320 - (t*(t/15 + 1/3)*(t - 1)*(t + 1)*(t - 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/10080 - (t/5 + 1)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t
11、- 4)*(t + 4)*(t - 5)/2880 + (t*(t/20 + 1/4)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t - 4)*(t + 4)*(t - 5)/40320 + (t*(t/12 + 5/12)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/2880 - (t*(t/17 + 5/17)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t - 5)/362880 + (t*(t/
12、26 + 2/13)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t - 5)/3628800 - (t*(t/35 + 1/7)*(t - 1)*(t + 1)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/4320 + (t*(t/80 + 1/16)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/10080 - (t*(t/153 + 5/153)*(t - 1)*(t + 1)*(t -
13、 2)*(t + 2)*(t - 3)*(t + 3)*(t + 4)*(t - 5)/40320 + (t*(t/260 + 1/52)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)/362880注：Matlab中用t表示变量x通过matlab绘图的插值后多项式曲线为：2切比雪夫零点插值：在-1,1上有11个不同零点： ，=（0.9898 0.9096 0.7557 0.5406 0.2817 0.0000 -0.2817 -0.5406 -0.7557 -0.9096 -0.9898）将-1,1区间到-5,5
14、区间进行转化得：x=5*=（4.9491 4.5482 3.7787 2.7032 1.4087 0.0000 -1.4087 -2.7032 -3.7787 -4.5482 -4.9491）y=1./(1+x.2)=（ 0.0392 0.0461 0.0654 0.1204 0.3351 1.0000 0.3351 0.1204 0.0654 0.0461 0.0392）由于插值后的多项式过于复杂，这里没有给出具体形式，直接给出函数曲线为：原函数、拉格朗日插值、切比雪夫零点插值后的曲线对比如下：3.误差分析：从上图可以看出，在区间两端高次等间距的朗格朗日插值与原函数之间的偏差很大，切比雪夫零
15、点插值多项式在整个区间上与原函数之间的偏差都很小，所以直观上可以看出切比雪夫插值多项式更趋近原函数。下面利用无穷范数对两种插值方法做误差分析，计算结果如下：对比上面计算结果易知用切比雪夫多项式的零点做插值节点得到的多项式于原函数的拟合度更高。4.收获与体会：利用matlab的强大数学计算功能，加深了对拉格朗日插值方法的认识。深刻体会到对于一些函数，等间距节点的高次插值多项式的偏差会很大，通过计算可以看到，利用相应的切比雪夫多项式零点做拉格朗日插值会消弱这一“龙格”现象，而且效果非常明显。题目三：非线性方程求根利用改进Newton法求解方程，其中迭代条件分别为：(1) (2) (3) 对不同条件
16、下获得的结果进行分析比较。一.数学原理：牛顿法解非线性方程的近似解，已知的近似解，通过下式得到更精确地零点近似根。 该题的迭代方式采用的牛顿下山法，是基于牛顿法的改进，即在下山法保证函数值稳定下降的前提下，用牛顿法加快收敛速度，得出下面的迭代计算公式 其中称为下山因子。由数值计算知识可知，虽然(线性方程的重根数)大于1但上式对也是二阶收敛的，所以计算时的初始值取，迭代终止条件为。二.计算过程：（1），；程序运行结果如下：x=0.50175381y=0.50175381 k=10001（k为迭代次数）由于设定了迭代次数最大为10001，因此迭代10001次时还没有使，。（2）程序运行结果如下：x
17、=0.50000000y=0.50000000 k=25可见迭代至25次即找到了最优解。(3) x=0.99999999y=0.99999999 k=12可见迭代至12次即找到了最优解。三结果分析：由方程式可知其根为，。(1)、(2)中的初始值靠近，迭代时向收敛，但（2）中r=1.5可知函数在(2)初始条件下得到的解更精确； (2)迭代25次满足精度要求，可知函数在(2)初始条件下收敛的更快。(3)的初始值为0.85靠近，迭代时向收敛，迭代12次满足精度要求，对比可知在(3)初始条件下，函数的收敛速度最快，得到的解最精确。四收获与体会：牛顿下山法，不仅收敛速度快而且精度高。在所给三个不同初始条
18、件下，迭代的次数和近似解的精确度均不一样，初始重根数越接近1，收敛速度更快。若非线性方程有不同根，给的初始零点值不同时，函数迭代过程会趋向靠近的真实根。数值计算是与计算机密切相关的一门课程，通过书本与matlab相结合很好地理解了其中的原理与方法，为以后学习提供了强有力地工具和方法。附录 程序设计题目一： 解线性方程组的直接法 clear A=1 1 1 1 1 1;1 1.1 1.12 1.13 1.14 1.15;1 1.2 1.22 1.23 1.24 1.25;1 1.3 1.32 1.33 1.34 1.35;1 1.4 1.42 1.43 1.44 1.45;1 1.5 1.52
19、1.53 1.54 1.55A = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.1000 1.2100 1.3310 1.4641 1.6105 1.0000 1.2000 1.4400 1.7280 2.0736 2.4883 1.0000 1.3000 1.6900 2.1970 2.8561 3.7129 1.0000 1.4000 1.9600 2.7440 3.8416 5.3782 1.0000 1.5000 2.2500 3.3750 5.0625 7.5938 X=1;1;1;1;1;1X = 1 1 1 1 1 1 B
20、=A*XB = 6.0000 7.7156 9.9299 12.7560 16.3238 20.7813 Aug=A B n,m=size(Aug)n = 6m = 7 for k = 1:n-1 piv,r = max(abs(Aug(k:n,k); %找列主元所在子矩阵的行r r = r + k - 1; % 列主元所在大矩阵的行 if rk temp=Aug(k,:); Aug(k,:)=Aug(r,:); Aug(r,:)=temp; end if Aug(k,k)=0, error(对角元出现0), end % 把增广矩阵消元成为上三角 for p = k+1:n Aug(p,:)=
21、Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k); end end AugAug = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 6.0000 0 0.5000 1.2500 2.3750 4.0625 6.5938 14.7813 0 0 -0.0600 -0.2220 -0.5514 -1.1492 -1.9826 0 0 0 -0.0080 -0.0408 -0.1306 -0.1794 0 0 0 0 -0.0012 -0.0074 -0.0086 0 0 0 0 0 0.0001 0.0001 % 解上三角方程组 A =
22、Aug(:,1:n); b = Aug(:,n+1); x(n) = b(n)/A(n,n); for k = n-1:-1:1 x(k)=b(k); for p=n:-1:k+1 x(k) = x(k)-A(k,p)*x(p); end x(k)=x(k)/A(k,k); end xx = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000%* %（1)A不变，对B的元素b6加一个扰动10e-4，求解方程组； A=1 1 1 1 1 1;1 1.1 1.12 1.13 1.14 1.15;1 1.2 1.22 1.23 1.24 1.25;1 1.3 1.32
23、 1.33 1.34 1.35;1 1.4 1.42 1.43 1.44 1.45;1 1.5 1.52 1.53 1.54 1.55A = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.1000 1.2100 1.3310 1.4641 1.6105 1.0000 1.2000 1.4400 1.7280 2.0736 2.4883 1.0000 1.3000 1.6900 2.1970 2.8561 3.7129 1.0000 1.4000 1.9600 2.7440 3.8416 5.3782 1.0000 1.5000 2.250
24、0 3.3750 5.0625 7.5938 B =6.0000 7.7156 9.9299 12.7560 16.3238 20.7813+0.0001B = 6.0000 7.7156 9.9299 12.7560 16.3238 20.7814 Aug=A B n,m=size(Aug)n = 6m = 7 for k = 1:n-1 piv,r = max(abs(Aug(k:n,k); %找列主元所在子矩阵的行r r = r + k - 1; % 列主元所在大矩阵的行 if rk temp=Aug(k,:); Aug(k,:)=Aug(r,:); Aug(r,:)=temp; end
25、 if Aug(k,k)=0, error(对角元出现0), end % 把增广矩阵消元成为上三角 for p = k+1:n Aug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k); end end AugAug = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 6.0000 0 0.5000 1.2500 2.3750 4.0625 6.5938 14.7814 0 0 -0.0600 -0.2220 -0.5514 -1.1492 -1.9827 0 0 0 -0.0080 -0.0408 -0.1306 -0.17
26、95 0 0 0 0 -0.0012 -0.0074 -0.0087 0 0 0 0 0 0.0001 0.0001 % 解上三角方程组 A = Aug(:,1:n); b = Aug(:,n+1); x(n) = b(n)/A(n,n); for k = n-1:-1:1 x(k)=b(k); for p=n:-1:k+1 x(k) = x(k)-A(k,p)*x(p); end x(k)=x(k)/A(k,k); end xx = 0.5997 2.6920 -1.8500 3.3917 0.0000 1.1667%* %(2)B不变，对A的元素a22和a66分别加一个扰动10e-6，求解
27、方程组； clear format long A=1 1 1 1 1 1;1 1.1+0.000001 1.12 1.13 1.14 1.15;1 1.2 1.22 1.23 1.24 1.25;1 1.3 1.32 1.33 1.34 1.35;1 1.4 1.42 1.43 1.44 1.45;1 1.5 1.52 1.53 1.54 1.55+0.000001 B =6.000000 7.715600 9.929900 12.756000 16.323800 20.781300B = 6.000000000000000 7.715600000000000 9.92990000000000
28、0 12.756000000000000 16.323799999999999 20.781300000000002 Aug=A B n,m=size(Aug); for k = 1:n-1 piv,r = max(abs(Aug(k:n,k); %找列主元所在子矩阵的行r r = r + k - 1; % 列主元所在大矩阵的行 if rk temp=Aug(k,:); Aug(k,:)=Aug(r,:); Aug(r,:)=temp; end if Aug(k,k)=0, error(对角元出现0), end % 把增广矩阵消元成为上三角 for p = k+1:n Aug(p,:)=Aug
29、(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k); end % 解上三角方程组 A = Aug(:,1:n); b = Aug(:,n+1); x(n) = b(n)/A(n,n); for k = n-1:-1:1 x(k)=b(k); for p=n:-1:k+1 x(k) = x(k)-A(k,p)*x(p); end x(k)=x(k)/A(k,k); end xx = 0.825832305593523 1.742109746543727 -0.259537880024995 2.064585114452380 0.551832152872303 1.07517856
30、0563062题目二：多项式插值在区间上对龙格函数做插值，分别用等距节点(节点步长)和切比雪夫多项式的零点做插值节点，画出原函数和两个插值函数的图像进行比较，并利用对两种插值方法做误差分析。原函数： x=-5:0.1:5; y=1./(1+x.2); plot(x,y)2. 利用拉格朗日差分(等距节点)： x=-5:1:5x = -5 -4 -3 -2 -1 0 1 2 3 4 5 y=1./(1+x.2)y =0.0385 0.0588 0.1000 0.2000 0.5000 1.0000 0.5000 0.2000 0.1000 0.0588 0.0385 syms t l;if(len
31、gth(x) = length(y) n = length(x);else disp(x和y的维数不相等！); return; %检错end p=sym(0);for (i=1:n) l=sym(y(i); for(k=1:i-1) l=l*(t-x(k)/(x(i)-x(k); end; for(k=i+1:n) l=l*(t-x(k)/(x(i)-x(k); end; p=p+l;end pp =(t*(t/8 + 5/8)*(t - 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/4320 - (t*(t/15 + 1
32、/3)*(t - 1)*(t + 1)*(t - 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/10080 - (t/5 + 1)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/2880 + (t*(t/20 + 1/4)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t - 4)*(t + 4)*(t - 5)/40320 + (t*(t/12 + 5/12)*(t + 1)*(t - 2)*(t + 2)*(t - 3
33、)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/2880 - (t*(t/17 + 5/17)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t - 5)/362880 + (t*(t/26 + 2/13)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t - 5)/3628800 - (t*(t/35 + 1/7)*(t - 1)*(t + 1)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)*(t - 5
34、)/4320 + (t*(t/80 + 1/16)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t + 3)*(t - 4)*(t + 4)*(t - 5)/10080 - (t*(t/153 + 5/153)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t + 4)*(t - 5)/40320 + (t*(t/260 + 1/52)*(t - 1)*(t + 1)*(t - 2)*(t + 2)*(t - 3)*(t + 3)*(t - 4)*(t + 4)/362880 x0=-5:0.1:5; f = subs
35、(p,t,x0); %计算插值点的函数值 plot(x0,f)2切比雪夫零点插值：在-1,1上有11个不同零点： ， k=1:11k = 1 2 3 4 5 6 7 8 9 10 11 s=cos(2.*k-1).*pi./22)s =0.9898 0.9096 0.7557 0.5406 0.2817 0.0000 -0.2817 -0.5406 -0.7557 -0.9096 -0.9898 x=5*sx = 4.9491 4.5482 3.7787 2.7032 1.4087 0.0000 -1.4087 -2.7032 -3.7787 -4.5482 -4.9491 y=1./(1+x
36、.2)y = 0.0392 0.0461 0.0654 0.1204 0.3351 1.0000 0.3351 0.1204 0.0654 0.0461 0.0392 syms t l;if(length(x) = length(y) n = length(x);else disp(x和y的维数不相等！); return; %检错endp=sym(0);for (i=1:n) l=sym(y(i); for(k=1:i-1) l=l*(t-x(k)/(x(i)-x(k); end; for(k=i+1:n) l=l*(t-x(k)/(x(i)-x(k); end; p=p+l;end x0=-
37、5:0.1:5; f = subs (p,t,x0); %计算插值点的函数值plot(x0,f) 题目三：非线性方程求根 (1) clear;clc;x=0.55;r=2;k=1;while k10e-9;x=y;else breakendk=k+1;endfprintf(n%s%.8ft%s%.8f t%s%d,x=,x, y=,y,k=,k)x=0.50175381y=0.50175381 k=10001f=2*y3-5*y2+4*y-1f = 8.7076e-004(2) clear;clc;x=0.55;r=1.5;k=1;while k10e-9;x=y;else breakendk=k+1;endfprintf(n%s%.8ft%s%.8f t%s%d,x=,x, y=,y,k=,k)x=0.50000000y=0.50000000 k=25 f=2*y3-5*y2+4*y-1f = -8.5649e-010(3) clear;clc;x=0.85;r=1.5;k=1;while k10e-9;x=y;else breakendk=k+1;endfprintf(n%s%.8ft%s%.8f t%s%d,x=,x, y=,y,k=,k)x=0.99999999y=0.99999999 k=12 f=2*y3-