Newton迭代法算法报告.doc

实验题目4 Newton迭代法
摘要
为初始猜测，则由递推关系

产生逼近解的迭代序列，这个递推公式就是Newton法。当距较近时，很快收敛于。但当选择不当时，会导致发散。故我们事先规定迭代的最多次数。若超过这个次数，还不收敛，则停止迭代另选初值。
前言
利用牛顿迭代法求的根

程序设计流程
否
是
否
是
是
定义
输入
开 始

输出迭代
失败标志
输出
输出奇
异标志
结 束

否

问题1
（1） 程序运行如下：
r = NewtSolveOne('fun1_1',pi/4,1e-6,1e-4,10)
r = 0.7391
（2） 程序运行如下：
r = NewtSolveOne('fun1_2',0.6,1e-6,1e-4,10)
r = 0.5885
问题2
（1） 程序运行如下：
r = NewtSolveOne('fun2_1',0.5,1e-6,1e-4,10)
r = 0.5671
（2） 程序运行如下：
r = NewtSolveOne('fun2_2',0.5,1e-6,1e-4,20)
r = 0.5669
问题3
（1） 程序运行如下：
① p = LegendreIter(2)
p = 1.0000 0 -0.3333

p = LegendreIter(3)
p = 1.0000 0 -0.6000 0

p = LegendreIter(4)
p = 1.0000 0 -0.8571 0 0.0857

p = LegendreIter(5)
p = 1.0000 0 -1.1111 0 0.2381 0

② p = LegendreIter(6)
p = 1.0000 0 -1.3636 0 0.4545 0 -0.0216
r = roots(p)'
r = -0.932469514203150 -0.661209386466265 0.932469514203153
0.661209386466264 -0.238619186083197 0.238619186083197
用二分法求根为：
r = BinSolve('LegendreP6',-1,1,1e-6)
r = -0.932470204878826 -0.661212531887755 -0.238620057397959
0.238600127551020 0.661192602040816 0.932467713647959

（2） 程序运行如下：
① p = ChebyshevIter(2)
p = 1.0000 0 -0.5000

p = ChebyshevIter(3)
p = 1.0000 0 -0.7500 0

p = ChebyshevIter(4)
p = 1.0000 0 -1.0000 0 0.1250

p = ChebyshevIter(5)
p = 1.0000 0 -1.2500 0 0.3125 0

② p = ChebyshevIter(6)
p = 1.0000 0 -1.5000 0 0.5625 0 -0.0313
r = roots(p)'
r = -0.965925826289067 -0.707106781186548 0.965925826289068
0.707106781186547 -0.258819045102521 0.258819045102521

用二分法求根为：
r = BinSolve('ChebyshevT6',-1,1,1e-6)
r = -0.965929926658163 -0.707110969387755 -0.258828922193878
0.258818957270408 0.707105986926020 0.965924944196429

与下列代码结果基本一致，只是元素顺序稍有不同：
j = 0:5;
x = cos((2*j+1)*pi/2/(5+1))
x = 0.965925826289068 0.707106781186548 0.258819045102521
-0.258819045102521 -0.707106781186547 -0.965925826289068

（3） 程序运行如下：
① p = LaguerreIter(2)
p = 1 -4 2

p = LaguerreIter(3)
p = 1 -9 18 -6

p = LaguerreIter(4)
p = 1 -16 72 -96 24

p = LaguerreIter(5)
p =1.0000 -25.0000 200.0000 -600.0000 600.0000 -120.000

② p = LaguerreIter(5)
p =1.0000 -25.0000 200.0000 -600.0000 600.0000 -120.000
r = roots(p)'
r = 12.640800844275732 7.085810005858891 3.596425771040711
1.413403059106520 0.263560319718141

用二分法求根为：
r = BinSolve('LaguerreL5',0,13,1e-6)
r = 0.263560314567722 1.413403056105789 3.596425765631150
7.085810005360720 12.640800843813590

（4） 程序运行如下：
① p = HermiteIter(2)
p = 1.0000 0 -0.5000

p = HermiteIter(3)
p = 1.0000 0 -1.5000 0

p = HermiteIter(4)
p = 1.0000 0 -3.0000 0 0.7500

p = HermiteIter(5)
p = 1.0000 0 -5.0000 0 3.7500 0

