数值分析大作业（牛顿下山法，拉格朗日法，切比雪夫法）及Matlab程序

课程名称：设计题目：

学 号：姓 名：

完成时间：

程 设 计

数值分析

2014.11.18

课

题目一： 解线性方程组的直接法 设方程组Ax?b，其中

?1x0?1x1A?????1x5?2x05?x0?x15?， ?5?x5??x122x5矩阵中xk?1?0.1k(k?0,1,,5)，b由相应的矩阵元素计算，使解向量

x?(1,1,,1)T。

(1) A不变，对b的元素b6加一个扰动10?4，求解方程组；

(2) b不变，对A的元素a22和a66分别加一个扰动10?6，求解方程组； (3) 对上述两种扰动方程组的解做误差分析。

一.数学原理：

本计算采用直接法中的列主元高斯消元法，高斯列主元消元法原理如下： 1、设有n元线性方程组如下：

?a11???a?n1a1n???ann???x1???x?n??b1???＝???b??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)1 a(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)1

0 a(1)22 . . . a(1)2n x2 b(1)2 . . . . . . . = . 0 0 . . a(n-1)n-1n xn-1 b(n-1)n-1 0 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.0000 7.7156 9.9299 12.7560 16.3238 20.7813

（1） A不变，对b的元素b6加一个扰动10，求解方程组。

?4B =[6.0000 7.7156 9.9299 12.7560 16.3238 20.7813+0.0001]'

解得x =[0.5997 2.6920 -1.8500 3.3917 0.0000 1.1667]’

（2）b不变，对A的元素a22和a66分别加一个的扰动，求解方程组。

A =

1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000

1.000000000000000 1.100001000000000 1.210000000000000 1.331000000000000 1.464100000000001 1.610510000000001

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.841599999999999 5.378239999999998

1.000000000000000 1.500000000000000 2.250000000000000 3.375000000000000 5.062500000000000 7.593751000000000

x =[0.825832305593523 1.742109746543727 -0.259537880024995 2.064585114452380 0.551832152872303 1.075178560563062]’

三、误差分析：

从上面计算结果可以看出，当系数矩阵或右端向量发生极小的扰动，方程组的解也会产生很大的误差，产生的原因是范德蒙阵为变态阵。

由数值计算知识可知

?xx??b?cond?A???A? ???1??b?1?A?A?A其中cond?A??A?A?1 为条件数，从上式看到，当A的条件数很大时，解的相对误差也很大，此时的对应的线性方程为病态线性方程组。

计算条件数时，取矩阵的无穷范数，经计算得矩阵A、受c,d扰动后的矩阵D和E的条件数为

cond?A??1.0668e+07；cond?B??6.9639e+06；cond?E??2.3434e+07； 可以看到三个矩阵的条件数非常大，即使系数矩阵或右端向量发生很小的变化，也会导致解产生很大的误差。

四、收获与体会：

运用matlab编程解决数学问题很方便，病态阵的条件数非常大，给系数矩阵或者右端向量一个很微小的扰动，方程组的解也会产生很大的变动。通过做这个题目，我对让课上抽象病态现象有了直观的认识。

题目二：多项式插值

在区间[?5,5]上对龙格函数R(x)?1做插值，分别用等距节点(节点步长h?1)1?x2和切比雪夫多项式T11(x)的零点做插值节点，画出原函数和两个插值函数的图像进行比较，并利用||?||?对两种插值方法做误差分析。 一.数学原理： 1.拉格朗日等距节点插值 拉格朗日插值多项式为

其中 余项为

Ln?x???lk?x?yk=?ykk?0k?0nn?n?1?x??x?xk??n?1??x?

?n?1?x?=?xk?x0??xk?xk?1??xk?xk?1??xk?xn?

f?n?1????Rn?x??f?x??Ln?x???n?1?x?

?n?1?!?与x0xn,x有关。

2.拉格朗日以切比雪夫零点为节点插值 切比雪夫多项式

Tn?x??cos?narccosx?,x?1.

在区间??1,1? 上有n 个零点，为

xk?cos2k?1?,k?0,1,2,32n,n?1

由上式得到切比雪夫多项式Tn+1?x?的n?1个零点，并作为插值节点做拉格朗日插值，得到多项式Pn?x?。 一.计算过程：

原函数为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.2000 0.5000 1.0000 0.5000 0.2000 0.1000 0.0588 0.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.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]区间进行转化得：

x=5*xk=（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.误差分析：

从上图可以看出，在区间两端高次等间距的朗格朗日插值与原函数之间的偏差很大，切比雪夫零点插值多项式在整个区间上与原函数之间的偏差都很小，所以直观上可以看出切比雪夫插值多项式更趋近原函数。下面利用无穷范数对两种插值方法做误差分析，计算结果如下：

f?x??P11?x?f?x??L11?x???maxf?x??P11?x?=0.1092

