通过GSL2.7解非线性方程组的多维根

本文翻译自：Multidimensional Root-Finding
本章介绍多维求根函数（求解具有 n 个未知数的 n 个方程的非线性系统）。该库为各种迭代求解器和收敛测试提供了低级组件。用户可以将这些组合起来以实现所需的解决方案，并可以完全访问迭代的中间步骤。每一个类方法都使用相同的框架，因此您可以在运行时在求解器之间切换，而无需重新编译程序。求解器的每个实例都跟踪自己的状态，从而允许在多线程程序中使用求解器。求解器基于原始 Fortran 库 MINPACK。

头文件 gsl_multiroots.h 包含多维求根函数的原型和相关声明。

多维求根问题需要在 n 个变量 $x_i$ 中同时解 n 个方程 $f_i$

$f_i(x_1, x_2, ...,x_n)=0, i =1,2...,n$

一般来说，没有可用于 n 维系统的闭式方法，也无法知道是否存在任何解决方案。所有算法都使用牛顿迭代的变体从初始猜测开始，

$x = x'-J^{-1}f(x)$

其中 $x, f$ 是向量， $J$ 是雅格比矩阵， $J_{ij}=\frac{\partial{f_i}}{\partial{x_j}}$ 。可以使用其他策略来扩大收敛区域。其中包括要求降低范数 $∣ f ∣$ 在牛顿方法提出的每一步上，或在 $∣ f ∣$ 的负梯度方向上采取最陡下降步骤。

在一个框架中可以使用多种寻根算法。用户为算法提供高级驱动程序，库为每个步骤提供所需的单独函数。迭代过程有三个主要阶段。步骤是：

为算法 T 初始化求解器状态 s
使用迭代 T 更新 s
测试 s 的收敛性，并在必要时重复迭代

雅可比矩阵的评估可能会出现问题，要么是因为对导数进行编程是困难的，要么是因为矩阵的 $n^2$ 项的计算变得过于昂贵。由于这些原因，该库提供的算法根据导数是否有效分为两类。

具有解析雅可比矩阵的求解器的状态保存在 gsl_multiroot_fdfsolver 结构中。更新过程需要用户提供函数及其导数。
不使用解析雅可比矩阵的求解器的状态保存在 gsl_multiroot_fsolver 结构中。更新过程仅使用函数评估（而不是导数）。该算法通过近似方法估计矩阵 $J$ 或 $J^{-1}$ 。

初始化求解器

以下函数初始化一个多维求解器，有或没有导数。求解器本身仅取决于问题的维度和算法，并且可以重复用于不同的问题。

type gsl_multiroot_fsolver

这是一个无需导数的多维求根工作区

type gsl_multiroot_fdfsolver

这是一个使用导数进行多维求根的工作空间

gsl_multiroot_fsolver *gsl_multiroot_fsolver_alloc(const gsl_multiroot_fsolver_type *T, size_t n)

该函数返回一个指针，指向一个新分配的 T 类型求解器实例，用于 n 维系统。例如，以下代码创建一个混合求解器的实例，以求解一个 3 维方程组：

const gsl_multiroot_fsolver_type * T = gsl_multiroot_fsolver_hybrid;
gsl_multiroot_fsolver * s = gsl_multiroot_fsolver_alloc (T, 3);

如果没有足够的内存来创建求解器，则该函数返回一个空指针，并使用错误代码 GSL_ENOMEM 调用错误处理程序。

gsl_multiroot_fdfsolver *gsl_multiroot_fdfsolver_alloc(const gsl_multiroot_fdfsolver_type *T, size_t n)

此函数返回一个指针，指向一个新分配的 T 类型导数求解器实例，用于 n 维系统。例如，以下代码为二维方程组创建 Newton-Raphson求解器的实例：

const gsl_multiroot_fdfsolver_type * T = gsl_multiroot_fdfsolver_newton;
gsl_multiroot_fdfsolver * s = gsl_multiroot_fdfsolver_alloc (T, 2);