② p = HermiteIter(6)
p = 1.0000 0 -7.5000 0 11.2500 0 -1.8750
r = roots(p)'
r = -2.350604973674487 2.350604973674488 -1.335849074013696
1.335849074013698 -0.436077411927617 0.436077411927616

用二分法求根为：
r = BinSolve('HermiteH6',-3,3,1e-6)
r = -2.350604981792216 -1.335849100229691 -0.436077818578603
0.436077351472816 1.335848983453244 2.350604952598104

所用到的函数
function r = NewtSolveOne(fun, x0, ftol, dftol, maxit)
% NewtSolveOne 用Newton法解方程f(x)=0在x0附近的一个根
%
% Synopsis: r = NewtSolveOne(fun, x0)
% r = NewtSolveOne(fun, x0, ftol, dftol)
%
% Input: fun = (string) 需要求根的函数及其导数
% x0 = 猜测根，Newton法迭代初始值
% ftol = (optional)误差，默认为5e-9
% dftol = (optional)导数容忍最小值，小于它表明Newton法失败，默认为5e-9
% maxit = (optional)迭代次数，默认为25
%
% Output: r = 在寻根区间内的根或奇点
if nargin < 3
ftol = 5e-9;
end
if nargin < 4
dftol = 5e-9;
end
if nargin < 5
maxit = 25;
end

x = x0; %设置初始迭代位置为x0
k = 0; %初始化迭代次数为0
while k <= maxit
k = k + 1;
[f,dfdx] = feval(fun,x); %fun返回f(x)和f'(x)的值

if abs(dfdx) < dftol %如果导数小于dftol，Newton法失败，返回空值
r = [];
warning('dfdx is too small!');
return;
end

dx = f/dfdx; %x(n+1) = x(n) - f( x(n) )/f'( x(n) )，这里设dx = f( x(n) )/f'( x(n) )
x = x - dx;

if abs(f) < ftol %如果误差小于ftol，返回当前x为根
r = x;
return;
end
end

r = []; %如果牛顿法未收敛，返回空值

function p = LegendreIter(n)
% LegendreIter 用递推的方法计算n次勒让德多项式的系数向量 Pn+2(x) = (2*i+3)/(i+2) * x*Pn+1(x) - (i+1)/(i+2) * Pn(x)
%
% Synopsis: p = LegendreIter(n)
%
% Input: n = 勒让德多项式的次数
%
% Output: p = n次勒让德多项式的系数向量
if round(n) ~= n | n < 0
error('n必须是一个非负整数');
end

if n == 0 %P0(x) = 1
p = 1;
return;
elseif n == 1 %P1(x) = x
p = [1 0];
return;
end

pBk = 1; %初始化三项递推公式后项为P0
pMid = [1 0]; %初始化三项递推公式中项为P1
for i = 0:n-2
pMidCal = zeros(1,i+3); %构造用于计算的x*Pn+1
pMidCal(1:i+2) = pMid;
pBkCal = zeros(1,i+3); %构造用于计算的Pn
pBkCal(3:i+3) = pBk;

pFwd = (2*i+3)/(i+2) * pMidCal - (i+1)/(i+2) * pBkCal; %勒让德多项式三项递推公式Pn+2(x) = (2*i+3)/(i+2) * x*Pn+1(x) - (i+1)/(i+2) * Pn(x)

pBk = pMid; %把中项变为后项进行下次迭代
pMid = pFwd; %把前项变为中项进行下次迭代
end

p = pFwd/pFwd(1); %把勒让德多项式最高次项系数归一化

function p = ChebyshevIter(n)
% ChebyshevIter 用递推的方法计算n次勒让德-切比雪夫多项式的系数向量 Tn+2(x) = 2*x*Tn+1(x) - Tn(x)
%
% Synopsis: p = ChebyshevIter(n)
%
% Input: n = 勒让德-切比雪夫多项式的次数
%
% Output: p = n次勒让德-切比雪夫多项式的系数向量
if round(n) ~= n | n < 0
error('n必须是一个非负整数');
end

if n == 0 %T0(x) = 1
p = 1;
return;
elseif n == 1 %T1(x) = x
p = [1 0];
return;
end

