MATLAB数值分析学习笔记：牛顿迭代法和割线法

牛顿迭代法 

在工程上所应用到的求根公式中，牛顿-拉弗森方法（The Newton-Raphson method）是使用的较多的一种方法：首先给定初始值，那么过（）作一条切线，其与轴的交点代表方程的数值解。 

使用斜率公式可以得到：

整理得到：

给定容差和最大迭代次数后，就可以迭代求解数值解了

• 代码实现

function [root,ea,iter] = Newt_raph(func,dfunc,xi,es,maxit)
%求根公式之牛顿-拉弗森方法
%%输入:
%func=将要解的函数方程
%dfunc=将要解的函数的导函数
%xi=（initial x）初始条件
%es=允许的容差
%maxit = 允许的最大迭代次数
%%输出：
%root=所求的根
%ea=相对误差（ea<=es时即可输出）
%iter=迭代次数
if nargin<3,error('最少需要输入三个参数——函数，导数与初始条件'), end
%当输入的参数少于3个时，即函数，导数与初始条件没有齐全时，报错！
if nargin<4|isempty(es),es=0.001;end
%当输入的参数少于4个时（没有输入容差），默认为0.001
if nargin<5|isempty(maxit),maxit=50;end
%当输入的参数少于5个时（没有输入最大迭代次数），默认为50次
iter = 0;
%实际迭代次数
while(1)
    x_old=xi;
    xi=xi - func(xi)/dfunc(xi);
    iter=iter+1;
    if xi ~= 0,ea = abs((xi-x_old)/xi); end
    if ea<=es | iter>=maxit,break,end
end
root = xi;
end

解释：

输入：func=将要解的函数方程，dfunc=将要解的函数的导函数，xi=（initial x）初始条件，es=允许的容差，maxit = 允许的最大迭代次数

输出：root=所求的根，ea=相对误差（ea<=es时即可输出），iter=迭代次数

• 问题求解

求函数的一个根：

f=@(x) 2^x-x^2;
%函数
df=@(x) log(2)/log(exp(1))*(2^x)-2*x;
%导函数
[root,ea,iter] = Newt_raph(f,df,1,1e-5,20)

解释：

由于MATLAB没有专门的导数函数，在这里考虑用解析法求导带入牛顿迭代法中

（或许可以使用diff函数）

• 结果：

>> Newt_raphtest
root =
     2
ea =
   1.7833e-08
iter =
     5

---

割线法（牛顿-拉弗森方法的改进）

正如上文所说，牛顿法的缺点是需要知道函数具体的导数，而有的时候函数的导数求解不便，可以用有限差分代替导数

注意和划界法的区别：割线法不需要估计值的函数值异号。

另一种计算导数的方法是给自变量增加扰动量Delte，即

得到如下的迭代方程

称之为改进割线法

• 代码实现

function [root,ea,iter] = secant(func,xi,delta,es,maxit)
%求根公式之割线法方法
%%输入:
%func=将要解的函数方程
%xi=（initial x）初始条件
%delta=扰动量
%es=允许的容差
%maxit = 允许的最大迭代次数
%%输出：
%root=所求的根
%ea=相对误差（ea<=es时即可输出）
%iter=迭代次数
if nargin<3,error('最少需要输入三个参数——函数，导数与初始条件'), end
%当输入的参数少于3个时，即函数，导数与初始条件没有齐全时，报错！
if nargin<4|isempty(es),es=0.001;end
%当输入的参数少于4个时（没有输入容差），默认为0.001
if nargin<5|isempty(maxit),maxit=50;end
%当输入的参数少于5个时（没有输入最大迭代次数），默认为50次
iter = 0;
%实际迭代次数
while(1)
    x_old=xi;
    xi=xi - delta*xi*func(xi)/(func(xi+delta*xi)-func(xi));
    iter=iter+1;
    if xi ~= 0,ea = abs((xi-x_old)/xi); end
    if ea<=es | iter>=maxit,break,end
end
root = xi;
end

 问题求解

设半径为r的圆环导体均匀带电为Q，轴线上的电荷q离圆环中心距离为多少时静电力为1.25N。

%%求解半径为r的均匀带电为Q的圆环导体对轴线上电荷q静电力为1.25N时，电荷离圆环中心距离x
%公式：F=q*Q*x/(4*pi*e0*(x^2+r^2)^1.5)
%其中e0=8.9e-12
%r=0.85,q=Q=2e-5
e0=8.9e-12;
f=@(x,q,Q,r) q*Q*x/(4*pi*e0*(x^2+r^2)^1.5)-1.25;
f(0,2e-5,2e-5,0.85)
f(0.5,2e-5,2e-5,0.85)
[root1,ea1,iter1] = secant(@(x) f(x,2e-5,2e-5,0.85),0.5,1e-6,1e-4,20)
[root2,ea2,iter2] = secant(@(x) f(x,2e-5,2e-5,0.85),0.25,1e-6,1e-4,20)
[root3,ea3,iter3] = secant(@(x) f(x,2e-5,2e-5,0.85),2,1e-6,1e-4,20)
[root4,ea4,iter4] = secant(@(x) f(x,2e-5,2e-5,0.85),1,1e-6,1e-4,20)

 解释：

只求解x>0的情况，注意初值分别取了0.5,0.25,2,1。

• 结果

>> secanttest
ans =
   -1.2500
ans =
    0.6146
root1 =
    0.2410
ea1 =
   5.4211e-06
iter1 =
     6
root2 =
    0.2410
ea2 =
   1.1046e-07
iter2 =
     3
root3 =
    1.2913
ea3 =
   8.0075e-06
iter3 =
     4
root4 =
    1.2913
ea4 =
   1.8957e-05
iter4 =
     3

发现有x>0时两个解0.2410和1.2913 

 声明：文章来源于笔者学习【美】Steven C. CHapra所著，林赐译 《工程于科学数值方法的MATLAB实现》（第4版）的笔记，如有谬误或想深入了解，请翻阅原书。