如果没有足够的内存来创建求解器，则该函数返回一个空指针，并使用错误代码 GSL_ENOMEM 调用错误处理程序。

int gsl_multiroot_fsolver_set(gsl_multiroot_fsolver *s, gsl_multiroot_function *f, const gsl_vector *x)
int gsl_multiroot_fdfsolver_set(gsl_multiroot_fdfsolver *s, gsl_multiroot_function_fdf *fdf, const gsl_vector *x)

这些函数设置或重置现有求解器 s 以使用函数 f 或函数和导数 fdf，以及初始猜测 x。请注意，初始位置是从 x 复制的，此参数不会被后续迭代修改。

void gsl_multiroot_fsolver_free(gsl_multiroot_fsolver *s)
void gsl_multiroot_fdfsolver_free(gsl_multiroot_fdfsolver *s)

这些函数释放与求解器相关的所有内存。

const char *gsl_multiroot_fsolver_name(const gsl_multiroot_fsolver *s)
const char *gsl_multiroot_fdfsolver_name(const gsl_multiroot_fdfsolver *s)

这些函数返回一个指向求解器名称的指针。例如：

printf ("s is a '%s' solver\n", gsl_multiroot_fdfsolver_name (s));

将会打印：s is a 'newton' solver

提供求解函数

你必须为根查找器提供 n 个变量的 n 个函数才能对其进行操作。为了允许通用参数，函数由以下数据类型定义：

type gsl_multiroot_function

这种数据类型定义了一个带有参数的通用函数系统

int (* f) (const gsl_vector * x, void * params, gsl_vector * f)

此函数应将向量结果 f(x,params) 存储在 f 中，用于参数 x 和参数 params，如果无法计算函数，则返回适当的错误代码。
size_t n

系统的维数，即向量 x 和 f 的分量数。
void * params

指向函数参数的指针。

这是一个使用Powell测试函数的例子

$f_1(x)=Ax_0x_1-1-1 \\ f_2(x)=exp(-x_0)+exp(-x_1)-(1+1/A)$

其中， $A=10^4$ 。以下代码定义了一个 gsl_multiroot_function 系统 F，您可以将其传递给求解器。

struct powell_params {
    
     double A; };

int
powell (gsl_vector * x, void * p, gsl_vector * f) {
    
    
   struct powell_params * params
     = (struct powell_params *)p;
   const double A = (params->A);
   const double x0 = gsl_vector_get(x,0);
   const double x1 = gsl_vector_get(x,1);

   gsl_vector_set (f, 0, A * x0 * x1 - 1);
   gsl_vector_set (f, 1, (exp(-x0) + exp(-x1)
                          - (1.0 + 1.0/A)));
   return GSL_SUCCESS
}

gsl_multiroot_function F;
struct powell_params params = {
    
     10000.0 };

F.f = &powell;
F.n = 2;
F.params = &params;

type gsl_multiroot_function_fdf

该数据类型定义了具有参数的一般函数系统和相应的导数Jacobian矩阵

int (* f) (const gsl_vector * x, void * params, gsl_vector * f)

此函数应将向量结果f(x，params)存储在f中的参数x和参数参数中，如果无法计算该函数，则返回适当的错误代码。
int (* df) (const gsl_vector * x, void * params, gsl_matrix * J)

此功能应存储 $n\times n$ 矩阵结果

$J_{ij} = \frac{\partial{f_i(x,params)}}{\partial{x_j}}$
J表示参数x和参数params，如果无法计算该函数，则返回适当的错误代码。

int (* fdf) (const gsl_vector * x, void * params, gsl_vector * f, gsl_matrix * J)

对于参数x和参数params，此函数应将f和J的值设置为上述。此函数提供了对F(x)和J(x)的单独函数的优化，即同时计算函数及其导数始终更快。
size_t n

系统的维数，即向量 x 和 f 的分量数。
void * params

