型別

牛頓法可用於找到函式 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]