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

课程设计
课程名称：数值分析
设计题目：
学号：
姓名：
完成时间：2014.11.18
题目一： 解线性方程组的直接法 设方程组Ax b =，其中
250
002
51
112
55
5
5111x x x x x x A x x x ⎡⎤⎢
⎥⎢
⎥=⎢⎥⎢⎥⎢⎥⎣⎦
， 矩阵中10.1(0,1,
,5)k x k k =+=，b 由相应的矩阵元素计算，使解向量
(1,1,
,1)T x =。

(1) A 不变，对b 的元素6b 加一个扰动410-，求解方程组；
(2) b 不变，对A 的元素22a 和66a 分别加一个扰动610-，求解方程组； (3) 对上述两种扰动方程组的解做误差分析。

一.数学原理：
本计算采用直接法中的列主元高斯消元法，高斯列主元消元法原理如下： 1、设有n 元线性方程组如下：
1111n n nn a a a a ⎛⎫ ⎪ ⎪ ⎪⎝⎭
1
n
x x ⎛⎫ ⎪ ⎪ ⎪⎝⎭＝1n
b b ⎛⎫ ⎪ ⎪ ⎪⎝⎭
2、
第一步：如果a11!=0, 令
l i1= ai1/a11, I= 2,3,……,n
用（-li1）乘第一个方程加到第i 个方程上,得同解方程组：
a (1)11 a (1)12 . . . a (1)1n
x 1 b (1)1 a (1)21 a (1)22 . . . a (1)2n x 2 b (1)2 . . . . . . . = . a (1)n-11 a (1)n-12 . . a (1)n-1n x n-1 b (1)n-1 a (1)n1 a (1)n2 . . . a (1)nn x n b (1)n
简记为：
A (2) x = b (2) 其中
a (2)ij = a (1)ij – l i1 * a (1)1j , I ,j = 2,3,..,n
b (2)I = b (1)I – l i1 * b (1)1 , I = 2,3,...,n 第二步：如果a (2)22 != 0,令
l i2= a （2）i2/a （2）22, I= 3,……,n
依据同样的原理，对矩阵进行化间（省略），依次下去，直到完成！最后，得到上三角方程组：
a(1)
11 a(1)
12
. . . a(1)
1n
x
1
b(1)
1
0 a(1)
22 . . . a(1)
2n
x
2
b(1)
2
. . . . . . . = .
0 0 . . a(n-1)
n-1n x
n-1
b(n-1)
n-1
0 0 . . . a(n)
nn x
n
b(n)
n
简记为：
A(n) x = b(n)
最后从方程组的最后一个方程进行回代求解为：
X
n = b(n) / a(n)
nn
X
i = ( b(k)
k
- ∑ a(k)kj x j ) / a(k)kk
二．解题过程：
1.由题中所给条件可求出b。

B =
6.0000
7.7156
9.9299
12.7560
16.3238
20.7813
（1）A不变，对b的元素6b加一个扰动4
10-，求解方程组。

B =[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的元素22a和66a分别加一个的扰动，求解方程组。

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]’
三、误差分析：
从上面计算结果可以看出，当系数矩阵或右端向量发生极小的扰动，方程组的解也会产生很大的误差，产生的原因是范德蒙阵为变态阵。

由数值计算知识可知
()1cond 1x
A
b A x
A
b A A δδδδ-⎛⎫≤
+ ⎪ ⎪-⎝⎭ 其中()1cond A A A -=⋅ 为条件数，从上式看到，当A 的条件数很大时，解的相对误差也很大，此时的对应的线性方程为病态线性方程组。

计算条件数时，取矩阵的无穷范数，经计算得矩阵A 、受,c d 扰动后的矩阵D 和E 的条件数为
()cond 1.0668e+07A =；()cond 6.9639e+06B =；()cond 2.3434e+07E =； 可以看到三个矩阵的条件数非常大，即使系数矩阵或右端向量发生很小的变化，也会导致解产生很大的误差。

四、收获与体会：
运用matlab 编程解决数学问题很方便，病态阵的条件数非常大，给系数矩阵或者右端向量一个很微小的扰动，方程组的解也会产生很大的变动。

通过做这个题目，我对让课上抽象病态现象有了直观的认识。

题目二：多项式插值
在区间[5,5]-上对龙格函数2
1
()1R x x
=
+做插值，分别用等距节点(节点步长1h =)和切比雪夫多项式11()T x 的零点做插值节点，画出原函数和两个插值函数的图像进行比较，并利用||||∞⋅对两种插值方法做误差分析。

一.数学原理：
1.拉格朗日等距节点插值 拉格朗日插值多项式为
()()()
()()
10
1=n n
n n k k k
k k k n x L x l x y y x x x ωω+==+='-∑∑
其中 ()()()()()1011=n k k k k k k n x x x x x x x x x ω+-+----
余项为
()()()()()
()()111!
n n n n f R x f x L x x n ξω++=-=
+
ξ与0
,n x x x 有关。