指向函数参数的指针。

Powell在上述定义的测试功能的示例，可以扩展到使用以下代码来包括分析导数

int
powell_df (gsl_vector * x, void * p, gsl_matrix * J)
{
    
    
   struct powell_params * params
     = (struct powell_params *)p;
   const double A = (params->A);
   const double x0 = gsl_vector_get(x,0);
   const double x1 = gsl_vector_get(x,1);
   gsl_matrix_set (J, 0, 0, A * x1);
   gsl_matrix_set (J, 0, 1, A * x0);
   gsl_matrix_set (J, 1, 0, -exp(-x0));
   gsl_matrix_set (J, 1, 1, -exp(-x1));
   return GSL_SUCCESS
}

int
powell_fdf (gsl_vector * x, void * p,
            gsl_matrix * f, gsl_matrix * J) {
    
    
   struct powell_params * params
     = (struct powell_params *)p;
   const double A = (params->A);
   const double x0 = gsl_vector_get(x,0);
   const double x1 = gsl_vector_get(x,1);

   const double u0 = exp(-x0);
   const double u1 = exp(-x1);

   gsl_vector_set (f, 0, A * x0 * x1 - 1);
   gsl_vector_set (f, 1, u0 + u1 - (1 + 1/A));

   gsl_matrix_set (J, 0, 0, A * x1);
   gsl_matrix_set (J, 0, 1, A * x0);
   gsl_matrix_set (J, 1, 0, -u0);
   gsl_matrix_set (J, 1, 1, -u1);
   return GSL_SUCCESS
}

gsl_multiroot_function_fdf FDF;

FDF.f = &powell_f;
FDF.df = &powell_df;
FDF.fdf = &powell_fdf;
FDF.n = 2;
FDF.params = 0;

请注意，在计算Jacobian时，功能Powell_fDf能够从函数中重复使用现有项，从而节省了时间。

迭代

以下函数驱动每个算法的迭代。每个函数执行一次迭代以更新相应类型的任何求解器的状态。相同的函数适用于所有求解器，因此可以在运行时替换不同的方法，而无需修改代码。

int gsl_multiroot_fsolver_iterate(gsl_multiroot_fsolver *s)
int gsl_multiroot_fdfsolver_iterate(gsl_multiroot_fdfsolver *s)

这些函数执行求解器 s 的单次迭代。如果迭代遇到意外问题，则会返回错误代码。

GSL_EBADFUNC

迭代遇到了一个奇异点，函数或其导数计算为 Inf 或 NaN。
GSL_ENOPROG

迭代没有取得任何进展，从而阻止算法继续。

求解器始终保持对根 s->x 及其函数值 s->f 的当前最佳估计。可以使用以下辅助功能访问此信息。

gsl_vector *gsl_multiroot_fsolver_root(const gsl_multiroot_fsolver *s)
gsl_vector *gsl_multiroot_fdfsolver_root(const gsl_multiroot_fdfsolver *s)

这些函数返回求解器 s 的根的当前估计值，由 s->x 给出。

gsl_vector *gsl_multiroot_fsolver_f(const gsl_multiroot_fsolver *s)
gsl_vector *gsl_multiroot_fdfsolver_f(const gsl_multiroot_fdfsolver *s)

这些函数在求解器 s 的根的当前估计处返回函数值 f(x)，由 s->f 给出。

gsl_vector *gsl_multiroot_fsolver_dx(const gsl_multiroot_fsolver *s)
gsl_vector *gsl_multiroot_fdfsolver_dx(const gsl_multiroot_fdfsolver *s)

这些函数返回求解器 s 采取的最后一步 dx，由 s->dx 给出。

搜索停止参数

当以下条件之一为真时，根查找过程应停止：

已发现多维根在用户指定的精度范围内
已达到用户指定的最大迭代次数
发生了错误

这些条件的处理由用户控制。以下功能允许用户以多种标准方式测试当前结果的精度。