pBk = 1; %初始化三项递推公式后项为T0
pMid = [1 0]; %初始化三项递推公式中项为T1
for i = 0:n-2
pMidCal = zeros(1,i+3); %构造用于计算的x*Tn+1
pMidCal(1:i+2) = pMid;
pBkCal = zeros(1,i+3); %构造用于计算的Pn
pBkCal(3:i+3) = pBk;

pFwd = 2*pMidCal - pBkCal; %勒让德-切比雪夫多项式三项递推公式Tn+2(x) = 2*x*Tn+1(x) - Tn(x)

pBk = pMid; %把中项变为后项进行下次迭代
pMid = pFwd; %把前项变为中项进行下次迭代
end

p = pFwd/pFwd(1); %把勒让德-切比雪夫多项式最高次项系数归一化

function p = LaguerreIter(n)
% LaguerreIter 用递推的方法计算n次拉盖尔多项式的系数向量 Ln+2(x) = (2*n+3-x)*Ln+1(x) - (n+1)*Ln(x)
%
% Synopsis: p = LaguerreIter(n)
%
% Input: n = 拉盖尔多项式的次数
%
% Output: p = n次拉盖尔多项式的系数向量
if round(n) ~= n | n < 0
error('n必须是一个非负整数');
end

if n == 0 %L0(x) = 1
p = 1;
return;
elseif n == 1 %L1(x) = -x+1
p = [-1 1];
return;
end

pBk = 1; %初始化三项递推公式后项为L0
pMid = [-1 1]; %初始化三项递推公式中项为L1
for i = 0:n-2
pMidCal1 = zeros(1,i+3); %构造用于计算的x*Ln+1(x)
pMidCal1(1:i+2) = pMid;
pMidCal2 = zeros(1,i+3); %构造用于计算的Ln+1(x)
pMidCal2(2:i+3) = pMid;
pBkCal = zeros(1,i+3); %构造用于计算的Ln(x)
pBkCal(3:i+3) = pBk;

pFwd =( (2*i+3)*pMidCal2 - pMidCal1 - (i+1)*pBkCal )/ (i+2); %拉盖尔多项式三项递推公式Ln+2(x) = (2*n+3-x)*Ln+1(x) - (n+1)^2*Ln(x)

pBk = pMid; %把中项变为后项进行下次迭代
pMid = pFwd; %把前项变为中项进行下次迭代
end

p = pFwd/pFwd(1); %把拉盖尔多项式最高次项系数归一化

function p = HermiteIter(n)
% HermiteIter 用递推的方法计算n次埃尔米特多项式的系数向量 Hn+2(x) = 2*x*Hn+1(x) - 2*(n+1)*Hn(x)
%
% Synopsis: p = HermiteIter(n)
%
% Input: n = 埃尔米特多项式的次数
%
% Output: p = n次埃尔米特多项式的系数向量
if round(n) ~= n | n < 0
error('n必须是一个非负整数');
end

if n == 0 %H0(x) = 1
p = 1;
return;
elseif n == 1 %H1(x) = 2*x
p = [2 0];
return;
end

pBk = 1; %初始化三项递推公式后项为L0
pMid = [2 0]; %初始化三项递推公式中项为L1
for i = 0:n-2
pMidCal = zeros(1,i+3); %构造用于计算的x*Hn+1(x)
pMidCal(1:i+2) = pMid;
pBkCal = zeros(1,i+3); %构造用于计算的Hn(x)
pBkCal(3:i+3) = pBk;

pFwd =2*pMidCal - 2*(i+1)*pBkCal; %埃尔米特多项式三项递推公式Hn+2(x) = 2*x*Hn+1(x) - 2*(n+1)*Hn(x)

pBk = pMid; %把中项变为后项进行下次迭代
pMid = pFwd; %把前项变为中项进行下次迭代
end

p = pFwd/pFwd(1); %把拉盖尔多项式最高次项系数归一化

function r = BinSolve(fun, a, b, tol)
% BinSolve 用二分法解方程f(x)=0在区间[a,b]的根
%
% Synopsis: r = BinSolve(fun, a, b)
% r = BinSolve(fun, a, b, tol)
%
% Input: fun = (string) 需要求根的函数
% a,b = 寻根区间上下限
% tol = (optional)误差，默认为5e-9
%
% Output: r = 在寻根区间内的根
if nargin < 4
tol = 5e-9;
end
Xb = RootBracket(fun, a, b); %粗略寻找含根区间
[m,n] = size(Xb);
r = [];
nr = 1; %初始化找到的根的个数为1

