常州大学数值分析作业—第三章

第一章：9.设2
cos 1)(x
x
x f -=，给出计算函数值)012.0(f 的一个合适算法，并在字长m 给定的，十进制计算机上给出数值计算结果。

解：由 )24
21(242)2421(1)cos(1224242x x x x x x x -=-=+-
-≈- 得 )24
21(cos 1)(22x x x x f -≈-=
10. 字长为5的十进制计算机上计算
)015.0(f 和)015.0(g ，并与)015.0(f 的精确值
1.0075376410479比较，说明差异存在理由，其中x e x f x 1)(-=，24
621)(3
2x x x x g +++=。

clear
f=@(x)1/2-x^2/24; f(0.012)
ans =
0.5000
解：字长为5时的误差很大，这是因为设置的字长有限，就不可避免的使舍入误差不断积累。

把字长改为9时，误差已经大幅度减小。

这说明，加大字长可以显著减小误差。

11. 举例介绍数组矩阵常见运算。

解：举例如下
clear
f=@(x)digit(digit(exp(x)-1,5)/x,5);
g=@(x)digit(digit(1,5)+digit(x/2,5)+digit... (digit(x^2,5)/6,5)+digit(digit(x^3,5)/24,5),5); exc=1.0075376410479; f(0.015) g(0.015)
err1=f(0.015)-exc err2=g(0.015)-exc
ans =
1.0075 ans =
1.0075 err1 =
-3.7641e-05 err2 =
-3.7641e-05
clear
A=[1:4;5:8;9:12;13:16]
B=[1,1,1,1;2,2,2,2;3,3,3,3;4,4,4,4] A ’ A =
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 B =
1 1 1 1
2 2 2 2
3 3 3 3
4 4 4 4
ans =
1 5 9 13
2 6 10 14
3 7 11 15
4 8 12 16
A*B
ans = 30 30 30 30 70 70 70 70 110 110 110 110 150 150 150 150 A.*B
ans =
1 2 3 4
10 12 14 16
27 30 33 36
52 56 60 64
A^2
ans = 90 100 110 120 202 228 254 280 314 356 398 440 426 484 542 600 A.^2 ans =
1 4 9 16 25 36 49 64 81 100 121 144 169 196 225 256
%%编写m 文件使用digit 函数设置字长%% function y=digit(x,m) k=max(size(x)); y=x;
for i=1:k if x(i)<0 sign=-1; else
sign=1; end
x(i)=abs(x(i)); p=0;
if x(i)<0.1&x(i)>eps while x(i)<0.1 x(i)=x(i)*10; p=p-1; end end
if x(i)>=1
while x(i)>=1 x(i)=x(i)/10; p=p+1; end end
y(i)=round(x(i)*10^m)/10^m; y(i)=sign*y(i)*10^p; end return
f=@(x)digit(digit(exp(x)-1,9)/x,9);
g=@(x)digit(digit(1,9)+digit(x/2,9)+digit... (digit(x^2,9)/6,9)+digit(digit(x^3,9)/24,9),9); err1=f(0.015)-exc err2=g(0.015)-exc
err1 =
-1.0479e-09 err2 =
-1.0479e-09
12.对任意给定的实数a 、b 、c 、试编写Matlab 程序,求方程02
=++c bx ax 的根。

解：利用教材例11的方法： 当b>0时，a ac b b x 2421---=，b
ac b c x +--=4222。

当b<0时，b
ac b c
x --=4221，a ac
b b x 2422-+-=。

13.利用1,7
53arctan 7
53<+-+-=x x x x x x 及()
3/3arctan 6=π，给出一个计算
π的方法，根据此方法编写程序，给出π的至少有10位有效数字的近似值。

解：根据题中所给公式，容易得到：
()
1
2)3/3(16)3/3arctan(61
21
1
--≈=-=+∑i i n
i i π
A+B ans = 2 3 4 5 7 8 9 10 12 13 14 15 17 18 19 20
A-B ans =
0 1 2 3 3 4 5 6 6 7 8 9 9 10 11 12
clear
a=input('输入a='); b=input('输入b='); c=input('输入c='); d=b^2-4*a*c if b>=0
x1=(-b-d^0.5)/2*a x2=-2*c/(d^0.5+b) else b<0
x1=-2*c/(d^0.5-b) x2=(-b+d^0.5)/2*a end
clear
x=3^(-0.5);
y=atan(x) ; %精确值% s=0 ; %计算值% for k=1:2:100;
s=s+(-1)^((k+1)/2)*(x^k)/k;
err=y-s; m=abs(err) if m<=1e-11 break end end pi=6*s
输入a=2
输入b=2 输入c=1 d = -4 x1 =
-2.0000 - 2.0000i x2 =
-0.5000 + 0.5000i 输入a=1 输入b=2 输入c=1 d = 0 x1 = -1 x2 = -1 输入a=1 输入b=5 输入c=2 d = 17 x1 =
-4.5616 x2 =
-0.4384
14.分别利用下式给出计算ln2的近似方法，编写相应程序并比较算法运行情况。

11,32)
1()1ln(321
1
≤<-+++-=-=+∑∞
=+x n
x x x x n x x n
n n n 11),1253(21
2211ln 1253112<<-+-++++=-=-+-∞
=-∑x n x x x x n x x x n n n 解：
由运行结果可知，
方法二的绝对误差比方法一的误差要小得多。

这是因为方法一给出的计算公式含有相近数相减项，损失了有效数字。

而方法二给出的计算公式避免了相近数相减，具有较好的精度。

20.分别用Jacobi 迭代法、Gauss-seidel 迭代法解方程组
⎪⎩⎪
⎨⎧=++=++=-+.
122,1,122321
321321x x x x x x x x x 解：Jacobi 迭代法收敛⎥⎥⎥
⎦
⎤
⎢⎢⎢⎣⎡-=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡133321x x x ，Gauss-seidel 迭代法不收敛。

27.编写LU 分解法、改进平方根法、追赶法的Matlab 程序，并进行相关数值实验。

clear x=1/3; s=0;
for k=1:2:100; s=s+((x)^k)/k; end
y=2*s %计算值% g=log(2) %真实值% err=y-g %绝对误差% y =
0.6931 g =
0.6931 err =
-2.2204e-16
clear x=1; s=0;
for k=1:100;
s=s+((-1)^(k+1))*(x^k)/k; end
y=s %计算值% g=log(2) %真实值% err=y-g %绝对误差% y =
0.6882 g =
0.6931 err =
-0.0050
3.将矩阵⎥
⎥
⎥⎥⎦
⎤
⎢⎢⎢
⎢⎣⎡-=110011021110
02
01A 进行Doolittle 和Crout 分解 解：Doolittle 分解：结果如下，程序见后面。

⎥⎥
⎥⎥⎦⎤⎢
⎢⎢⎢⎣⎡-⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢
⎣⎡-==2.10001500
1110
02
0112.00001020010
00
01
LU A Crout 分解：结果如下，程序见后面。

⎥
⎥
⎥⎥⎦⎤⎢⎢⎢
⎢⎣⎡-⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢
⎣⎡-==100
02.0100
11100201
2.110
05020010
0001LU A 7.用改进平方根法解方程组⎥⎥⎥⎥
⎦
⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢
⎢⎣⎡----000142002511013101144321x x x x 解：结果如下，程序见后面。

⎥⎥⎥⎥
⎦
⎤⎢⎢⎢
⎢⎣⎡--=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡0256.00513.00769.02821.04321x x x x 8（2）.用追赶法求解方程组⎥⎥⎥⎥⎥
⎥⎦⎤⎢
⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⋅⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡--------00001210001210001210001210001254321x x x x x
解：结果如下，程序见后面。

⎥⎥⎥⎥
⎦⎤⎢⎢⎢
⎢⎣
⎡---=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡1429.02413.02857.03571.04321x x x x
function [x,iternum,flag]=jacobi(A,b,x0,delta,max1) %检验输入参数，初始化
if nargin<2,error('more augments are needed');end if nargin<3,x0=zeros(size(b));end if nargin<4,delta=1e-13;end if nargin<5,max1=100;end
if nargin>5,error('incorrect number of input');end n=length(b);x=0*b;flag=0;iternum=0;
%用Jacobi 迭代法解方程组 for k=1:max1
iternum=iternum+1; for i=1:n
if abs(A(i,i))<eps
error('A(i,i) equal to zero,divided by zero'); end
function [x,iternum,flag]=gseid(A,b,x0,delta,max1) %检验输入参数，初始化
if nargin<2,error('more augments are needed');end if nargin<3,x0=zeros(size(b));end if nargin<4,delta=1e-13;end if nargin<5,max1=100;end
if nargin>5,error('incorrect number of input');end n=length(b);flag=0;iternum=0;
%用Gauss-Seidel 迭代法解方程组 for k=1:max1
iternum=iternum+1; x=x0;
for i=1:n
if abs(A(i,i))<eps
error('A(i,i) equal to zero,divided by zero');
A=[1 2 -2;1 1 1;2 2 1] b=[1;1;1] A =
1 2 -2 1 1 1 2 2 1 b = 1 1 1 jacobi(A,b) The Jacobi method converges. ans = -3 3 1
gseid(A,b)
The Gauss-seidel method does not converge with 100 iterations
ans =
1.0e+31 *
-9.1905 9.2222 -0.0634
A=[1 2 -2;1 1 1;2 2 1] b=[1;1;1] A =
1 2 -2 1 1 1 2 2 1 b = 1 1 1 function [L,U]=Crout(A)
%检验输入参数，初始化 b=size(A);n=b(1); if b(1)~=b(2)
error('Matrix A should be a Square matrix.'); end
if n~=rank(A)
error('Matrix A should be FULL RANK.'); end
L=zeros(n,n);U=zeros(n,n); for i=1:n U(i,i)=1; end
%将矩阵A 进行Crout 分解 for k=1:n
A=[1 0 2 0;0 1 1 1;2 0 -1 1;0 0 1 1]
A =
1 0
2 0
0 1 1 1
2 0 -1 1
0 0 1 1 [L,U]=Crout(A)
L =
1.0000 0 0 0
0 1.0000 0 0
2.0000 0 -5.0000 0
0 0 1.0000 1.2000 U =
1.0000 0
2.0000 0
0 1.0000 1.0000 1.0000
0 0 1.0000 -0.2000
0 0 0 1.0000
function [L,U]=Crout(A)
%检验输入参数，初始化
b=size(A);n=b(1);
if b(1)~=b(2)
error('Matrix A should be a Square matrix.'); end
if n~=rank(A)
error('Matrix A should be FULL RANK.');
end
L=zeros(n,n);U=zeros(n,n);
for i=1:n
U(i,i)=1;
end
%将矩阵A进行Doolittle分解
for k=1:n
A=[1 0 2 0;0 1 1 1;2 0 -1 1;0 0 1 1]
A =
1 0
2 0
0 1 1 1
2 0 -1 1
0 0 1 1 [L,U]=Doolittle(A)
L =
1.0000 0 0 0
0 1.0000 0 0
2.0000 0 1.0000 0
0 0 -0.2000 1.0000 U =
1.0000 0
2.0000 0
0 1.0000 1.0000 1.0000
0 0 -5.0000 1.0000
0 0 0 1.2000
function [x]=improvecholesky(A,b,n)
%检验输入参数，初始化
L=zeros(n,n);D=diag(n,0);S=L*D;
for i=1:n
L(i,i)=1;
end
for i=1:n
for j=1:n
if (eig(A)<=0)|(A(i,j)~=A(j,i))
disp('Matrix A should be a Positive-definite matrix.');
break;
end
end
end
D(1,1)=A(1,1);
function [L,U,x]=pursue(a,b,c,f) n=length(b);u(1)=b(1); for i=2:n
if (u(i-1)~=0)
l(i-1)=a(i-1)/u(i-1); u(i)=b(i)-l(i-1)*c(i-1); else
break ; end end
L=eye(n)+diag(l,-1); U=diag(u)+diag(c,1);
x=zeros(n,1);y=x;y(1)=f(1); for i=2:n
y(i)=f(i)-l(i-1)*y(i-1);
A=[4 1 -1 0;1 3 -1 0;-1 -1 5 2;0 0 2 4]
b=[1;0;0;0]
A =
4 1 -1 0
1 3 -1 0
-1 -1 5 2
0 0 2 4
b =
1
n=4 [x]=improvecholesky(A,b,n) n = 4 x = 0.2821 -0.0769 0.0513 -0.0256
第三章：1.设节点x 0=0，x 1=π/8，x 2=π/4，x 3=3π/8，x 4=π/2，试适当选取上述节点，用Lagrange 插值法分别构造cos x 在区间[0，π/2]上的一次、二次、四次差值多项式P 1（x ），P 2（x ）和P 4（x ），并分别计算P 1（π/3），P 2（π/3）和P 4（π/3）。

function [A,Y]=lagrange(x,y,x0) %检验输入参数
if nargin <2 || nargin > 3
error('Incorrect Number of Inputs'); end
if length(x)~=length(y)
error('The length of x must be equal to it of y'); end
a=[1 -1 -1 -1];b=[2 2 2 2 2];c=[-1 -1 -1 -1]; f=[1 0 0 0 0 ];
[L,U,x]=pursue(a,b,c,f) L =
1.0000 0 0 0 0 0.5000 1.0000 0 0 0 0 -0.4000 1.0000 0 0 0 0 -0.6250 1.0000 0 0 0 0 -0.7273 1.0000 U =
2.0000 -1.0000 0 0 0 0 2.5000 -1.0000 0 0 0 0 1.6000 -1.0000 0 0 0 0 1.3750 -1.0000 0 0 0 0 1.2727 x =
0.3571 -0.2857 -0.2143 -0.1429
7.根据列表函数，选取适当的节点，用逐次线性插值法给出三次差值多项式在2.8处的值。

x 0.0 1.0 2.0 3.0 4.0 f （x ）
0.50
1.25
2.75
3.50
2.75
x = [0,pi/4,pi/2];
y= cos(x);
x0 = pi/3;
[A,Y] = lagrange(x,y,x0);
P2=vpa(poly2sym(A),3)
Y
P2 =
x^2 - 0.109*x - 0.336
Y =
0.5174 x = [pi/8,3*pi/8]; y = cos(x); x0 = pi/3; [A,Y] = lagrange(x,y,x0); P1 = vpa(poly2sym(A),3) Y P1 = 1.19*x - 0.689 Y = 0.4729 function [T]=aitken(x,y,x0,T0) if nargin == 3 T0=[]; end
n0=size(T0,1);m=max(size(x));n=n0+m;T=zeros(n,n+1); T(1:n0,1:n0+1)=T0;T(n0+1:n,1)=x;T(n0+1:n,2)=y; if n0==0 i0=2; x = [0,pi/8,pi/4,3*pi/8,pi/2]; y= cos(x); x0 = pi/3; [A,Y]=lagrange(x,y,x0); P4=vpa(poly2sym(A),3) Y P4 = x^4 + 0.00282*x^3 - 0.514*x^2 + 0.0232*x + 0.0287 Y = 0.5001
x=[0 1];
y=[0.5 1.25];
x0=2.8;
T0=aitken(x,y,x0);
T=T0;
x=[3.0,4.0]';y=[3.5,2.75]';x0=2.8;T=aitken(x,y,x0,T0);
n=max(size(x))+size(T0,1);
for i=1:n
for j=1:i+1
fprintf('%10.4f',T(i,j));
end
fprintf('\n');
end
return
0.0000 0.5000
1.0000 1.2500
2.6000
3.0000 3.5000 3.3000 3.2300
4.0000 2.7500 2.0750 2.2850 3.4190
16.选取适当的函数y=f（x）和插值节点，编写Matlab程序，分别给出利用Lagrange插值
方法、Newton插值方法确定的插值多项式，并将函数y=f（x）的插值多项式及插值余项的图形画在同一坐标系中，观测节点变化对插值余项的影响。

function [C,D,Y]=newpoly(x0,y0,x)
if nargin < 2 | nargin> 3
error('Incorrect Number of Inputs');
end
if length(x0)~=length(y0)
error('The length of x0 must be equal to it of y0');
end
n=length(x0);
fl(x)=0.5*x^4 - 0.312*x^3 + 1.47*x^2 - 0.438*x + 0.0312 fn(x)=0.5*x^4 - 0.312*x^3 + 1.47*x^2 - 0.438*x + 0.0312 第四章：6.已知列表函数，用最小二乘法求形如
bx axe y 的拟合函数。

x 1.0 2.0 3.0 4.0 5.0 f （x ）
1.222
2.984
5.466
8.902
13.592
2.
>> a = [1 2 3 4 5]; >> x = log(a);
>> y = [2 2.693 3.099 3.386 3.609]; >> n = 1;
>> [C]=lspoly(x,y,n);
x = [0 1 2 3 4 ];
y = [0.5,1.25,2.75,3.5,2.75];
[A,Y]=lagrange(x,y,x0)
x0 = [0 1 2 3 4 ];
y0 = [0.5,1.25,2.75,3.5,2.75];
[C,D,X]=newpoly(x0,y0,x)
plot(x,Y ,'b-',x0,X,'r:')
A =
0.5000 -0.3125 1.4687 -0.4375 0.0313
Y =
0.5000 1.2500 2.7500 3.5000 2.7500
C =
0.0313 -0.4375 1.4688 -0.3125 0.5000
D =
0.5000 0 0 0 0
1.2500 0.7500 0 0 0
2.7500 1.5000 0.3750 0 0
3.5000 0.7500 -0.3750 -0.2500 0
2.7500 -0.7500 -0.7500 -0.1250 0.0313
X =
0.5000 1.2500 2.7500 3.5000 2.7500
>> y = vpa(poly2sym(C),3) 结果如下：
y = 1.0*x + 2.0 表示成题目中拟合函数格式即为 ： y = 2 + ln(x) 4. 解： >> x = a; >> y = a./b; >> n = 1;
>> [C]=lspoly(x,y,n);
>> y = vpa(poly2sym(C),3)
有结果为：y = 2.0*x + 1.0 写成题目中的拟合函数形式则为 ： y = x / (2.0*x + 1.0) 5
对题中函数进行变形：
原式 → y/x = a* exp(b*x) → ln(y/x) = ln(a) + b*exp(x) 化为线性形式 计算：
>> a = [1 2 3 4 5];
>> b = [1.222 2.984 5.466 8.902 13.592]; >> x = exp(a);
>> y = log(b)-log(a); >> n = 1;
>> [C]=lspoly(x,y,n);
>> y = vpa(poly2sym(C),3) 结果如下：
y = 0.00464*x + 0.384
写成题中拟合函数的形式即为： y = 1.4679*x*exp(0.00464*x)
7.已知人体表面积S 和身高h ，体重w 有近似关系210a a w h a S
=，试根据身高、体重及
相应的人体表面积的一组观测值（hi ，wi ，Si ）（i=0,1,2，，n ）来估计参数a0a1a2的大小。

8.学习matlab 内部函数lsqcurvefit 的用法，并设计数值实验使用函数lsqcurvefit 。

第五章：4.分别用梯形公式、simpson 公式计算积分⎰
+dx x )12(10
，并与精确值比较。

function [s,t]=tiandsim(f,a,b) fa= feval(f,a); fb= feval(f,b);
f2= feval(f,(a+b)/2); t=(1/2)*(b-a)*(fa+fb);
s=(1/6)*(b-a)*(fa+fb+4*f2); return a=0;b=1; f=@(x) 2*x+1; %计算精确解 s0=int(f,x,a,b) %计算数值解
[s,t]=tiandsim(f,a,b) %计算误差 err1=s-s0 err2=t-s0
输出结果为： s0 =2 s =2 t =2
err1 =0 err2 =0
7.取n=2,3,4，用复合梯形公式、复合simpson 公式计算积分
dx xe x 2
10
，并与精确值比较。

13.根据列表函数，试给出至少三种不同的数值方法来计算函数y=lnx 在x=2.0处的导数，
并将各种方法的计算结果与精确值相比较。

x 1.8 1.9 2.0 2.1 2.2 ln （x ）
0.5878
0.6419
0.6931
0.7419
0.7885
function s=trapezoidal(f,a,b,n) if nargin == 3 n=20; end
if nargin<3||nargin>4||... a>b||n<1||n~=round(n)
error('Incorrect input arguments'); end
h=(b-a)/n;s=0; for k=1:(n-1)
x=a+h*k;s=s+feval(f,x); end
s=h*(feval(f,a)+feval(f,b))/2+h*s; return
function s=simpson(f,a,b,n)
if nargin == 3
n=20;
end
if nargin<3||nargin>4||...
a>b||n<1||n~=round(n)
error('Incorrect input arguments');
end
h=(b-a)/(2*n);s1=0;s2=0;
for k=1:n
x=a+h*(2*k-1);s1=s1+feval(f,x);
end
for k=1:(n-1)
x=a+h*2*k;s2=s2+feval(f,x);
end
s=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;
return format long
a=0;b=1; f=@(x) x*exp(x^2); %计算精确解 s0=int(f,x,a,b) %计算数值解 for n=2:1:4;
s1=trapezoidal(f,a,b,n) s2=simpson(f,a,b,n) err1=vpa(s1-s0,6) err2=vpa(s2-s0,4) end
输出结果为：
s0 = 0.859140914229523
s1 = 1.000576811286697 s2 = 0.860997139578795 err1 = 0.141436 err2 = 0.001856 s1 = 0.923798756293777 s2 = 0.859533825596209 err1 =0.0646578 err2 =0.0003929 s1 =0.895892057505771 s2 =0.859268455239111 err1 = 0.0367511
err2 =0.0001275 function [Dc,err]=dfDc(f,x0,h) d0=1/x0; Dc=(f(x0+h)-f(x0-h))/(2*h); err=Dc-d0; return
function [Sc,err]=dfSc(f,x0,h) f=@(x) log(x); x0=2;h=0.1;
[Dc,err]=dfDc(f,x0,h) [Sc,err]=dfSc(f,x0,h) [Cc,err]=dfCc(f,x0,h)
Dc =
第六章：3.编写matlab 程序，分别用二分法和试位法求方程0222
3
=--+x x x 的根，并给出各自达到精度要求所需计算函数值f(x)的次数，这里设2105.0-⨯=ε。

function [c,err,yc,k]=bisect(f,a,b,epsilon) yb=f(b);ya=f(a);max1=1+round((log(b-a)-log(epsilon))/log(2)); flag=1;k=0; while flag==1 k=1:max1; c=(a+b)/2; yc=f(c); if yc==0 function [c,err,yc,k]=fapo(f,a,b,epsilon,n) ya=f(a);yb=f(b);
flag=1;k=0; while flag==1 k=k+1; c=b-yb*(b-a)/(yb-ya); a=-2.5;b=-1.5; n=20;epsilon=0.5e-2; f=@(x)(x^3+2*x^2-x-2); [c,err,yc,k]=bisect(f,a,b,epsilon) function [Cc,err]=dfCc(f,x0,h) d0=1/x0; Cc=16/15*dfSc(f,x0,h/2)... -1/15*dfSc(f,x0,h); err=Cc-d0; return
4.用迭代方法求解方程05ln 22
=--x x 在给定x 0=3附近的根，要求误差不超过10-4。

10.分别用（1）Newton 法，取x0=-2，（2）割线法，取x0=-2，x1=-2.1，求方程0433
=+-x x 在x0=-2附近的根，并比较各算法的数值表现。

function [k,x,err,X]=fixpt(f,x0,tol,max1) X(1)=x0; for k=2:max1 X(k)=feval(f,X(k-1)); err=abs(X(k)-X(k-1)); relerr=err/(abs(X(k))+eps); x=X(k); if (err<tol)|(relerr<tol) break ; end end if k==max1 disp('maximum number exceeded') end X=X'; f=@(x) (2*log(x)+5)^(1/2); [k,x,err,X]=fixpt(f, 3,1e-4,20) 输出结果为： k = 6 x = 2.633810444323049 err = 1.240357194669528e-04 X = 3.000000000000000 2.682764353672573 2.640775544482588 2.634795111907815 2.633934480042516 2.633810444323049
function [x,err,k,y]=newton(f,df,x0,epsilon,max1) flag=0;
for k=1:max1
x=x0-f(x0)/df(x0); err=abs(x-x0); x0=x;y=f(x0);
if (err<epsilon)||(abs(y)<epsilon) flag=1;
function [x,err,k,y]=secant(f,x0,x1,epsilon,max1)
flag=0;
for k=1:max1
x=x0-f(x0)*(x1-x0)/(f(x1)-f(x0));
err=abs(x-x1);
x0=x1; x1=x; y=f(x0);
if(err<epsilon)||(abs(y)<epsilon)
flag=1;
break;
end
end
if flag ~= 1
disp(['secant Method failed to get the zero within ... 'num2str(max1) 'iterations'])
end
return
x0=2; x1=-2.1;
epsilon=1e-6; max1=100;
f=@(x) (x^3-3*x+4);
df=@(x) (3*x^2-3);
%用牛顿法
[x,err,k,y]=newton(f, df,x0,epsilon,max1) %用割线法
[x,err,k,y]=secant(f, x0,x1,epsilon,max1) 输出结果为：
x =
-2.1958
err =
3.1965e-04
k =
42
y =
-6.7313e-07
输出结果为：
x =
-2.1958
err =
1.4741e-10
k =
7
y =
1.6901e-09。