int gsl_multiroot_test_delta(const gsl_vector *dx, const gsl_vector *x, double epsabs, double epsrel)

此函数通过将最后一步 dx 与绝对误差 epsabs 和相对误差 epsrel 与当前位置 x 进行比较来测试序列的收敛性。如果满足以下条件，则测试返回 GSL_SUCCESS

$dx_i|<epsabs+epsrel|x_i|$

对于 x 的每个分量，否则返回 GSL_CONTINUE。

int gsl_multiroot_test_residual(const gsl_vector *f, double epsabs)

此函数根据绝对误差界限 epsabs 测试残差值 f。如果满足以下条件，则测试返回 GSL_SUCCESS，

$\sum_i|f_i|< epsabs$

否则返回 GSL_CONTINUE。此标准适用于根的精确位置 x 不重要的情况，只要可以找到残差足够小的值。

使用导数的算法

本节中描述的寻根算法利用了该函数及其导数。他们需要对根的位置进行初始猜测，但是没有绝对保证的收敛性，即该功能必须适合此技术，并且初始猜测必须足够接近根以使其工作。当满足条件时，收敛是二次的。

type gsl_multiroot_fdfsolver_type

gsl_multiroot_fdfsolver_type *gsl_multiroot_fdfsolver_hybridsj
gsl_multiroot_fdfsolver_type *gsl_multiroot_fdfsolver_hybridj
gsl_multiroot_fdfsolver_type *gsl_multiroot_fdfsolver_newton
gsl_multiroot_fdfsolver_type *gsl_multiroot_fdfsolver_gnewton

不使用导数的算法

本节中描述的算法不需要用户提供任何导数信息。所需的任何导数都通过有限差异近似。请注意，如果这些例程选择的有限差异步长不合适，则可以始终将明确的用户提供的数值导数与上一节中描述的算法一起使用。

type gsl_multiroot_fsolver_type

gsl_multiroot_fsolver_type *gsl_multiroot_fsolver_hybrids
gsl_multiroot_fsolver_type *gsl_multiroot_fsolver_hybrid
gsl_multiroot_fsolver_type *gsl_multiroot_fsolver_dnewton
gsl_multiroot_fsolver_type *gsl_multiroot_fsolver_broyden

例子

多维求解器的使用方式与一维求根算法类似。第一个示例演示了不需要导数的 HYBRIDS 比例混合算法。该程序求解 Rosenbrock 方程组

$f_1(x,y) = a(1-x) \\ f_2(x,y)=b(y-x^2)$

其中， $a = 1, b = 10$ 。该系统的解位于狭窄山谷 (x,y) = (1,1) 处。

程序的第一阶段是定义方程组

#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>

struct rparams
{
    
    
    double a;
    double b;
};

int
rosenbrock_f(const gsl_vector* x, void* params,
    gsl_vector* f)
{
    
    
    double a = ((struct rparams*)params)->a;
    double b = ((struct rparams*)params)->b;

    const double x0 = gsl_vector_get(x, 0);
    const double x1 = gsl_vector_get(x, 1);

    const double y0 = a * (1 - x0);
    const double y1 = b * (x1 - x0 * x0);

    gsl_vector_set(f, 0, y0);
    gsl_vector_set(f, 1, y1);

    return GSL_SUCCESS;
}

主程序首先创建函数对象 f，其参数为 (x,y) 和参数 (a,b)。使用 gsl_multiroot_fsolver_hybrids 方法初始化求解器 s 以使用此函数

