求解多元非线性方程组F(x)=0的Newton-Raphson方法及其MATLAB实现

牛顿迭代法可以推广到多元非线性方程组 \(\boldsymbol{F}(\boldsymbol{x})=\boldsymbol{0}\)的情况，称为牛顿-- 拉夫逊方法 (Newton-Raphson method). 当 \(\boldsymbol{F}(\boldsymbol{x})\) 关于 \(\boldsymbol{x}\) 的 Jacobi 矩阵 \(\boldsymbol{J}(\boldsymbol{x}) = (\cfrac{\partial \boldsymbol{F}}{\partial \boldsymbol{x}})\) 可逆时, 有

\[\boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}- \boldsymbol{J}^{-1}(\boldsymbol{x}^{(k)})\boldsymbol{F}(\boldsymbol{x}^{(k)}), \]

求解非线性方程组的Newton-Raphson方法：
1、 取初始点 \(\boldsymbol{x}^{(0)}\),最大迭代次数 \(N\) 和精度要求 \(\varepsilon\), 置 \(k=0\);
2、 求解线性方程组 \(\boldsymbol{J}(\boldsymbol{x}^{(k)})\boldsymbol{d} = -\boldsymbol{F}(\boldsymbol{x}^{(k)})\);
3、 若 $|\boldsymbol{d}|<\varepsilon $, 则停止计算；否则，置
\(\boldsymbol{x}^{(k+1)} = \boldsymbol{x}^{(k)}+\boldsymbol{d}^{(k)}\);
4、 若 \(k=N\), 则停止计算；否则，置 \(k = k+1\), 转(2).

function [x_star, it] = NewtonRapshon(fun,dfun,x0,ep,it_max)
if nargin <6 it_max = 100; end
if nargin<5 ep =1e-5; end
k = 1; err = 0; n = length(x0);
while k<it_max
    d = -feval(dfun,x0)\feval(fun,x0);
    x1 = x0 + d;
    err = norm(d,inf);
    if err<ep break; end
    x0 = x1; k = k+1;
end
x_star = x1; it = k;

用 Newton 法求非线性方程组$$\begin{cases}
x_1^2+x_2^2-5=0,\
(x_1+1)x_2-(3x_1+1)=0,
\end{cases}$$
取初始点\(\boldsymbol{x}^{(0)}=[1,1]^T\), 精度要求 \(\varepsilon = 10^{-3}\).
解：$$\boldsymbol{F}(\boldsymbol{x}) =\begin{bmatrix}
f_1(x_1,x_2)\
f_2(x_1,x_2)\
\end{bmatrix}= \begin{bmatrix}
x_1^2+x_2^2-5\
(x_1+1)x_2-(3x_1+1)\
\end{bmatrix}$$

\[\boldsymbol{J}(\boldsymbol{x}) = \begin{bmatrix} \dfrac{\partial f_1}{\partial x_1}&\dfrac{\partial f_1}{\partial x_2}\\ \dfrac{\partial f_2}{\partial x_2}&\dfrac{\partial f_2}{\partial x_2}\\ \end{bmatrix}= \begin{bmatrix} 2x_1&2x_2\\ x_2-3&x_1+1\\ \end{bmatrix}\]

所以有

\[\boldsymbol{x}^{(0)}=\begin{bmatrix} 1\\1\\ \end{bmatrix}, \boldsymbol{F}(\boldsymbol{x}^{(0)})=\begin{bmatrix} -3\\-2\\ \end{bmatrix},\boldsymbol{J}(\boldsymbol{x}^{(0)})=\begin{bmatrix} 2&2\\-2&2\\ \end{bmatrix},\]

解方程$$\begin{bmatrix}
2&2\-2&2\
\end{bmatrix}\begin{bmatrix}
d_1\d_2\
\end{bmatrix}=-\begin{bmatrix}
-3\-2\
\end{bmatrix}$$
得到 \(\boldsymbol{d}^{(0)}=[0.25,1.25]^T\), \(\boldsymbol{x}^{(1)} = \boldsymbol{x}^{(0)}+\boldsymbol{d}^{(0)}=[1.25,2.25]^T\),进入下一轮循环。
在 MATLAB 中，

>> [x_star, it] = NewtonRapshon(Fun,JFun,[1;1])

x_star =
    1.0000
    2.0000
it =
     5