【运筹优化】牛顿法详解 + Matlab代码实现
1 牛顿法简介
牛顿迭代法（Newton’s method）又称为牛顿-拉夫逊（拉弗森）方法（Newton-Raphson method），它是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法。
多数方程不存在求根公式，因此求精确根非常困难，甚至不可解，从而寻找方程的近似根就显得特别重要。方法使用函数 f ( x ) f(x) f(x) 的泰勒级数的前面几项来寻找方程 f ( x ) = 0 f(x)=0 f(x)=0 的根。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程 f ( x ) = 0 f(x)=0 f(x)=0 的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根，此时线性收敛，但是可通过一些方法变成超线 性收㪉。另外该方法广泛用于计算机编程中。
2 牛顿法原理
设
r
r
r 是
f
(
x
)
=
0
f(x)=0
f(x)=0 的根，选取
x
0
x_{0}
x0 作为
r
r
r 的初始近似值，过点
(
x
0
,
f
(
x
0
)
)
\left(x_{0}, f\left(x_{0}\right)\right)
(x0,f(x0)) 做曲线
y
=
f
(
x
)
y=f(x)
y=f(x) 的切线
L
L
L ，
L
:
y
=
f
(
x
0
)
+
f
′
(
x
0
)
(
x
−
x
0
)
L: y=f\left(x_{0}\right)+f^{\prime}\left(x_{0}\right)\left(x-x_{0}\right)
L:y=f(x0)+f′(x0)(x−x0) ，则
L
L
L 与
x
x
x 轴交点的横坐标
x
1
=
x
0
−
f
(
x
0
)
f
′
(
x
0
)
x_{1}=x_{0}-\frac{f\left(x_{0}\right)}{f^{\prime}\left(x_{0}\right)}
x1=x0−f′(x0)f(x0) ，称
x
1
x_{1}
x1 为
r
r
r 的一次近似值。过点
(
x
1
,
f
(
x
1
)
)
\left(x_{1}, f\left(x_{1}\right)\right)
(x1,f(x1)) 做曲线
y
=
f
(
x
)
y=f(x)
y=f(x) 的切线，并求该切线与
×
\times
× 轴交点的横坐标
x
2
=
x
1
−
f
(
x
1
)
f
′
(
x
1
)
x_{2}=x_{1}-\frac{f\left(x_{1}\right)}{f^{\prime}\left(x_{1}\right)}
x2=x1−f′(x1)f(x1) ，称
x
2
x_{2}
x2 为
r
\mathrm{r}
r 的二次近似值。重曷 以上过程，得
r
r
r 的近似值序列，其中，
x
n
+
1
=
x
n
−
f
(
x
n
)
f
′
(
x
n
)
x_{n+1}=x_{n}-\frac{f\left(x_{n}\right)}{f^{\prime}\left(x_{n}\right)}
xn+1=xn−f′(xn)f(xn) 称为
r
r
r 的
n
+
1
n+1
n+1 次近似值，上式称为牛顿迭代公式。
用牛顿迭代法解非线性方程，是把非线性方程 f ( x ) = 0 f(x)=0 f(x)=0 线性化的一种近似方法。把 f ( x ) f(x) f(x) 在点 x 0 x_{0} x0 的桌邻域内展开成泰勒 级数 f ( x ) = f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) + f ′ ′ ( x 0 ) ( x − x 0 ) 2 2 ! + ⋯ + f ( n ) ( x 0 ) ( x − x 0 ) n n ! + R n ( x ) f(x)=f\left(x_{0}\right)+f^{\prime}\left(x_{0}\right)\left(x-x_{0}\right)+\frac{f^{\prime \prime}\left(x_{0}\right)\left(x-x_{0}\right)^{2}}{2 !}+\cdots+\frac{f^{(n)}\left(x_{0}\right)\left(x-x_{0}\right)^{n}}{n !}+R_{n}(x) f(x)=f(x0)+f′(x0)(x−x0)+2!f′′(x0)(x−x0)2+⋯+n!f(n)(x0)(x−x0)n+Rn(x) ，取其线性部分 (即泰勒展开的前两项)，并令其等于 0 ，即 f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) = 0 f\left(x_{0}\right)+f^{\prime}\left(x_{0}\right)\left(x-x_{0}\right)=0 f(x0)+f′(x0)(x−x0)=0 ，以此作为非线性方程 f ( x ) = 0 f(x)=0 f(x)=0 的近似方程， 若 f ′ ( x 0 ) ≠ 0 f^{\prime}\left(x_{0}\right) \neq 0 f′(x0)=0 ，则其解为 x 1 = x 0 − f ( x 0 ) f ′ ( x 0 ) x_{1}=x_{0}-\frac{f\left(x_{0}\right)}{f^{\prime}\left(x_{0}\right)} x1=x0−f′(x0)f(x0) ，这样，得到牛顿迭代法的一个朱代关系式: x n + 1 = x n − f ( x n ) f ′ ( x n ) x_{n+1}=x_{n}-\frac{f\left(x_{n}\right)}{f^{\prime}\left(x_{n}\right)} xn+1=xn−f′(xn)f(xn) 。
已经证明，如果是连续的，并且待求的零点是孤立的，那么在零点周围存在一个区域，只要初始值位于这个邻近区域内，那 么牛顿法必定收敛。并且，如果不为 0 , 那么牛顿法将具有平方收敛的性能. 粗略的说，这意味着每造代一次，牛顿法结果的有效 数字将增加一倍。
3 牛顿法推导
4 Matlab代码实现
下面用Matlab代码求解上面的示例。
clear;clc;
% 定义原函数
syms xx yy
fy(xx,yy) = 0.5 * xx^2 + 2 * yy^2;
% 确定迭代次数
n = 10
% 确定初始点
x0 = 1
y0 = 1
% 求初始点函数值
fy(x0,y0)
% 求函数梯度
xf = -5:0.2:5;
yf = xf';
ff = 0.5 * xf.^2 + 2 * yf.^2;
surf(xf,yf,ff)
xlabel('x')
ylabel('y')
zlabel('z')
view([119.1 40.8])
[fx,fy] = gradient(ff,0.2);
% 提取点初始点处的梯度值
t = (xf == x0) & (yf == y0);
indt = find(t);
f_grad = [fx(indt) fy(indt)]
% 求海森矩阵
syms x y
f(x,y) = 0.5 * x^2 + 2 * y^2;
H = hessian(f,[x,y])
% 迭代
for i=1:n
% 判断是否可以跳出(如果梯度向量都接近0，就跳出)
b = 0;
for j = 1:length(f_grad)
if f_grad(j) > 0.000001
b = 1;
break
end
end
if b==0
break
end
% 确定下降方向
d = -inv(H)*(f_grad)';
dk = d(x0,y0);
% 确定步长，牛顿法步长为1
a = 1;
% 获取下一状态的点
newX = [x0,y0] + dk' .* a
x0 = newX(1);
y0 = newX(2);
% 更新梯度信息
t = (xf == x0) & (yf == y0);
indt = find(t);
f_grad = [fx(indt) fy(indt)];
end
5 低版本Matlab报错
最近有朋友向我反应代码运行会报错，具体报错内容如下：
他使用的matlab版本是2016a，推测可能是低版本不支持(151)的矩阵和(511)的矩阵直接做运算，如果大家有遇到这样的报错的话，可以试一下将原代码的16、17行删去，换成以下代码应该就可以了：
yf = xf;
s = size(xf,2);
ff = zeros(s,s);
for i = 1 : s
for j = 1 : s
obj = 0.5 * xf(i)^2 + 2*yf(j)^2;
ff(i,j) = obj;
end
end