maxit = 50; %最大二分迭代次数为50
for i = 1:m
a = Xb(i,1); %初始化第i个寻根区间下限
b = Xb(i,2); %初始化第i个寻根区间上限
err = 1; %初始化误差

k = 0;
while k < maxit
fa = feval(fun, a); %计算下限函数值
fb = feval(fun, b); %计算上限函数值
m = (a+b)/2;
fm = feval(fun, m);
err = abs(fm);
if sign(fm) == sign(fb) %若中点处与右端点函数值同号，右端点赋值为中点
b = m;
else %若中点处与左端点函数值同号或为0，左端点赋值为中点
a = m;
end

if err < tol %如果在a处函数值小于tol
r(nr) = a; %一般奇点不符合该条件，这样可以去除奇点
nr = nr + 1; %找到根的个数递增
k = maxit; %改变k值跳出循环
end

k = k + 1; %二分迭代次数递增
end
end

function X = powerX(x,a,b)
% powerX 对给定向量(x1, x2,..., xn)返回增幂矩阵(x1^a, x2^a,..., xn^a; x1^a+1, x2^a+1,..., xn^a+1; ...; x1^b, x2^b,..., xn^b;)
%
% Synopsis: X = powerX(x,a,b)
%
% Input: x = 需要返回增幂矩阵的向量
% a,b = 寻根区间上下限
%
% Output: X = 增幂矩阵(x1^a, x2^a,..., xn^a; x1^a+1, x2^a+1,..., xn^a+1; ...; x1^b, x2^b,..., xn^b;)
if round(a) ~= a | round(b) ~= b
error('a,b must be integers');
elseif a >= b
error('a must be smaller than b!');
end
x = x(:)';

row = b-a+1;
col = length(x);
X = zeros(row, col);
for i = b:-1:a
X(b-i+1,:) = x.^i;
End

function [f, dfdx] = fun1_1(x)
f = cos(x) - x;
dfdx = -sin(x) - 1;

function [f, dfdx] = fun1_2(x)
f = exp(-x) - sin(x);
dfdx = -exp(-x) - cos(x);

function [f, dfdx] = fun2_1(x)
f = x - exp(-x);
dfdx = 1 + exp(-x);

function [f, dfdx] = fun2_2(x)
f = x.^2 - 2*x*exp(-x) + exp(-2*x);
dfdx = 2*x - 2*exp(-x) + 2*x*exp(-x) - 2*exp(-2*x);

function y = LegendreP6(x)
p = LegendreIter(6);
X = powerX(x,0,6);
y = p*X;

function y = ChebyshevT6(x)
p = ChebyshevIter(6);
X = powerX(x,0,6);
y = p*X;

function y = LaguerreL5(x)
p = LaguerreIter(5);
X = powerX(x,0,5);
y = p*X;

function y = HermiteH6(x)
p = HermiteIter(6);
X = powerX(x,0,6);
y = p*X;

思考题
（1） 由于Newton法具有局部收敛性，所以在实际问题中，当实际问题本身能提供接近于根的初始近似值时，就可保证迭代序列收敛，但当初值难以确定时，迭代序列就不一定收敛。实际计算时应先用比较稳定的算法，如二分法，计算根的近似值，再将该近似值作为牛顿法的初值，以保证迭代序列的收敛性。
（2） 实验2中两个方程根其实相同，只是第二个方程为重根，通过比较迭代次数，第一个方程迭代了3次得出结果，第二个方程迭代了8次得出结果，且第二个方程的结果不如第一个准确，这是由于第二个方程在根处导数为0，在根的领域内导数很小使Newton法收敛速度变慢，精度变低。
（3） 我们来看下这几个多项式的图形：
Legendre P6 Chebyshev T6
Laguerre L5 Hermite H6
我们发现，这些多项式在比较小的区间内有多个根，这就致使其导数也会有多个根，因此如果用Newton法寻根的话初值非常不好估计，所以我们要用最稳定的二分法找它们的根。