拟Newton迭代法

本页目录

上面的Newton迭代法具有较好的收敛性，但是每一步都要计算\({F}^{'}({x}^{k})\)，很不方便，所以可以用比较简单的矩阵Ak代替Newton迭代法中的\({F}^{'}({x}^{k})\)，迭代公式为

\[{x}^{k+1}={x}^{k}-{A}_{k}^{-1}F({x}^{k})\]

下一步是确定Ak+1的值，由公式Ak+!=Ak+ΔAk可以求得Ak+1的值。根据求ΔAk的方法的不同，拟Newton迭代法可以被分为秩1的拟Newton迭代法（即Broyden迭代法）和逆Broyden秩1的迭代法。

1．秩1的拟Newton迭代法

秩1的拟Newton迭代法的迭代公式为

\[{x}^{k+1}={x}^{k}-{A}_{k}^{-1}F({x}^{k})\begin{matrix} {p}^{k}={x}^{k+1}-{x}^{k} & {q}^{k}=F({x}^{k+1})-F({x}^{k}) \end{matrix}\]
\[{A}_{k+1}={A}_{k}+\frac{({q}^{k}-{A}_{k}{p}^{k})({p}^{k}{)}^{T}}{{\left\| {p}^{k} \right\|}_{2}^{2}}\]

据此，可以编写秩1的拟Newton迭代法的M文件BroydenIterate.m，如下所示。

code.matlab 
function s=BroydenIterate(x,eps)
  % 用秩1的拟Newton法求非线性方程组
  % x为迭代初值，eps为允许误差值
  if nargin==1
      eps=1.0e-6;
  elseif nargin<1
      error
      return
  end
  % 第一次迭代
  x1=fx2(x);
  b1=inv(dfx2(x));
  p=-b1*x1';
  x3=x'+p;
  d=fx2(x3');
  q=(d-x1)';
  b1=b1+(q-b1*p)*p'/norm(p);
  while(norm(p)>=eps)%循环迭代
     % b1=b1+b3;
       x1=fx2(x3');
       p=-inv(b1)*(x1)';
       x3=x3+p;
       q=((fx2(x3'))-x1)';
       b1=b1+(q-b1*p)*p'/norm(p);
  end
  s=x3;
  return

【例8】用秩1的拟Newton迭代函数求解下列方程组

\[\left\{ \begin{aligned} \&{x}_{1}^{2}-{x}_{2}-1=0 \\ \&({x}_{1}-2{)}^{2}+({x}_{2}-0.5{)}^{2}-1=0 \end{aligned} \right.\]

下面用MATLAB实现，首先编写上述非线性方程组的M文件fx2.m：

code.matlab 
function y=fx2(x)
  y(1)=x(1)x(1)-x(2)-1;
  y(2)=(x(1)-2)(x(1)-2)+(x(2)-0.5)(x(2)-0.5)-1;
  y=[y(1) y(2)];

然后，再编写上述非线性方程组导数的M文件dfx2.m，如下所示。

code.matlab 
function y=dfx2(x)
  y(1)=2x(1);
  y(2)=-1;
  y(3)=2x(1)-4;
  y(4)=2x(2)-1;
  y=[y(1) y(2);y(3) y(4)];

最后，用编写的秩1的拟Newton迭代函数BroydenIterate进行求解。

code.matlab 
>> BroydenIterate([0 0])
ans =
    1.5463
    1.3912

2．逆Broyden秩1迭代法

逆Broyden秩1迭代法的迭代公式为

\[{x}^{k+1}={x}^{k}-{B}_{K}F({x}^{k})\begin{matrix} {p}^{k}={x}^{k+1}-{x}^{k} & \end{matrix}\]
\[{q}^{k}=F({x}^{k+1})-F({x}^{k})\]
\[{B}_{k+1}={B}_{k}+\frac{({p}^{k}-{B}_{k}{q}^{k})({p}^{k}{)}^{T}{B}_{k}}{({p}^{k}{)}^{T}{B}_{k}{q}^{k}}\]

据此，可以编写秩1的拟Newton法的M文件CBroydenIterate.m，如下所示。

code.matlab 
function s=CBroydenIterate(x,eps)
  % 用逆Broyden秩1迭代法求非线性方程组
  %x为迭代初值，eps为允许误差值
  if nargin==1
      eps=1.0e-6;
  elseif nargin<1
      error
      return
  end
  x1=fx3(x); %第一次迭代
  b1=inv(dfx3(x));
  p=-b1*x1';
  x3=x'+p;
  d=fx3(x3');
  q=(d-x1)';
  b1=b1+(p-b1*q)*p'*b1*inv(p'*b1*q);
  while(norm(p)>=eps) %循环迭代
     % b1=b1+b3;
       x1=fx3(x3');
       p=-b1*x1';
       x3=x3+p;
       q=((fx3(x3'))-x1)';
       b1=b1+(p-b1*q)*p'*b1*inv(p'*b1*q);
   end
   s=x3;
   return

【例9】用逆Broyden秩1迭代函数求解下列方程组

\[\left\{ \begin{matrix} {x}_{1}^{2}-10{x}_{1}+{x}_{2}^{2}+8=0 \\ {x}_{1}{x}_{2}^{2}+{x}_{1}-10{x}_{2}+8=0 \end{matrix} \right.\]

首先，编写上述非线性方程组的M文件fx3.m，如下所示。

code.matlab 
function y=fx3(x)
  y(1)=x(1)*x(1)-10*x(1)+x(2)*x(2)+8;
  y(2)=x(1)*x(2)*x(2)+x(1)-10*x(2)+8;
  y=[y(1) y(2)];

然后，再编写上述非线性方程组导数的M文件dfx.m，如下所示。

code.matlab 
function y=dfx3(x)
  y(1)=2*x(1)-10;
  y(2)=2*x(2);
  y(3)=x(2)*x(2)+1;
  y(4)=2*x(1)*x(2)-10;
  y=[y(1) y(2);y(3) y(4)];

最后，用编写的Broyden秩1迭代函数CBroydenIterate进行求解。

code.matlab 
>> CBroydenIterate([0 0])
ans =
    1.0000
    1.0000