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.3333p = LegendreIter(3)p = 1.0000 0 -0.6000 0p = LegendreIter(4)p = 1.0000 0 -0.8571 0 0.0857p = 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.0216r = 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.5000p = ChebyshevIter(3)p = 1.0000 0 -0.7500 0p = ChebyshevIter(4)p = 1.0000 0 -1.0000 0 0.1250p = 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.0313r = 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 2p = LaguerreIter(3)p = 1 -9 18 -6p = LaguerreIter(4)p = 1 -16 72 -96 24p = 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.000r = 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.5000p = HermiteIter(3)p = 1.0000 0 -1.5000 0p = HermiteIter(4)p = 1.0000 0 -3.0000 0 0.7500p = 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.8750r = 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 = 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;Endfunction 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 P6Chebyshev T6Laguerre L5Hermite H6我们发现，这些多项式在比较小的区间内有多个根，这就致使其导数也会有多个根，因此如果用Newton法寻根的话初值非常不好估计，所以我们要用最稳定的二分法找它们的根。