?5?x5??maxf?x??L11?x?=1.9156

?5?x5对比上面计算结果易知用切比雪夫多项式T11?x?的零点做插值节点得到的多项式于原函数的拟合度更高。 4.收获与体会：

利用matlab的强大数学计算功能，加深了对拉格朗日插值方法的认识。深刻体会到对于一些函数，等间距节点的高次插值多项式的偏差会很大，通过计算可以看到，利用相应的切比雪夫多项式零点做拉格朗日插值会消弱这一“龙格”现象，而且效果非常明显。 题目三：非线性方程求根

利用改进Newton法求解方程(x?1)2(2x?1)?0，其中迭代条件分别为： (1) r?2,(2) r?1.5,x0?0.55;

x0?0.55;

(3) r?1.5,x0?0.85;

对不同条件下获得的结果进行分析比较。 一.数学原理：

牛顿法解非线性方程的近似解，已知f?x??0的近似解xk，通过下式得到

f?x?更精确地零点近似根。

xk+1=xk?f?xk? ?f?xk?该题的迭代方式采用的牛顿下山法，是基于牛顿法的改进，即在下山法保证函数值稳定下降的前提下，用牛顿法加快收敛速度，得出下面的迭代计算公式

xk+1=xk??f?xk?,k?0,1,2,?f?xk?.

其中??0＜??1?称为下山因子。由数值计算知识可知，虽然r(线性方程的重根数)大于1但上式对r也是二阶收敛的，所以计算时?的初始值取r，迭代终止条件为xk?1?xk?10?9。 二.计算过程：

（1）r?2，x0?0.55；

程序运行结果如下：

x=0.50175381 y=0.50175381 k=10001（k为迭代次数）

由于设定了迭代次数最大为10001，因此迭代10001次时还没有使

xk?1?xk?10?9（2）r?1.5,，。。

x0?0.55;

程序运行结果如下：

x=0.50000000 y=0.50000000 k=25 可见迭代至25次即找到了最优解。 (3) r?1.5,x0?0.85;

x=0.99999999 y=0.99999999 k=12 可见迭代至12次即找到了最优解。 三．结果分析：

由方程式可知其根为x1??1，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.1^2 1.1^3 1.1^4 1.1^5;1 1.2 1.2^2 1.2^3 1.2^4 1.2^5;1 1.3 1.3^2 1.3^3 1.3^4 1.3^5;1 1.4 1.4^2 1.4^3 1.4^4 1.4^5;1 1.5 1.5^2 1.5^3 1.5^4 1.5^5] A =

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;1]’ X = 1 1 1 1 1 1 >> B=A*X B =

6.0000 7.7156 9.9299 12.7560 16.3238 20.7813 >> Aug=[A B]

>> [n,m]=size(Aug) n =

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:11 k =

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.7557 -0.9096 -0.9898

>> x=5*s x =

4.9491 4.5482 3.7787 2.7032 1.4087 0.0000 -2.7032 -3.7787 -4.5482 -4.9491

>> y=1./(1+x.^2) y =

0.0392 0.0461 0.0654 0.1204 0.3351 1.0000 0.1204 0.0654 0.0461 0.0392

-0.5406 -1.4087 0.3351 >> syms t l;

if(length(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

>> 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 k<=10000;

y=x-r*(2*x^3-5*x^2+4*x-1)/(6*x^2-10*x+4); if abs(y-x)>10e-9; x=y;

else break end k=k+1; end

fprintf('\\n%s%.8f\\t%s%.8f \\t%s%d','x=',x, 'y=',y,'k=',k) x=0.50175381 y=0.50175381 k=10001>> f=2*y^3-5*y^2+4*y-1 f =

8.7076e-004 (2) r?1.5,x0?0.55;

clear; clc; x=0.55; r=1.5; k=1;

while k<=10000;

y=x-r*(2*x^3-5*x^2+4*x-1)/(6*x^2-10*x+4); if abs(y-x)>10e-9; x=y;

else break end k=k+1; end

fprintf('\\n%s%.8f\\t%s%.8f \\t%s%d','x=',x, 'y=',y,'k=',k) x=0.50000000 y=0.50000000 k=25>> >> f=2*y^3-5*y^2+4*y-1 f =

-8.5649e-010 (3) r?1.5,clear; clc; x=0.85; r=1.5; k=1;

x0?0.85;

while k<=10000;

y=x-r*(2*x^3-5*x^2+4*x-1)/(6*x^2-10*x+4); if abs(y-x)>10e-9; x=y;

else break end k=k+1; end

fprintf('\\n%s%.8f\\t%s%.8f \\t%s%d','x=',x, 'y=',y,'k=',k) x=0.99999999 y=0.99999999 >> f=2*y^3-5*y^2+4*y-1 f =

0

k=12>>

本文来源：https://www.bwwdw.com/article/31fd.html