类型

牛顿法可用于找到函数 f 的一个根（零）x 的越来越好的近似值

x:f(x)=0\,.

如果想找到另一个根，必须使用其他初始点再次应用该方法。

如何找到所有根？ ^[3]^[4]

牛顿法可应用于

实值函数 ^[5]
复值函数
方程组

非线性方程组

也可以使用牛顿法求解^[6]

2 个（非线性）方程
包含 2 个复变量： $x_{1},x_{2}$ ^[7]

 ${\begin{cases}F_{1}(x_{1},x_{2})=0\\F_{2}(x_{1},x_{2})=0\end{cases}}$

可以使用向量表示更简洁地表达

 $\mathbf {F} (\mathbf {x} )=0$

其中 $\mathbf {x}$ 是一个包含 2 个复变量的向量

 $\mathbf {x} ={\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}$

而 $\mathbf {F}$ 是一个函数向量（或是一个向量变量 x 的向量值函数）

  $\mathbf {F} ={\begin{bmatrix}F_{1}\\F_{2}\end{bmatrix}}$

可以通过迭代求解：^[8]

 $\mathbf {x} ^{(k+1)}=\mathbf {x} ^{(k)}+\mathbf {s} ,\ k=0,1,2,\ldots .$

其中 s 是一个增量（现在是一个包含 2 个分量的向量）

  $\mathbf {s} =\mathbf {x} ^{(k+1)}-\mathbf {x} ^{(k)}$

增量可以通过以下公式计算：

 $\mathbf {s} =J^{-1}(\mathbf {x} ^{(k)})\mathbf {F} (\mathbf {x} ^{(k)})$

