当 Jacobi 矩阵不易求解时，可以采用拟 Newton 法。拟 Newton 法可以认为是割线法在高维空间中的推广。求解非线性方程组的拟 Newton 方法如下：

<1> 取初始点 $boldsymbol{x}^{(0)}$, 初始矩阵$boldsymbol{A}^{(0)}$, 最大迭代次数 $N$ 和精度要求 $varepsilon$, 置 $k=0$;

<2> 求解线性方程组 $boldsymbol{A}_kboldsymbol{s} = -boldsymbol{F}(boldsymbol{x}^{(k)})$;

<3> 若 $|boldsymbol{s}^{(k)}|<varepsilon $, 则停止计算；否则，置       $boldsymbol{x}^{(k+1)} = boldsymbol{x}^{(k)}+boldsymbol{s}^{(k)}$;

<4> 令 $boldsymbol{y}^{(k)} = boldsymbol{F}(boldsymbol{x}^{(k+1)})-boldsymbol{F}(boldsymbol{x}^{(k)})$,置$$ boldsymbol{A}_{k+1} = boldsymbol{A}_k + dfrac{(boldsymbol{y}_{(k)} -boldsymbol{A}_k boldsymbol{s}^{(k)})(boldsymbol{s}^{(k)})^T}{(boldsymbol{s}^{(k)})^Tboldsymbol{s}^{(k)}};$$

<5> 若 $k=N$, 则停止计算；否则，置$k = k+1$, 转(2).

MATLAB程序为

function [x_star, it] = quasinewton(fun,x0,ep,it_max)

if nargin <5 it_max = 100; end

if nargin<4 ep =1e-5; end

k = 1; err = 0; n = length(x0);A0 = eye(n);

while k<it_max

   s = -A0feval(fun,x0);

   x1 = x0 + s;

   err = norm(s,inf);

   if err<ep break; end

   y = feval(fun,x1)-feval(fun,x0);

   A1 = A0 + (y-A0*s)*s'/(s'*s);

   x0 = x1; A0 = A1;

   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}

$$ $\begin{cases} x_1^2+x_2^2-5=0,\\ (x_1+1)x_2-(3x_1+1)=0, \end{cases}$ $\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}$，写出MATLAB程序，并与用fsolve()得到的结果比较。

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

x_star =

1.0000

2.0000

it =

10

>> [x_star, err] = fsolve(Fun,[1;1])

Equation solved.

fsolve completed because the vector of function values is near zero

as measured by the default value of the function tolerance, and

the problem appears regular as measured by the gradient.

<stopping criteria details>

x_star =

   1.0000

   2.0000

err =

  1.0e-09 *

   0.8088

   0.3489