数值分析大作业(牛顿下山法拉格朗日法切比雪夫法)及Matlab程序.pdf

课程名称：设计题目：学号：姓名：完成时间：程设计数值分析2014.11.18课题目一：题目一：解线性方程组的直接法解线性方程组的直接法设方程组Axb，其中2x0 x125x0 x15，5x51x01x1A1x5矩阵中xk10.1k(k 0,1,x (1,1,1)T。2x5,5)，b由相应的矩阵元素计算，使解向量(1)A不变，对b的元素b6加一个扰动104，求解方程组；(2)b不变，对A的元素a22和a66分别加一个扰动106，求解方程组；(3)对上述两种扰动方程组的解做误差分析。一.数学原理：本计算采用直接法中的列主元高斯消元法，高斯列主元消元法原理如下：1、设有 n 元线性方程组如下：a11an1a1nannx1xnb1bn2、第一步：如果 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)n-1(1)(1)(1)(1)an1 an2 .ann xn bn简记为：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.=.0 0 .a(n-1)n-1n xn-1 b(n-1)n-10 0 .a(n)nn xn b(n)n简记为：(n)(n)Ax=b最后从方程组的最后一个方程进行回代求解为：Xn=b(n)/a(n)nn Xi=(b(k)k-a(k)kjxj)/a(k)kk二解题过程：1.由题中所给条件可求出 b。B=6.00007.71569.929912.756016.323820.7813（1）A不变，对b的元素b6加一个扰动10，求解方程组。4B=6.0000 7.7156 9.9299 12.7560 16.3238 20.7813+0.0001解得 x=0.59972.6920-1.85003.39170.00001.1667（2）b不变，对A的元素a22和a66分别加一个的扰动，求解方程组。A=1.0000000000000001.0000000000000001.0000000000000001.0000000000000001.0000000000000001.0000000000000001.0000000000000001.1000010000000001.2100000000000001.3310000000000001.4641000000000011.6105100000000011.0000000000000001.2000000000000001.4400000000000001.7280000000000002.0736000000000002.4883199999999991.0000000000000001.3000000000000001.6900000000000002.1970000000000002.8561000000000013.7129300000000011.0000000000000001.4000000000000001.9600000000000002.7439999999999993.8415999999999995.3782399999999981.0000000000000001.5000000000000002.2500000000000003.3750000000000005.0625000000000007.593751000000000 x=0.8258323055935231.742109746543727-0.2595378800249952.0645851144523800.5518321528723031.075178560563062三、误差分析：从上面计算结果可以看出，当系数矩阵或右端向量发生极小的扰动，方程组的解也会产生很大的误差，产生的原因是范德蒙阵为变态阵。由数值计算知识可知xxb condAA1b1 AAA其中condA A A1为条件数，从上式看到，当 A 的条件数很大时，解的相对误差也很大，此时的对应的线性方程为病态线性方程组。计算条件数时，取矩阵的无穷范数，经计算得矩阵A、受c,d扰动后的矩阵D和E的条件数为condA1.0668e+07；condB6.9639e+06；condE 2.3434e+07；可以看到三个矩阵的条件数非常大，即使系数矩阵或右端向量发生很小的变化，也会导致解产生很大的误差。四、收获与体会：运用matlab编程解决数学问题很方便，病态阵的条件数非常大，给系数矩阵或者右端向量一个很微小的扰动，方程组的解也会产生很大的变动。通过做这个题目，我对让课上抽象病态现象有了直观的认识。题目二：多项式插值题目二：多项式插值在区间5,5上对龙格函数R(x)1做插值，分别用等距节点(节点步长h1)21 x和切比雪夫多项式T11(x)的零点做插值节点，画出原函数和两个插值函数的图像进行比较，并利用|对两种插值方法做误差分析。一.数学原理：1.拉格朗日等距节点插值拉格朗日插值多项式为其中余项为Lnxlkxyk=ykk0k0nnn1xx xkn1xn1x=xkx0 xkxk1xkxk1 xkxnfn1Rnx fxLnxxn,x有关。n1!n1x与x02.拉格朗日以切比雪夫零点为节点插值切比雪夫多项式Tnxcosnarccosx,x 1.在区间1,1上有n个零点，为xk cos2k 1,k 0,1,2,32n,n1由上式得到切比雪夫多项式Tn+1x的n1个零点，并作为插值节点做拉格朗日插值，得到多项式Pnx。一.计算过程：原函数为 y=1./(1+x.2)，通过 matlab 绘图得到原函数图像为：1.利用拉格朗日差分(等距节点)：x=(-5-4-3-2-1012345)y=1./(1+x.2)=(0.03850.05880.10000.20000.50001.00000.50000.20000.10000.05880.0385)拉格朗日插值多项式为：l10(x)y0l0(x)y1l1(x)y2l2(x)y3l3(x)y4l4(x)y5l5(x)y6l6(x)y7l7(x)y8l8(x)y9l9(x)y10l10(x)=(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-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/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-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切比雪夫零点插值：T11(x)在-1,1上有 11 个不同零点：xk cos(2k 1)(k 1,2,，11)22xk=（0.98980.90960.75570.54060.28170.0000-0.2817-0.5406-0.7557-0.9096-0.9898）将-1,1区间到-5,5区间进行转化得：x=5*xk=（4.94914.54823.77872.70321.40870.0000-1.4087-2.7032-3.7787-4.5482-4.9491）y=1./(1+x.2)=（0.03920.04610.06540.12040.33511.00000.33510.12040.06540.04610.0392）由于插值后的多项式过于复杂，这里没有给出具体形式，直接给出函数曲线为：原函数、拉格朗日插值、切比雪夫零点插值后的曲线对比如下：3.误差分析：从上图可以看出，在区间两端高次等间距的朗格朗日插值与原函数之间的偏差很大，切比雪夫零点插值多项式在整个区间上与原函数之间的偏差都很小，所以直观上可以看出切比雪夫插值多项式更趋近原函数。下面利用无穷范数对两种插值方法做误差分析，计算结果如下：fxP11xfx L11x max fxP11x=0.10925x5 max fxL11x=1.91565x5对比上面计算结果易知用切比雪夫多项式T11x的零点做插值节点得到的多项式于原函数的拟合度更高。4.收获与体会：利用 matlab 的强大数学计算功能，加深了对拉格朗日插值方法的认识。深刻体会到对于一些函数，等间距节点的高次插值多项式的偏差会很大，通过计算可以看到，利用相应的切比雪夫多项式零点做拉格朗日插值会消弱这一“龙格”现象，而且效果非常明显。题目三：非线性方程求根题目三：非线性方程求根利用改进Newton法求解方程(x1)2(2x1)0，其中迭代条件分别为：(1)r 2,(2)r 1.5,x0 0.55;x0 0.55;(3)r 1.5,x0 0.85;对不同条件下获得的结果进行分析比较。一.数学原理：牛顿法解非线性方程的近似解，已知fx0的近似解xk，通过下式得到fx更精确地零点近似根。xk+1=xkfxk fxk该题的迭代方式采用的牛顿下山法，是基于牛顿法的改进，即在下山法保证函数值稳定下降的前提下，用牛顿法加快收敛速度，得出下面的迭代计算公式xk+1=xkfxk,k 0,1,2,fxk.其中01称为下山因子。由数值计算知识可知，虽然r(线性方程的重根数)大于 1 但上式对r也是二阶收敛的，所以计算时的初始值取r，迭代终止条件为xk1xk109。二.计算过程：（1）r 2，x0 0.55；程序运行结果如下：x=0.50175381y=0.50175381k=10001（k 为迭代次数）由 于 设 定 了 迭 代 次 数 最 大 为 10001，因 此 迭 代 10001 次时 还 没 有 使xk1xk109（2）r 1.5,x0 0.55;，。程序运行结果如下：x=0.50000000y=0.50000000k=25可见迭代至 25 次即找到了最优解。(3)r 1.5,x0 0.85;x=0.99999999y=0.99999999k=12可见迭代至 12 次即找到了最优解。三结果分析：由方程式可知其根为x11，x2 0.5。(1)、(2)中的初始值x0靠近x2，迭代时向x2收敛，但（2）中 r=1.5 可知函数在(2)初始条件下得到的解更精确；(2)迭代 25 次满足精度要求，可知函数在(2)初始条件下收敛的更快。(3)的初始值为0.85 靠近x1，迭代时向x1收敛，迭代 12 次满足精度要求，对比可知在(3)初始条件下，函数的收敛速度最快，得到的解最精确。四收获与体会：牛顿下山法，不仅收敛速度快而且精度高。在所给三个不同初始条件下，迭代的次数和近似解的精确度均不一样，初始重根数越接近 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.331.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.00001.00001.00001.00001.00001.00001.00001.10001.21001.33101.46411.61051.00001.20001.44001.72802.07362.48831.00001.30001.69002.19702.85613.71291.00001.40001.96002.74403.84165.37821.00001.50002.25003.37505.06257.5938 X=1;1;1;1;1;1X=111111 B=A*XB=6.00007.71569.929912.756016.323820.7813 Aug=A B n,m=size(Aug)n=6m=7 for k=1:n-1piv,r=max(abs(Aug(k:n,k);%找列主元所在子矩阵的行rr=r+k-1;%列主元所在大矩阵的行if rktemp=Aug(k,:);Aug(k,:)=Aug(r,:);Aug(r,:)=temp;endif Aug(k,k)=0,error(对角元出现 0),end%把增广矩阵消元成为上三角for p=k+1:nAug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k);endend AugAug=1.00001.00001.00001.00001.00001.00006.000000.50001.25002.37504.06256.593814.781300-0.0600-0.2220-0.5514-1.1492-1.9826000-0.0080-0.0408-0.1306-0.17940000-0.0012-0.0074-0.0086000000.00010.0001%解上三角方程组A=Aug(:,1:n);b=Aug(:,n+1);x(n)=b(n)/A(n,n);for k=n-1:-1:1x(k)=b(k);for p=n:-1:k+1x(k)=x(k)-A(k,p)*x(p);endx(k)=x(k)/A(k,k);end xx=1.00001.00001.00001.00001.00001.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 1.331.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.00001.00001.00001.00001.00001.00001.00001.10001.21001.33101.46411.61051.00001.20001.44001.72802.07362.48831.00001.30001.69002.19702.85613.71291.00001.40001.96002.74403.84165.37821.00001.50002.25003.37505.06257.5938 B=6.0000 7.7156 9.9299 12.7560 16.3238 20.7813+0.0001B=6.00007.71569.929912.756016.323820.7814Aug=A B n,m=size(Aug)n=6m=7 for k=1:n-1piv,r=max(abs(Aug(k:n,k);%找列主元所在子矩阵的行rr=r+k-1;%列主元所在大矩阵的行if rktemp=Aug(k,:);Aug(k,:)=Aug(r,:);Aug(r,:)=temp;endif Aug(k,k)=0,error(对角元出现 0),end%把增广矩阵消元成为上三角for p=k+1:nAug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k);endendAugAug=1.00001.00001.00001.00001.00001.00006.000000.50001.25002.37504.06256.593814.781400-0.0600-0.2220-0.5514-1.1492-1.9827000-0.0080-0.0408-0.1306-0.17950000-0.0012-0.0074-0.0087000000.00010.0001%解上三角方程组A=Aug(:,1:n);b=Aug(:,n+1);x(n)=b(n)/A(n,n);for k=n-1:-1:1x(k)=b(k);for p=n:-1:k+1x(k)=x(k)-A(k,p)*x(p);endx(k)=x(k)/A(k,k);end xx=0.59972.6920-1.85003.39170.00001.1667%*%(2)B 不变，对 A 的元素 a22 和 a66 分别加一个扰动 10e-6，求解方程组；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.31.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.0000000000000007.7156000000000009.92990000000000012.75600000000000016.32379999999999920.781300000000002 Aug=A B n,m=size(Aug);for k=1:n-1piv,r=max(abs(Aug(k:n,k);%找列主元所在子矩阵的行rr=r+k-1;%列主元所在大矩阵的行if rktemp=Aug(k,:);Aug(k,:)=Aug(r,:);Aug(r,:)=temp;endif Aug(k,k)=0,error(对角元出现 0),end%把增广矩阵消元成为上三角for p=k+1:nAug(p,:)=Aug(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:1x(k)=b(k);for p=n:-1:k+1x(k)=x(k)-A(k,p)*x(p);endx(k)=x(k)/A(k,k);end xx=0.8258323055935231.742109746543727-0.2595378800249952.0645851144523800.5518321528723031.075178560563062题目二：多项式插值题目二：多项式插值在区间5,5上对龙格函数R(x)1做插值，分别用等距节点(节点步长h 1)1 x2和切比雪夫多项式T11(x)的零点做插值节点，画出原函数和两个插值函数的图像进行比较，并利用|对两种插值方法做误差分析。原函数：x=-5:0.1:5;y=1./(1+x.2);plot(x,y)2.利用拉格朗日差分(等距节点)：x=-5:1:5x=-5-4-3-2-1012345 y=1./(1+x.2)y=0.03850.05880.10000.20000.50001.00000.50000.20000.10000.05880.0385 syms t l;if(length(x)=length(y)n=length(x);elsedisp(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/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)*(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)/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(p,t,x0);%计算插值点的函数值 plot(x0,f)2切比雪夫零点插值：T(2k 1)11(x)在-1,1上有 11 个不同零点：xk cos22(k 1,2,，11)k=1:11k=1234567891011s=cos(2.*k-1).*pi./22)s=0.98980.90960.75570.54060.28170.0000-0.2817-0.7557-0.9096-0.9898 x=5*sx=4.94914.54823.77872.70321.40870.0000-2.7032-3.7787-4.5482-4.9491 y=1./(1+x.2)y=0.03920.04610.06540.12040.33511.00000.12040.06540.04610.0392-0.5406-1.40870.3351syms t l;if(length(x)=length(y)n=length(x);elsedisp(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=-5:0.1:5;f=subs(p,t,x0);%计算插值点的函数值plot(x0,f)题目三：非线性方程求根题目三：非线性方程求根(1)r 2,clear;clc;x=0.55;r=2;x0 0.55;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.50175381k=10001f=2*y3-5*y2+4*y-1f=8.7076e-004(2)r 1.5,x0 0.55;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.50000000k=25 f=2*y3-5*y2+4*y-1f=-8.5649e-010(3)r 1.5,clear;clc;x=0.85;r=1.5;k=1;x0 0.85;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.99999999f=2*y3-5*y2+4*y-1f=0k=12