int
main(void)
{
    
    
    const gsl_multiroot_fsolver_type* T;
    gsl_multiroot_fsolver* s;

    int status;
    size_t i, iter = 0;

    const size_t n = 2;
    struct rparams p = {
    
     1.0, 10.0 };
    gsl_multiroot_function f = {
    
     &rosenbrock_f, n, &p };

    double x_init[2] = {
    
     -10.0, -5.0 };
    gsl_vector* x = gsl_vector_alloc(n);

    gsl_vector_set(x, 0, x_init[0]);
    gsl_vector_set(x, 1, x_init[1]);

    T = gsl_multiroot_fsolver_hybrids;
    s = gsl_multiroot_fsolver_alloc(T, 2);
    gsl_multiroot_fsolver_set(s, &f, x);

    print_state(iter, s);

    do
    {
    
    
        iter++;
        status = gsl_multiroot_fsolver_iterate(s);

        print_state(iter, s);

        if (status)   /* check if solver is stuck */
            break;

        status =
            gsl_multiroot_test_residual(s->f, 1e-7);
    } while (status == GSL_CONTINUE && iter < 1000);

    printf("status = %s\n", gsl_strerror(status));

    gsl_multiroot_fsolver_free(s);
    gsl_vector_free(x);
    return 0;
}

请注意，检查每个求解器步骤的返回状态很重要，以防算法卡住。如果检测到错误条件，表明算法无法继续，则可以将错误报告给用户，选择新的起点或使用不同的算法。

解的中间状态由以下函数显示。求解器状态包含作为当前位置的向量 s->x，以及具有相应函数值的向量 s->f

void
print_state (size_t iter, gsl_multiroot_fsolver * s)
{
    
    
  printf ("iter = %3u x = % .3f % .3f "
          "f(x) = % .3e % .3e\n",
          iter,
          gsl_vector_get (s->x, 0),
          gsl_vector_get (s->x, 1),
          gsl_vector_get (s->f, 0),
          gsl_vector_get (s->f, 1));
}

以下是运行程序的结果。该算法从远离解的 (-10,-5) 处开始。由于解决方案隐藏在一个狭窄的山谷中，最早的步骤遵循函数的梯度下坡，试图减少残差的大值。一旦根大致定位，在第 8 次迭代中，牛顿行为接管并且收敛非常迅速：

iter =   0 x = -10.000 -5.000 f(x) =  1.100e+01 -1.050e+03
iter =   1 x = -10.000 -5.000 f(x) =  1.100e+01 -1.050e+03
iter =   2 x = -3.976  24.827 f(x) =  4.976e+00  9.020e+01
iter =   3 x = -3.976  24.827 f(x) =  4.976e+00  9.020e+01
iter =   4 x = -3.976  24.827 f(x) =  4.976e+00  9.020e+01
iter =   5 x = -1.274 -5.680 f(x) =  2.274e+00 -7.302e+01
iter =   6 x = -1.274 -5.680 f(x) =  2.274e+00 -7.302e+01
iter =   7 x =  0.249  0.298 f(x) =  7.511e-01  2.359e+00
iter =   8 x =  0.249  0.298 f(x) =  7.511e-01  2.359e+00
iter =   9 x =  1.000  0.878 f(x) = -2.653e-10 -1.218e+00
iter =  10 x =  1.000  0.989 f(x) = -2.353e-11 -1.080e-01
iter =  11 x =  1.000  1.000 f(x) =  0.000e+00  0.000e+00
status = success

请注意，该算法不会在每次迭代时更新位置。在尝试发现发散的步骤之后，一些迭代用于调整信任区域参数，或者在检测到收敛行为不佳时重新计算雅可比。

下一个示例程序添加导数信息，以加速求解。有两个导数函数rosenbrock_df 和rosenbrock_fdf。后者同时计算函数及其导数。这允许优化任何常用术语。为简单起见，我们在下面的代码中替换了对单独 f 和 df 函数的调用：

int
rosenbrock_df (const gsl_vector * x, void *params,
               gsl_matrix * J)
{
    
    
  const double a = ((struct rparams *) params)->a;
  const double b = ((struct rparams *) params)->b;

  const double x0 = gsl_vector_get (x, 0);

  const double df00 = -a;
  const double df01 = 0;
  const double df10 = -2 * b  * x0;
  const double df11 = b;

  gsl_matrix_set (J, 0, 0, df00);
  gsl_matrix_set (J, 0, 1, df01);
  gsl_matrix_set (J, 1, 0, df10);
  gsl_matrix_set (J, 1, 1, df11);

  return GSL_SUCCESS;
}