这里 $J$ 是 $\mathbf {F}$ 在 $\mathbf {x} ^{(k)}$ 处的雅可比矩阵^[9]

  $J_{\mathbf {F} }(\mathbf {x} ^{(k)})=(F'(\mathbf {x} ))_{ij}={\dfrac {\partial F_{i}}{\partial x_{j}}}(\mathbf {x} )={\begin{bmatrix}{\dfrac {\partial F_{1}}{\partial x_{1}}}&{\dfrac {\partial F_{1}}{\partial x_{2}}}\\[1em]{\dfrac {\partial F_{2}}{\partial x_{1}}}&{\dfrac {\partial F_{2}}{\partial x_{2}}}\end{bmatrix}}$

而 $J^{-1}$ 是雅可比矩阵 $J$ 的逆矩阵。

在以下情况下

少量方程使用逆矩阵。
大量方程式：与其实际计算该矩阵的逆，不如通过求解线性方程组 $J_{F}(x_{n})(x_{n+1}-x_{n})=-F(x_{n})\,\!$ 来节省时间，从而求得未知数 x_n+1 − x_n。

伪代码算法：^[10]

从初始近似值（猜测） $\mathbf {x} ^{k};k=0$ 开始
重复直到收敛（迭代）
- 求解： $\mathbf {s} =J^{-1}(\mathbf {x} ^{(k)})\mathbf {F} (\mathbf {x} ^{(k)})$
- 计算： $\mathbf {x} ^{(k+1)}=\mathbf {x} ^{(k)}+\mathbf {s}$
- $k=k+1$

停止条件

绝对误差： $|s|<1.e-12$ ^[11]

需要什么？

函数 f （可微函数）
它的导数 f'
起点 x0 （应该在根的吸引域内）
停止规则（停止迭代的标准）^[12]

停止规则

$|x_{n+1}-x_{n}|<\epsilon \;$

陷阱或失败

为什么方法无法找到根？^[13]

方法有循环，永远不会收敛（
- 拐点^[14]
- 局部最小值
方法发散而不是收敛
- 局部最小值
多重根（根的重数> 1，例如：f(x) = f'(x) = 0）。缓慢逼近，导数趋于零，除法步骤存在问题（不要除以零）。使用修正牛顿法
收敛速度慢^[15]
- 初始近似值（远离根的初始猜测）不好
- 导数在根附近消失的函数，或二阶导数在根附近无界的函数

函数有拐点（方法有循环，永远不会收敛）
局部最小值

"对于大多数问题，在双精度下应用牛顿法有 2-3 步全局定向，然后 4-6 步二次收敛，直到双精度浮点数格式的精度被超过。因此，可以更大胆一些，如果在 10 步后迭代没有收敛，那么初始点就很糟糕，无论它是导致周期性轨道还是发散到无穷大。最有可能的是，附近的初始点将表现出类似的行为，因此在开始下一次迭代运行时对初始点进行非平凡的更改。"LutzL^[16]

伪代码

以下是使用牛顿法来帮助找到函数 f 的根的示例，该函数的导数为 fprime。

初始猜测将为 $x_{0}=1$ ，函数将为 $f(x)=x^{2}-2$ ，因此 $f'(x)=2x$ 。

牛顿法的每次新迭代都将用 x1 表示。在计算过程中，我们将检查分母 (yprime) 是否变得太小（小于 epsilon），如果 $f'(x_{n})\approx 0$ ，就会出现这种情况，因为否则会导致大量的误差。

%These choices depend on the problem being solved
x0 = 1                      %The initial value
f = @(x) x^2 - 2            %The function whose root we are trying to find
fprime = @(x) 2*x           %The derivative of f(x)
tolerance = 10^(-7)         %7 digit accuracy is desired
epsilon = 10^(-14)          %Don't want to divide by a number smaller than this

maxIterations = 20          %Don't allow the iterations to continue indefinitely
haveWeFoundSolution = false %Were we able to find the solution to the desired tolerance? not yet.

for i = 1 : maxIterations

    y = f(x0)
    yprime = fprime(x0)

    if(abs(yprime) < epsilon)                         %Don't want to divide by too small of a number
        fprintf('WARNING: denominator is too small\n')
        break;                                        %Leave the loop
    end

    x1 = x0 - y/yprime                                %Do Newton's computation
    
    if(abs(x1 - x0)/abs(x1) < tolerance)              %If the result is within the desired tolerance
        haveWeFoundSolution = true
        break;                                        %Done, so leave the loop
    end

    x0 = x1                                           %Update x0 to start the process again                  
    
end

if (haveWeFoundSolution) % We found a solution within tolerance and maximum number of iterations
    fprintf('The root is: %f\n', x1);
else %If we weren't able to find a solution to within the desired tolerance
    fprintf('Warning: Not able to find solution to within the desired tolerance of %f\n', tolerance);
    fprintf('The last computed approximate root was %f\n', x1)
end

代码

c ^[17]
matlab : newtonm 函数 ^[18]
Maxima CAS
scilab ^[19]
Mathematica^[20]
用于解释数学视频的 Python 动画引擎

误差分析

牛顿法在迭代复二次多项式中的应用

递推关系

递推关系是一个递归定义点序列（例如轨道）的方程（映射）。要计算序列，需要

起始点
映射（方程 = 函数 = 递推关系）

定义

定义迭代函数 $f^{n}$ 为

${\begin{aligned}f^{0}(z,c)&=z\\f^{1}(z,c)&=z^{2}+c\\f^{m+1}(z,c)&=f(f^{m}(z,c),c)\end{aligned}}$

请注意另一个函数

  $F(z)=f(z)-z$

用于

查找周期点

为迭代函数及其导数选择**任意名称**。这些导数可以通过递推关系计算。

描述	数学符号	任意名称	初始条件	递推步骤	常用名称
迭代函数	$f^{m}$	$A_{m}$	$A_{0}=z$	$A_{m+1}=A_{m}^{2}+c$	z 或 f
对 z 的一阶导数	${\dfrac {\partial }{\partial z}}f^{m}$	$B_{m}$	$B_{0}=1$	$B_{m+1}=2A_{m}B_{m}$	d（来自导数）或 p（来自素数）的 f'
对 z 的二阶导数	${\dfrac {\partial }{\partial z}}{\dfrac {\partial }{\partial z}}f^{m}$	$C_{m}$	$C_{0}=0$	$C_{m+1}=2(B_{m}^{2}+A_{m}C_{m})$
对c的一阶导数	${\dfrac {\partial }{\partial c}}f^{m}$	$D_{m}$	$D_{0}=0$	$D_{m+1}=2A_{m}D_{m}+1$
对c的二阶导数	${\dfrac {\partial }{\partial c}}{\dfrac {\partial }{\partial z}}f^{m}$	$E_{m}$	$E_{0}=0$	$E_{m+1}=2(A_{m}E_{m}+B_{m}D_{m})$

递推关系的推导

第一个基本规则

迭代函数 $f^{m+1}(z,c)=f(f^{m}(z,c),c)$ 是一个函数复合
应用链式法则对于复合函数，因此如果 $f(x)=h(g(x))$ 那么 $f'(x)=h'(g(x))\cdot g'(x).$

$A_{0}\quad {\xrightarrow {definition\ of\ A}}\ f^{0}(z,c)\quad {\xrightarrow {definition\ of\ f}}\ \quad z$

$B_{0}\quad {\xrightarrow {definition\ of\ B}}\ {\dfrac {\partial }{\partial z}}f^{0}(z,c){\xrightarrow {definition\ of\ f}}\ {\dfrac {\partial }{\partial z}}z{\xrightarrow {derivative}}\ 1$

$C_{0}\quad {\xrightarrow {definition\ of\ C}}\ {\dfrac {\partial }{\partial z}}{\dfrac {\partial }{\partial z}}f^{0}(z,c){\xrightarrow {definition\ of\ f}}\ {\dfrac {\partial }{\partial z}}{\dfrac {\partial }{\partial z}}z{\xrightarrow {derivative}}{\dfrac {\partial }{\partial z}}1{\xrightarrow {derivative}}0$

$D_{0}\quad {\xrightarrow {definition\ of\ D}}{\dfrac {\partial }{\partial c}}f^{0}(z,c){\xrightarrow {definition\ of\ f}}{\dfrac {\partial }{\partial c}}z{\xrightarrow {derivative}}0$

$E_{0}\quad {\xrightarrow {definition\ of\ E}}{\dfrac {\partial }{\partial c}}{\dfrac {\partial }{\partial z}}f^{0}(z,c)\quad {\xrightarrow {definitionoff}}{\dfrac {\partial }{\partial c}}{\dfrac {\partial }{\partial z}}z\quad {\xrightarrow {derivative}}{\dfrac {\partial }{\partial c}}1\quad {\xrightarrow {derivative}}0$

$A_{m+1}\quad {\xrightarrow {definitionofA}}f^{m+1}(z,c)\quad {\xrightarrow {definition\ of\ f}}f(f^{m}(z,c),c)\quad {\xrightarrow {definition\ of\ A}}f(A_{m},c)\quad {\xrightarrow {definition\ of\ f}}A_{m}^{2}+c$

$D_{m+1}\quad {\xrightarrow {definition\ of\ D}}{\dfrac {\partial }{\partial c}}A_{m+1}\quad {\xrightarrow {definition\ of\ A}}{\dfrac {\partial }{\partial c}}(A_{m}^{2}+c)\quad {\xrightarrow {distributivity}}{\dfrac {\partial }{\partial c}}A_{m}^{2}+{\dfrac {\partial }{\partial c}}c\quad {\xrightarrow {derivative}}{\dfrac {\partial }{\partial c}}A_{m}^{2}+1\quad {\xrightarrow {chain\ rule}}2A_{m}({\dfrac {\partial }{\partial c}}A_{m})+1\quad {\xrightarrow {definition\ of\ D}}2A_{m}D_{m}+1$

$E_{m+1}\quad {\xrightarrow {definition\ of\ E}}{\dfrac {\partial }{\partial c}}B_{m+1}\quad {\xrightarrow {definition\ of\ B}}{\dfrac {\partial }{\partial c}}(2A_{m}B_{m})\quad {\xrightarrow {linearity}}2{\dfrac {\partial }{\partial c}}(A_{m}B_{m})\quad {\xrightarrow {product\ rule}}2(A_{m}({\dfrac {\partial }{\partial c}}B_{m})+({\dfrac {\partial }{\partial c}}A_{m})B_{m})\quad {\xrightarrow {definitionof\ E}}2(A_{m}E_{m}+({\dfrac {\partial }{\partial c}}A_{m})B_{m})\quad {\xrightarrow {definition\ of\ D}}2(A_{m}E_{m}+D_{m}B_{m})$

计算

A = 0;
B = 1;
for (m = 0; m < mMax; m++){
  /* first compute B then A */
  B = 2*A*B; /* first derivative with respect to z */
  A = A*A + c; /* complex quadratic polynomial */
   }

或者不使用复数

/* Maxima CAS */
(%i1) a:ax+ay*%i;
(%o1)                                                                                                          %i ay + ax
(%i2) b:bx+by*%i;
(%o2)                                                                                                          %i by + bx
(%i3) bn:2*a*b;
(%o3)                                                                                                 2 (%i ay + ax) (%i by + bx)
(%i4) realpart(bn);
(%o4)                                                                                                      2 (ax bx - ay by)
(%i5) imagpart(bn);
(%o5)                                                                                                      2 (ax by + ay bx)
(%i6) c:cx+cy*%i;
(%o6)                                                                                                          %i cy + cx
(%i7) an:a*a+c;
                                                                                                                                2
(%o7)                                                                                                  %i cy + cx + (%i ay + ax)
(%i8) realpart(an);
                                                                                                                    2     2
(%o8)                                                                                                        cx - ay  + ax
(%i9) imagpart(an);
(%o9)                                                                                                         cy + 2 ax ay

所以

bx = 2 (ax bx - ay by)
by = 2 (ax by + ay bx)
ax = cx - ay^2 + ax^2
ay = cy + 2 ax ay

其中

a = ax + ay*i
b = bx + by*i
c = cx + cy*i

中心

在 GUI 程序中使用它

点击菜单项：核
使用鼠标标记双曲分量中心的矩形区域

一个中心（或核） $c=n$ 周期为 $p$ 满足

f^{p}(0,n)=0

应用单变量牛顿法

$n_{m+1}=N(n_{m})=n_{m}-{\frac {f^{p}(0,n_{m})}{{\dfrac {\partial }{\partial c}}f^{p}(0,n_{m})}}=n_{m}-{\frac {A_{p}(0,n_{m})}{D_{p}(0,n_{m})}}$

其中初始估计值 $n_{0}$ 是任意的（abs(n) <2.0）

停止规则： $A_{p}=0\Leftrightarrow n_{m+1}=n_{m}$

任意名称	初始条件	递推步骤
$A_{m}=f^{m}$	$A_{0}=z=0$	$A_{m+1}=A_{m}^{2}+c$
$D_{m}={\dfrac {\partial }{\partial c}}f^{m}$	$D_{0}=0$	$D_{m+1}=2A_{m}D_{m}+1$

/*

code from  unnamed c program ( book ) 
http://code.mathr.co.uk/book
see mandelbrot_nucleus.c

by Claude Heiland-Allen

COMPILE :

 gcc -std=c99 -Wall -Wextra -pedantic -O3 -ggdb  m.c -lm -lmpfr

usage:

./a.out bits cx cy period maxiters
    
output :  space separated complex nucleus on stdout
    
example

./a.out 53 -1.75 0 3 100
./a.out 53 -1.75 0 30 100
./a.out 53  0.471 -0.3541 14 100
./a.out 53 -0.12 -0.74 3 100
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h> // arbitrary precision
#include <mpfr.h>

extern int mandelbrot_nucleus(mpfr_t cx, mpfr_t cy, const mpfr_t c0x, const mpfr_t c0y, int period, int maxiters) {
  int retval = 0;
  mpfr_t zx, zy, dx, dy, s, t, u, v;
  mpfr_inits2(mpfr_get_prec(c0x), zx, zy, dx, dy, s, t, u, v, (mpfr_ptr) 0);
  mpfr_set(cx, c0x, GMP_RNDN);
  mpfr_set(cy, c0y, GMP_RNDN);

  for (int i = 0; i < maxiters; ++i) {
    // z = 0
    mpfr_set_ui(zx, 0, GMP_RNDN);
    mpfr_set_ui(zy, 0, GMP_RNDN);
    // d = 0
    mpfr_set_ui(dx, 0, GMP_RNDN);
    mpfr_set_ui(dy, 0, GMP_RNDN);
    for (int p = 0; p < period; ++p) {
      // d = 2 * z * d + 1;
      mpfr_mul(u, zx, dx, GMP_RNDN);
      mpfr_mul(v, zy, dy, GMP_RNDN);
      mpfr_sub(u, u, v, GMP_RNDN);
      mpfr_mul_2ui(u, u, 1, GMP_RNDN);
      mpfr_mul(dx, dx, zy, GMP_RNDN);
      mpfr_mul(dy, dy, zx, GMP_RNDN);
      mpfr_add(dy, dx, dy, GMP_RNDN);
      mpfr_mul_2ui(dy, dy, 1, GMP_RNDN);
      mpfr_add_ui(dx, u, 1, GMP_RNDN);
      // z = z^2 + c;
      mpfr_sqr(u, zx, GMP_RNDN);
      mpfr_sqr(v, zy, GMP_RNDN);
      mpfr_mul(zy, zx, zy, GMP_RNDN);
      mpfr_sub(zx, u, v, GMP_RNDN);
      mpfr_mul_2ui(zy, zy, 1, GMP_RNDN);
      mpfr_add(zx, zx, cx, GMP_RNDN);
      mpfr_add(zy, zy, cy, GMP_RNDN);
    }
    // check d == 0
    if (mpfr_zero_p(dx) && mpfr_zero_p(dy)) {
      retval = 1;
      goto done;
    }

    
    // st = c - z / d
    mpfr_sqr(u, dx, GMP_RNDN);
    mpfr_sqr(v, dy, GMP_RNDN);
    mpfr_add(u, u, v, GMP_RNDN);
    mpfr_mul(s, zx, dx, GMP_RNDN);
    mpfr_mul(t, zy, dy, GMP_RNDN);
    mpfr_add(v, s, t, GMP_RNDN);
    mpfr_div(v, v, u, GMP_RNDN);
    mpfr_mul(s, zy, dx, GMP_RNDN);
    mpfr_mul(t, zx, dy, GMP_RNDN);
    mpfr_sub(zy, s, t, GMP_RNDN);
    mpfr_div(zy, zy, u, GMP_RNDN);
    mpfr_sub(s, cx, v, GMP_RNDN);
    mpfr_sub(t, cy, zy, GMP_RNDN);
    // uv = st - c
    mpfr_sub(u, s, cx, GMP_RNDN);
    mpfr_sub(v, t, cy, GMP_RNDN);
    // c = st
    mpfr_set(cx, s, GMP_RNDN);
    mpfr_set(cy, t, GMP_RNDN);
    // check uv = 0
    if (mpfr_zero_p(u) && mpfr_zero_p(v)) {
      retval = 2;
      goto done;
    }
  } // for (int i = 0; i < maxiters; ++i)

 done: mpfr_clears(zx, zy, dx, dy, s, t, u, v, (mpfr_ptr) 0);

 return retval;
}

void DescribeStop(int stop)
{
  switch( stop )
  {
    case 0:
    printf(" method stopped because i = maxiters\n");
    break;
    
    case 1:
    printf(" method stopped because derivative == 0\n");
    break;
    
    //...
    case 2:
    printf(" method stopped because uv = 0\n");
    break;
    
   }
}

void usage(const char *progname) {
  fprintf(stderr,
    "program finds one center ( nucleus) of hyperbolic component of Mandelbrot set using Newton method"
    "usage: %s bits cx cy period maxiters\n"
    "outputs space separated complex nucleus on stdout\n"
    "example %s 53 -1.75 0 3 100\n",
    progname, progname);
}

int main(int argc, char **argv) {
  
  // check the input 
  if (argc != 6) { usage(argv[0]); return 1; }
  
  // read the values 
  int bits = atoi(argv[1]);
  mpfr_t cx, cy, c0x, c0y;
  mpfr_inits2(bits, cx, cy, c0x, c0y, (mpfr_ptr) 0);
  mpfr_set_str(c0x, argv[2], 10, GMP_RNDN);
  mpfr_set_str(c0y, argv[3], 10, GMP_RNDN);
  int period = atoi(argv[4]);
  int maxiters = atoi(argv[5]);
  int stop;

  //
  stop = mandelbrot_nucleus(cx, cy, c0x, c0y, period, maxiters);

   
  //
  printf(" nucleus ( center) of component with period = %s near c = %s ; %s is : \n ", argv[4], argv[2], argv[3]);
  mpfr_out_str(0, 10, 0, cx, GMP_RNDN);
  putchar(' ');
  mpfr_out_str(0, 10, 0, cy, GMP_RNDN);
  putchar('\n');
  //
  DescribeStop(stop) ;

  // clear memeory 
  mpfr_clears(cx, cy, c0x, c0y, (mpfr_ptr) 0);
  
  // 
  return 0;
}

边界

以 $n$ 为中心的成分的边界周期为 $p$ 的可以由内部角度参数化。

边界点 $c=b$ 在角度 $t$ （以圈为单位测量）满足 2 个方程组

${\begin{cases}f^{p}(w,b)-w&=0\\{\dfrac {\partial }{\partial z}}f^{p}(w,b)-e^{2\pi it}&=0\end{cases}}$

定义两个复变量的函数

${\begin{aligned}g{\begin{pmatrix}z\\c\end{pmatrix}}&={\begin{pmatrix}f^{p}(z,c)-z\\{\dfrac {\partial }{\partial z}}f^{p}(z,c)-e^{2\pi it}\end{pmatrix}}\end{aligned}}$

并在两个变量中应用牛顿法

${\begin{aligned}{\boldsymbol {v}}_{0}={\begin{pmatrix}n\\n\end{pmatrix}}\\J_{g}({\boldsymbol {v}}_{m})({\boldsymbol {v}}_{m+1}-{\boldsymbol {v}}_{m})=-g({\boldsymbol {v}}_{m})\end{aligned}}$

其中

${\begin{aligned}J_{g}={\begin{pmatrix}{\dfrac {\partial }{\partial z}}(f^{p}(z,c)-z)&{\dfrac {\partial }{\partial c}}(f^{p}(z,c)-z)\\{\dfrac {\partial }{\partial z}}({\dfrac {\partial }{\partial z}}f^{p}(z,c)-e^{2\pi it})&{\dfrac {\partial }{\partial c}}({\dfrac {\partial }{\partial z}}f^{p}(z,c)-e^{2\pi it})\end{pmatrix}}\end{aligned}}$

这可以使用递归关系表示为

${\begin{aligned}{\begin{pmatrix}w_{0}\\b_{0}\end{pmatrix}}&={\begin{pmatrix}n\\n\end{pmatrix}}\\{\begin{pmatrix}B_{p}-1&D_{p}\\C_{p}&E_{p}\end{pmatrix}}{\begin{pmatrix}w_{m+1}-w_{m}\\b_{m+1}-b_{m}\end{pmatrix}}&=-{\begin{pmatrix}A_{p}-w_{m}\\B_{p}-e^{2\pi it}\end{pmatrix}}\end{aligned}}$

其中 $\{A,B,C,D,E\}_{p}$ 在 $(w_{m},b_{m})$ 处计算。

例子

"该过程是为一次在某个分量上数值查找一个边界点，你不能一次性求解所有分量的整个边界。所以选择并固定一个特定核 n 和一个特定内角 t "(Claude Heiland-Allen ^[21])

让我们尝试找到边界点 (边界点)

$c=b$

的曼德勃罗集的双曲分量

周期 p=3
角 t=0
中心 (核) c = n = -0.12256116687665+0.74486176661974i

只有一个这样的点。

初始估计 (第一次猜测) 将是该分量的中心 (核) $c=n$ 。此中心 (复数) 将存储在一个包含该核两个副本的复数向量中

${\begin{aligned}{\boldsymbol {v}}_{0}={\begin{pmatrix}n\\n\end{pmatrix}}={\begin{pmatrix}-0.12256116687665+0.74486176661974i\\-0.12256116687665+0.74486176661974i\end{pmatrix}}\end{aligned}}$

在角度为 $t$ （以圈数为单位）的边界点满足 2 个方程组

${\begin{cases}f^{3}(z,c)-z&=0\\{\dfrac {\partial }{\partial z}}f^{3}(z,c)-e^{2\pi i0}&=0\end{cases}}$

所以

${\begin{cases}z^{8}+4cz^{6}+6c^{2}z^{4}+2cz^{4}+4c^{3}z^{2}+4c^{2}z^{2}-z+c^{4}+2c^{3}+c^{2}+c=0\\8z^{7}+24cz^{5}+24c^{2}z^{3}+8cz^{3}+8c^{3}z+8c^{2}z-1=0\end{cases}}$

然后函数 g

${\begin{aligned}g({\boldsymbol {v}})={\begin{cases}z^{8}+4cz^{6}+6c^{2}z^{4}+2cz^{4}+4c^{3}z^{2}+4c^{2}z^{2}-z+c^{4}+2c^{3}+c^{2}+c\\8z^{7}+24cz^{5}+24c^{2}z^{3}+8cz^{3}+8c^{3}z+8c^{2}z-1\end{cases}}\end{aligned}}$

雅可比矩阵 $J$ 是

$J_{g}({\boldsymbol {v}})={\begin{pmatrix}8z^{7}+24cz^{5}+24c^{2}z^{3}+8cz^{3}+8c^{3}z+8c^{2}z-1&4z^{6}+12cz^{4}+2z^{4}+12c^{2}z^{2}+8cz^{2}+4c^{3}+6c^{2}+2c+1\\56z^{6}+120cz^{4}+72c^{2}z^{2}+24cz^{2}+8c^{3}+8c^{2}&24z^{5}+48cz^{3}+8z^{3}+24c^{2}z+16cz\end{pmatrix}}$

$J_{g}^{-1}({\boldsymbol {v}})={\begin{pmatrix}{\frac {24z^{5}+48cz^{3}+8z^{3}+24c^{2}z+16cz}{d}}&{\frac {-4*z^{6}-12*c*z^{4}-2*z^{4}-12*c^{2}*z^{2}-8*c*z^{2}-4*c^{3}-6*c^{2}-2*c-1}{d}}\\{\frac {-56*z^{6}-120*c*z^{4}-72*c^{2}z^{2}-24*c*z^{2}-8*c^{3}-8*c^{2}}{d}}&{\frac {8*z^{7}+24*c*z^{5}+24*c^{2}*z^{3}+8*c*z^{3}+8*c^{3}*z+8*c^{2}*z-1}{d}}\end{pmatrix}}$

其中分母 d 为

$d=(24z^{5}+48cz^{3}+8z^{3}+24c^{2}z+16cz)(8z^{7}+24cz^{5}+24c^{2}z^{3}+8cz^{3}+8c^{3}z+8c^{2}z-1)+(-56z^{6}-120cz^{4}-72c^{2}z^{2}-24cz^{2}-8c^{3}-8c^{2})(4z^{6}+12cz^{4}+2z^{4}+12c^{2}z^{2}+8cz^{2}+4c^{3}+6c^{2}+2c+1)$

然后使用迭代法找到点 c=b 的更好的近似值

${\boldsymbol {v}}_{m+1}={\frac {{\boldsymbol {v}}_{m}-g({\boldsymbol {v}}_{m})}{J_{g}({\boldsymbol {v}}_{m})}}=({\boldsymbol {v}}_{m}-g({\boldsymbol {v}}_{m}))J_{g}^{-1}({\boldsymbol {v}}_{m})$

使用 ${\boldsymbol {v}}_{0}$ 作为起点

${\boldsymbol {v}}_{0}$

${\boldsymbol {v}}_{1}=({\boldsymbol {v}}_{0}-g({\boldsymbol {v}}_{0}))J_{g}^{-1}({\boldsymbol {v}}_{0})$

${\boldsymbol {v}}_{2}=({\boldsymbol {v}}_{1}-g({\boldsymbol {v}}_{1}))J_{g}^{-1}({\boldsymbol {v}}_{1})$

$...$

${\boldsymbol {v}}_{m+1}=({\boldsymbol {v}}_{m}-g({\boldsymbol {v}}_{m}))J_{g}^{-1}({\boldsymbol {v}}_{m})$

Maxima CAS 代码

C++ 代码

/*
cpp code by Claude Heiland-Allen 
http://mathr.co.uk/blog/
licensed GPLv3+

g++ -std=c++98 -Wall -pedantic -Wextra -O3 -o bonds bonds.cc

./bonds period nucleus_re nucleus_im internal_angle
./bonds 1 0 0 0
./bonds 1 0 0 0.5
./bonds 1 0 0 0.333333333333333
./bonds 2 -1 0 0
./bonds 2 -1 0 0.5
./bonds 2 -1 0 0.333333333333333
./bonds 3 -0.12256116687665 0.74486176661974 0
./bonds 3 -0.12256116687665 0.74486176661974 0.5
./bonds 3 -0.12256116687665 0.74486176661974 0.333333333333333

for b in $(seq -w 0 999)
do
  ./bonds 1 0 0 0.$b 2>/dev/null
done > period-1.txt
gnuplot
plot "./period-1.txt" with lines

-------------
for b in $(seq -w 0 999)
do
  ./bonds 3 -0.12256116687665 0.74486176661974 0.$b 2>/dev/null
done > period-3a.txt

gnuplot
plot "./period-3a.txt" with lines
*/
#include <complex>
#include <stdio.h>
#include <stdlib.h>

typedef std::complex<double> complex;

const double pi = 3.141592653589793;
const double eps = 1.0e-12;

struct param {
  int period;
  complex nucleus, bond;
};

struct recurrence {
  complex A, B, C, D, E;
};

void recurrence_init(recurrence &r, const complex &z) {
  r.A = z;
  r.B = 1;
  r.C = 0;
  r.D = 0;
  r.E = 0;
}

/*
void recurrence_step(recurrence &rr, const recurrence &r, const complex &c) {
  rr.A = r.A * r.A + c;
  rr.B = 2.0 * r.A * r.B;
  rr.C = 2.0 * (r.B * r.B + r.A * r.C);
  rr.D = 2.0 * r.A * r.D + 1.0;
  rr.E = 2.0 * (r.A * r.E + r.B * r.D);
}
*/

void recurrence_step_inplace(recurrence &r, const complex &c) {
  // this works because no component of r is read after it is written
  r.E = 2.0 * (r.A * r.E + r.B * r.D);
  r.D = 2.0 * r.A * r.D + 1.0;
  r.C = 2.0 * (r.B * r.B + r.A * r.C);
  r.B = 2.0 * r.A * r.B;
  r.A = r.A * r.A + c;
}

struct matrix {
  complex a, b, c, d;
};

struct vector {
  complex u, v;
};

vector operator+(const vector &u, const vector &v) {
  vector vv = { u.u + v.u, u.v + v.v };
  return vv;
}

vector operator*(const matrix &m, const vector &v) {
  vector vv = { m.a * v.u + m.b * v.v, m.c * v.u + m.d * v.v };
  return vv;
}

matrix operator*(const complex &s, const matrix &m) {
  matrix mm = { s * m.a, s * m.b, s * m.c, s * m.d };
  return mm;
}

complex det(const matrix &m) {
  return m.a * m.d - m.b * m.c;
}

matrix inv(const matrix &m) {
  matrix mm = { m.d, -m.b, -m.c, m.a };
  return (1.0/det(m)) * mm;
}

// newton stores <z, c> into vector as <u, v>

vector newton_init(const param &h) {
  vector vv = { h.nucleus, h.nucleus };
  return vv;
}

void print_equation(const matrix &m, const vector &v, const vector &k) {
  fprintf(stderr, "/  % 20.16f%+20.16fi  % 20.16f%+20.16fi  \\  /  z  -  % 20.16f%+20.16fi  \\  =  /  % 20.16f%+20.16fi  \\\n", real(m.a), imag(m.a), real(m.b), imag(m.b), real(v.u), imag(v.u), real(k.u), imag(k.u));
  fprintf(stderr, "\\  % 20.16f%+20.16fi  % 20.16f%+20.16fi  /  \\  c  -  % 20.16f%+20.16fi  /     \\  % 20.16f%+20.16fi  /\n", real(m.c), imag(m.c), real(m.d), imag(m.d), real(v.v), imag(v.v), real(k.v), imag(k.v));
}

vector newton_step(const vector &v, const param &h) {
  recurrence r;
  recurrence_init(r, v.u);
  for (int i = 0; i < h.period; ++i) {
    recurrence_step_inplace(r, v.v);
  }
  // matrix equation: J * (vv - v) = -g
  // solved by: vv = inv(J) * -g + v
  // note: the C++ variable g has the value of semantic variable -g
  matrix J = { r.B - 1.0, r.D, r.C, r.E };
  vector g = { v.u - r.A, h.bond - r.B };
  print_equation(J, v, g);
  return inv(J) * g + v;
}

int main(int argc, char **argv) {
  if (argc < 5) {
    fprintf(stderr, "usage: %s period nucleus_re nucleus_im internal_angle\n", argv[0]);
    return 1;
  }
  param h;
  h.period  = atoi(argv[1]);
  h.nucleus = complex(atof(argv[2]), atof(argv[3]));
  double angle = atof(argv[4]);
  if (angle == 0 && h.period != 1) {
    fprintf(stderr, "warning: possibly wrong results due to convergence into parent!\n");
  }
  h.bond = std::polar(1.0, 2.0 * pi * angle);
  fprintf(stderr, "initial parameters:\n");
  fprintf(stderr, "period = %d\n", h.period);
  fprintf(stderr, "nucleus = % 20.16f%+20.16fi\n", real(h.nucleus), imag(h.nucleus));
  fprintf(stderr, "bond = % 20.16f%+20.16fi\n", real(h.bond), imag(h.bond));
  vector v = newton_init(h);
  fprintf(stderr, "z =% 20.16f%+20.16fi\n", real(v.u), imag(v.u));
  fprintf(stderr, "c =% 20.16f%+20.16fi\n", real(v.v), imag(v.v));
  fprintf(stderr, "\n");
  fprintf(stderr, "newton iteration:\n");
  vector vv;
  do {
    vv = v;
    v = newton_step(vv, h);
    fprintf(stderr, "z =% 20.16f%+20.16fi\n", real(v.u), imag(v.u));
    fprintf(stderr, "c =% 20.16f%+20.16fi\n", real(v.v), imag(v.v));
    fprintf(stderr, "\n");
  } while (abs(v.u - vv.u) + abs(v.v - vv.v) > eps);
  fprintf(stderr, "bond location:\n");
  fprintf(stderr, "c =% 20.16f%+20.16fi\n", real(v.v), imag(v.v));
  // suitable for gnuplot
  printf("%.16f %.16f\n", real(v.v), imag(v.v));
  return 0;
}

大小

牛顿法（在初始估计足够好且函数表现良好的情况下）提供越来越精确的近似值。如何判断近似值足够好？一种方法是比较中心点与其边界上的点，并不断提高精度，直到此距离达到足够位数的精度。

算法

给定一个中心位置估计 $n_{m}$ 精确到 $P$ 位精度
使用中心点 $n_{m}$ 作为起始值，计算边界位置估计 $b_{m}$ ，并确保其精度与中心点相同。
计算 $|n_{m}-b_{m}|$ $|n_{m}-b_{m}|$ 以找到对该组件大小的估计。
1. 如果它不为零，并且它足够精确（可能为几十位），则结束。
2. 否则提高精度，将中心估计值细化到新的精度，从顶部重新尝试。

可以通过比较浮点数的指数和浮点数尾数的精度来测量两个非常接近的点之间距离的有效精度。

内部坐标

Mandelbrot 集双曲分量的内部坐标

https://mathr.co.uk/blog/2013-04-01_interior_coordinates_in_the_mandelbrot_set.html

可以将映射双曲分量 $\mathrm {H}$ 到以原点为中心的闭合单位圆盘 ${\bar {D}}_{1}$

$\lambda _{p}:\mathrm {H} _{p}\to {\bar {D}}_{1}$

$\lambda _{p}(c)=b$

${\bar {D}}_{1}=\{b\in \mathbb {C} :\left|b\right|\leq 1\}$ .

这种关系由两个方程组描述

 ${\begin{cases}F^{p}(z,c)=z\\{\frac {\partial }{\partial z}}F^{p}(z,c)=b\end{cases}}$

其中

p 是参数平面上目标双曲分量的周期
$\theta$ 是点 c 和 b 的所需内部角度
$r\leq 1$ 是内部半径
$b=r*e^{2\pi i\theta }$ 是单位圆上的一个点 = 内部坐标
c 是参数平面上的一个点
z 是动态平面上的一个点
$(F^{0}(z,c)=z$ 以及 $F^{q+1}(z,c)=F^{q}(F(z,c)^{2}+c)$

Claude Heiland-Allen 从 c 中找到内部坐标 b 的算法

选择 c
检查 c (动态平面上进行逃逸测试)。当 c 在 Mandelbrot 集之外时，现在放弃（或计算外部坐标）；
从周期一开始： $p=1$
当 $p<p_{max}$ $p<p_{max}$
- 找到周期点 $z_{0}$ 使得 $F^{p}(z_{0},c)=z_{0}$ 使用一个复变量中的牛顿方法；
- 通过计算 ${\frac {\partial }{\partial z}}F^{p}(z,c)$ 在 $z_{0}$ 处找到 b。
- 如果 $\left|b\right|\leq 1$ ，则返回 b。
- 否则继续下一个 p : $p=p+1$

要解决这样的方程组，可以使用牛顿法。

内部射线

/*
http://code.mathr.co.uk/mandelbrot-graphics/blob/HEAD:/c/bin/m-subwake-diagram-a.c
by Claude Heiland-Allen
*/

#include <mandelbrot-graphics.h>
#include <mandelbrot-numerics.h>
#include <mandelbrot-symbolics.h>
#include <cairo.h>

const double twopi = 6.283185307179586;

void draw_internal_ray(m_image *image, m_d_transform *transform, int period, double _Complex nucleus, const char *angle, double pt, m_pixel_t colour) {
  int steps = 128;
  mpq_t theta;
  mpq_init(theta);
  mpq_set_str(theta, angle, 10);
  mpq_canonicalize(theta);
  double a = twopi * mpq_get_d(theta);
  mpq_clear(theta);
  double _Complex interior = cos(a) + I * sin(a);

  double _Complex cl = 0, cl2 = 0;
  double _Complex c = nucleus;
  double _Complex z = c;
  cairo_surface_t *surface = m_image_surface(image);
  cairo_t *cr = cairo_create(surface);
  cairo_set_source_rgba(cr, m_pixel_red(colour), m_pixel_green(colour), m_pixel_blue(colour), m_pixel_alpha(colour));
  for (int i = 0; i < steps; ++i) {
    if (2 * i == steps) {
      cl = c;
    }
    if (2 * i == steps + 2) {
      cl2 = c;
    }
    double radius = (i + 0.5) / steps;
    m_d_interior(&z, &c, z, c, radius * interior, period, 64);
    double _Complex pc = c;
    double _Complex pdc = 1;
    m_d_transform_reverse(transform, &pc, &pdc);
    if (i == 0) {
      cairo_move_to(cr, creal(pc), cimag(pc));
    } else {
      cairo_line_to(cr, creal(pc), cimag(pc));
    }
  }
  cairo_stroke(cr);
  if (a != 0) {
    double t = carg(cl2 - cl);
    cairo_save(cr);
    double _Complex dcl = 1;
    m_d_transform_reverse(transform, &cl, &dcl);
    cairo_translate(cr, creal(cl), cimag(cl));
    cairo_rotate(cr, -t);
    cairo_translate(cr, 0, -pt/3);
    cairo_select_font_face(cr, "LMSans10", CAIRO_FONT_SLANT_NORMAL, CAIRO_FONT_WEIGHT_NORMAL);
    cairo_set_font_size(cr, pt);
    cairo_text_path(cr, angle);
    cairo_fill(cr);
    cairo_restore(cr);
  }
  cairo_destroy(cr);
}

外部射线

动态射线

//  qmnplane.cpp by Wolf Jung (C) 2007-2014.
// part of Mandel 5.11 http://mndynamics.com/indexp.html
// the GNU General   Public License

//Time ~ nmax^2 , therefore limited  nmax .
int QmnPlane::newtonRay(int signtype, qulonglong N, qulonglong D,
   double &a, double &b, int q, QColor color, int mode) //5, white, 1
{  uint n; int j, i, k, i0, k0; mndAngle t; t.setAngle(N, D, j);
   double logR = 14.0, x, y, u;
   u = exp(0.5*logR); y = t.radians();
   x = u*cos(y); y = u*sin(y);
   if (pointToPos(x, y, i0, k0) > 1) i0 = -10;
   stop(); QPainter *p = new QPainter(buffer); p->setPen(color);
   for (n = 1; n <= (nmax > 5000u ? 5000u : nmax + 2); n++)
   {  t.twice();
      for (j = 1; j <= q; j++)
      {  if (mode & 4 ? tricornNewton(signtype, n, a, b, x, y,
                exp(-j*0.69315/q)*logR, t.radians())
           : rayNewton(signtype, n, a, b, x, y,
                exp(-j*0.69315/q)*logR, t.radians()) )
         { n = (n <= 20 ? 65020u : 65010u); break; }
         if (pointToPos(x, y, i, k) > 1) i = -10;
         if (i0 > -10)
         {  if (i > -10) { if (!(mode & 8)) p->drawLine(i0, k0, i, k); }
            else { n = 65020u; break; }
         }
         i0 = i; k0 = k;
      }
   }
   //if rayNewton fails after > 20 iterations, endpoint may be accepted
   delete p; update(); if (n >= 65020u) return 2;
   if (mode & 2) { a = x; b = y; }
   if (n >= 65010u) return 1; else return 0;
} //newtonRay

int QmnPlane::rayNewton(int signtype, uint n, double a, double b,
   double &x, double &y, double rlog, double ilog)
{  uint k, l; double fx, fy, px, py, u, v = 0.0, d = 1.0 + x*x + y*y, t0, t1;
   for (k = 1; k <= 60; k++)
   {  if (signtype > 0) { a = x; b = y; }
      fx = cos(ilog); fy = sin(ilog);
      t0 = exp(rlog)*fx - 0.5*exp(-rlog)*(a*fx + b*fy);
      t1 = exp(rlog)*fy + 0.5*exp(-rlog)*(a*fy - b*fx);
      fx = x; fy = y; px = 1.0; py = 0.0;
      for (l = 1; l <= n; l++)
      {  u = 2.0*(fx*px - fy*py); py = 2.0*(fx*py + fy*px);
         px = u; if (signtype > 0) px++;
         u = fx*fx; v = fy*fy; fy = 2.0*fx*fy + b; fx = u - v + a;
         u += v; v = px*px + py*py; if (u + v > 1.0e100) return 1;
      }
      fx -= t0; fy -= t1; if (v < 1.0e-50) return 2;
      u = (fx*px + fy*py)/v; v = (fx*py - fy*px)/v;
      px = u*u + v*v; if (px > 9.0*d) return -1;
      x -= u; y += v; d = px; if (px < 1.0e-28 && k >= 5) break;
   }
   return 0;
} //rayNewton

参数射线

外部参数射线

等势线

势或逃逸时间水平集的边界（等势曲线）是多项式勒姆尼斯卡特。从符号方程计算出来的边界越接近曼德勃罗集或朱利亚集，其复杂度就呈指数级增长。

是否有更有效的方法来计算曼德勃罗集勒姆尼斯卡特？^[22]

$z\to z^{2}+c$ 的复迭代在 $O(n)$ 的工作量内，在 $c$ 中产生一个 $2^{n-1}$ 次多项式。

假设

c=x+iy

f_{c}(z)=z^{2}+c

f_{c}^{\circ (n+1)}(z)=f_{c}^{\circ n}(f_{c}(z))

那么

{\frac {d}{dc}}f_{c}^{\circ (n+1)}(z)=2f_{c}^{\circ n}(z){\frac {d}{dc}}f_{c}^{\circ n}(z)+1

这些可以在一个内部循环中一起计算（小心不要覆盖仍然需要的旧值）。

现在可以使用牛顿求根法来求解隐式形式 $f_{c}^{\circ n}(0)=re^{2\pi i\theta }=t$ 。

c_{m+1}=c_{m}-{\frac {f_{c_{m}}^{\circ n}(0)-t}{{\frac {d}{dc}}f_{c_{m}}^{\circ n}(0)}}

使用之前的 $c$ 作为下一个 $\theta$ 的初始猜测。 $\theta$ 的增量需要小于大约 ${\frac {1}{4}}$ ，才能使算法稳定。注意， $\theta$ （以圈数为单位测量）在 $c$ 回到起点之前，会在单位圆上绕 $2^{n-1}$ 圈。

这种方法的灵感来自于川平知樹的《Mandelbrot 集合外射线绘制算法》。

周期点

另请参阅

参考文献

[1] NEWTON’S METHOD AND FRACTALS by AARON BURTON

[2] Newton, Chebyshev, and Halley Basins of Attraction; A Complete Geometric Approach by Bart D. Stewart

[3] "How to find all roots of complex polynomials by Newton's method", by Hubbard, Schliecher, and Sutherland

[4] Finding the all roots of a polynomial by using Newton-Raphson method.

[5] Everything You Always Wanted to Ask About Newton's Method But Were Afraid to Know by Ernst W. Mayer

[6] th.stackexchange question: how-to-solve-simultaneous-equations-using-newton-raphsons-method

[7] Iterative Methods for Sparse Linear Systems Luca Bergamaschi

[8] Newton's method for nonlinear systems by Mark S. Gockenbach 2003-01-23

[9] wikipedia : 雅可比矩阵和行列式

[10] Iterative Methods for Sparse Linear Systems Luca Bergamaschi

[11] Multiple Nonlinear Equations using the Newton-Raphson Method by Bruce A. Finlayson

[12] Stopping Criterion by Prof. Alan Hood

[13] th.stackexchange question : why-does-the-newton-raphson-method-not-converge-for-some-functions

[14] Numerical Methods for Scientists and Engineers by Richard Hamming, page 69

[15] StackExchange : A function for which the Newton-Raphson method slowly converges?

[16] stackoverflow question: how-to-detect-a-cycle-in-newton-raphson-method-in-c

[17] Numerical Recipies In C

[18] Numerical_Methods_in_Civil_Engineering : 非线性方程 Matlab By Gilberto E. Urroz, September 2004

[19] Solving nonlinera equations with Scilab by G E Urroz

[20] Principles of Linear Algebra With Mathematica® 牛顿-拉夫森法 Kenneth Shiskowski and Karl Frinkle

[21] Home page of Claude Heiland-Allen

[22] th.stackexchange qestion: is-there-a-more-efficient-equation-for-the-mandelbrot-lemniscates

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]