2.拉格朗日以切比雪夫零点为节点插值 切比雪夫多项式
()()cos arccos , 1.n T x n x x =≤
在区间[]1,1- 上有n 个零点，为
21
cos
,0,1,2,3,12k k x k n n
π+==-
由上式得到切比雪夫多项式()+1n T x 的1n +个零点，并作为插值节点做拉格朗日插值，得到多项式()n P 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)
拉格朗日插值多项式为：
10001122334455667788991010()()()()()()()()()()()()
l x y l x y l x y l x y l x y l x y l x y l x y l x y l x y l x y l 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．切比雪夫零点插值：
11()
T x在[-1,1]上有11个不同零点：
(21)
cos
22
k
k
x
π
-
=(1,2,
k= (11)
k
x=（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*
k
x=（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.误差分析：
从上图可以看出，在区间两端高次等间距的朗格朗日插值与原函数之间的偏差很大，切比雪夫零点插值多项式在整个区间上与原函数之间的偏差都很小，所以直观上可以看出切比雪夫插值多项式更趋近原函数。

下面利用无穷范数对两种插值方法做误差分析，计算结果如下：
()()()()111155
max =0.1092x f x P x f x P x ∞
-≤-=-
()()
()()111155
max =1.9156x f x L x f x L x ∞
-≤-=-
对比上面计算结果易知用切比雪夫多项式()11T x 的零点做插值节点得到的多项式于原函数的拟合度更高。

4.收获与体会：
利用matlab 的强大数学计算功能，加深了对拉格朗日插值方法的认识。

深刻体会到对于一些函数，等间距节点的高次插值多项式的偏差会很大，通过计算可以看到，利用相应的切比雪夫多项式零点做拉格朗日插值会消弱这一“龙格”现象，而且效果非常明显。

题目三：非线性方程求根
利用改进Newton 法求解方程2(1)(21)0x x --=，其中迭代条件分别为： (1) 02,0.55;r x ==
(2) 01.5,
0.55;r x ==
(3) 01.5,
0.85;r x ==
对不同条件下获得的结果进行分析比较。

一.数学原理：
牛顿法解非线性方程的近似解，已知()0f x =的近似解k x ，通过下式得到
()f x 更精确地零点近似根。

()
()
+1=k k k k f x x x f x -
' 该题的迭代方式采用的牛顿下山法，是基于牛顿法的改进，即在下山法保证函数值稳定下降的前提下，用牛顿法加快收敛速度，得出下面的迭代计算公式
()
()
+1=,0,1,2,.k k k k f x x x k f x λ
-='
其中()01λλ≤＜称为下山因子。

由数值计算知识可知，虽然r (线性方程的重根数)大于1但上式对r 也是二阶收敛的，所以计算时λ的初始值取r ，迭代终止条件为9110k k x x -+-≤。

二.计算过程：
（1）2r =，00.55x =；
程序运行结果如下：
x=0.50175381 y=0.50175381 k=10001（k 为迭代次数）
由于设定了迭代次数最大为10001，因此迭代10001次时还没有使
9
110k k x x -+-≤，。

（2）01.5,
0.55;
r x ==
程序运行结果如下：
x=0.50000000 y=0.50000000 k=25 可见迭代至25次即找到了最优解。

(3) 01.5,
0.85;r x ==
x=0.99999999 y=0.99999999 k=12 可见迭代至12次即找到了最优解。

三．结果分析：
由方程式可知其根为11x *=，20.5x *=。

(1)、(2)中的初始值0x 靠近2x *，迭代时向2x *收敛，但（2）中r=1.5可知函数在(2)初始条件下得到的解更精确； (2)迭代25次满足精度要求，可知函数在(2)初始条件下收敛的更快。

(3)的初始值为0.85靠近1x *，迭代时向1x *收敛，迭代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 =
6
m =
7
>> for k = 1:n-1
[piv,r] = max(abs(Aug(k:n,k))); %找列主元所在子矩阵的行r
r = r + k - 1; % 列主元所在大矩阵的行
if r>k
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(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k);
end
end
>> Aug
Aug =
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 = 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
>> x
x =
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.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
>> B =[6.0000 7.7156 9.9299 12.7560 16.3238 20.7813+0.0001]'
B =
6.0000
7.7156
9.9299
12.7560
16.3238
20.7814
>> Aug=[A B]
>> [n,m]=size(Aug)
n =
6
m =
7
>> for k = 1:n-1
[piv,r] = max(abs(Aug(k:n,k))); %找列主元所在子矩阵的行r
r = r + k - 1; % 列主元所在大矩阵的行
if r>k
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(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k);
end
end
>> Aug
Aug =
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.1795
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
>> x
x =
0.5997 2.6920 -1.8500 3.3917 0.0000 1.1667
%***************************************************************************** >> %(2)B不变，对A的元素a22和a66分别加一个扰动10e-6，求解方程组；
>> clear
>> format long
>> A=[1 1 1 1 1 1;1 1.1+0.000001 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+0.000001] >> B =[6.000000 7.715600 9.929900 12.756000 16.323800 20.781300]'
B =
6.000000000000000
7.715600000000000
9.929900000000000
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 r>k
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(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
>> x
x =
0.825832305593523 1.742109746543727 -0.259537880024995 2.064585114452380 0.551832152872303 1.075178560563062
题目二：多项式插值
在区间[5,5]-上对龙格函数2
1()1R x x =+做插值，分别用等距节点(节点步长1h =)和切比雪夫多项式11()T x 的零点做插值节点，画出原函数和两个插值函数的图像进行比较，并利用||||∞⋅对两种插值方法做误差分析。

原函数：
>> x=-5:0.1:5;
>> y=1./(1+x.^2);
>> plot(x,y)
2. 利用拉格朗日差分(等距节点)：
>> x=-5:1:5
x =
-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(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
>> p
p =
(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．切比雪夫零点插值：
11()
T x在[-1,1]上有11个不同零点：
(21)
cos
22
k
k
x
π
-
=(1,2,
k= (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.5406 -0.7557 -0.9096 -0.9898
>> x=5*s
x =
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)
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; %检错
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) 02,
0.55;r x ==
clear;
clc;
x=0.55;
r=2;
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) 01.5,0.55;r x ==
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) 01.5,
0.85;r x ==
clear;
clc;
x=0.85;
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.99999999 y=0.99999999 k=12>>
>> f=2*y^3-5*y^2+4*y-1
f =。