int
rosenbrock_fdf (const gsl_vector * x, void *params,
                gsl_vector * f, gsl_matrix * J)
{
    
    
  rosenbrock_f (x, params, f);
  rosenbrock_df (x, params, J);

  return GSL_SUCCESS;
}

主程序现在调用函数的相应 fdfsolver 版本

int
main (void)
{
    
    
  const gsl_multiroot_fdfsolver_type *T;
  gsl_multiroot_fdfsolver *s;

  int status;
  size_t i, iter = 0;

  const size_t n = 2;
  struct rparams p = {
    
    1.0, 10.0};
  gsl_multiroot_function_fdf f = {
    
    &rosenbrock_f,
                                  &rosenbrock_df,
                                  &rosenbrock_fdf,
                                  n, &p};

  double x_init[2] = {
    
    -10.0, -5.0};
  gsl_vector *x = gsl_vector_alloc (n);

  gsl_vector_set (x, 0, x_init[0]);
  gsl_vector_set (x, 1, x_init[1]);

  T = gsl_multiroot_fdfsolver_gnewton;
  s = gsl_multiroot_fdfsolver_alloc (T, n);
  gsl_multiroot_fdfsolver_set (s, &f, x);

  print_state (iter, s);

  do
    {
    
    
      iter++;

      status = gsl_multiroot_fdfsolver_iterate (s);

      print_state (iter, s);

      if (status)
        break;

      status = gsl_multiroot_test_residual (s->f, 1e-7);
    }
  while (status == GSL_CONTINUE && iter < 1000);

  printf ("status = %s\n", gsl_strerror (status));

  gsl_multiroot_fdfsolver_free (s);
  gsl_vector_free (x);
  return 0;
}

向 gsl_multiroot_fsolver_hybrids 求解器添加导数信息不会对其行为产生任何显着影响，因为它能够以足够的精度在数值上逼近雅可比行列式。为了说明不同导数求解器的行为，我们切换到 gsl_multiroot_fdfsolver_gnewton。这是一个传统的牛顿求解器，其约束条件是，如果整步会导致“上坡”，它会缩减步长。这是 gsl_multiroot_fdfsolver_gnewton 算法的输出

iter =   0 x = -10.000 -5.000 f(x) =  1.100e+01 -1.050e+03
iter =   1 x = -4.231 -65.317 f(x) =  5.231e+00 -8.321e+02
iter =   2 x =  1.000 -26.358 f(x) = -8.882e-16 -2.736e+02
iter =   3 x =  1.000  1.000  f(x) =  1.998e-15  4.441e-15
status = success

收敛速度要快得多，但需要大幅度偏移(-4.23，-65.3)。这可能会导致算法在实际应用中误入歧途。混合算法更可靠地遵循解决方案的下坡路径。

参考文献

Powell 的以下文章中描述了 Hybrid 方法的原始版本

M.J.D. Powell, “A Hybrid Method for Nonlinear Equations” (Chap 6, p 87–114) and “A Fortran Subroutine for Solving systems of Nonlinear Algebraic Equations” (Chap 7, p 115–161), in Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz, editor. Gordon and Breach, 1970.

以下论文也与本节中描述的算法有关

J.J. Moré, M.Y. Cosnard, “Numerical Solution of Nonlinear Equations”, ACM Transactions on Mathematical Software, Vol 5, No 1, (1979), p 64–85
C.G. Broyden, “A Class of Methods for Solving Nonlinear Simultaneous Equations”, Mathematics of Computation, Vol 19 (1965), p 577–593
J.J. Moré, B.S. Garbow, K.E. Hillstrom, “Testing Unconstrained Optimization Software”, ACM Transactions on Mathematical Software, Vol 7, No 1 (1981), p 17–41