Fractals/Mathematics/Newton method
types
editNewton method can be used for finding successively better approximations to one root (zero) x of function f :
![](http://upload.wikimedia.org/wikipedia/commons/thumb/e/e0/NewtonIteration_Ani.gif/220px-NewtonIteration_Ani.gif)
If one wants to find another root must apply the method again with other initial point.
How to find all roots ? [3][4]
Newton method can be applied to :
- a real-valued function [5]
- a complex-valued function
- a system of equations
systems of nonlinear equations
editOne may also use Newton's method to solve systems of [6]
- 2 (non-linear) equations
- with 2 complex variables : [7]
It can be expressed more synthetically using vector notation :
where is a vector of 2 complex variables :
and is a vector of functions ( or is a vector-valued function of the vector variable x ) :
Solution can be find by iteration :[8]
where s is an increment ( which is now a vector of 2 components ):
Increment can be computed by :
Here is Jacobian matrix[9] of at :
and is inverse of the Jacobian matrix .
In case of
- a small number of equations use the inverse matrix.
- a big number of equations : rather than actually computing the inverse of this matrix, one can save time by solving the system of linear equations for the unknown xn+1 − xn.
Algorithm in pseudocode :[10]
- start with an initial approximation ( guess)
- repeat until convergence ( iterate)
- solve:
- compute :
Stop criteria :
- absolute error : [11]
What is needed ?
edit- function f ( a differentiable function )
- it's derivative f'
- starting point x0 ( which should be in the basin of root's attraction )
- stopping rule ( criteria to stop iteration )[12]
stopping rule
edit
Pitfalls or failure
editWhy method fails to find a root ?[13]
- method has a cycle and never converge (
- inflection point [14]
- local minimum
- method is diverging not converging
- local minimum
- multiple root ( multiplicity of root > 1 , for example : f(x) = f'(x) = 0). Slow approximation, derivative tends to zero, trouble in division step ( do not divide by zero ). Use modified Newton's method
- slow convergence[15]
- bad value of the initial approximation ( an initial guess far from the root )
- a function whose derivative vanishes near a root, or whose second derivative is unbounded near a root
-
Function has inflection point ( method has a cycle and never converge)
-
Local minimum
"For most problems, the application of Newton's method in double precision has 2-3 steps of global orientation and then 4-6 steps of quadratic convergence before the precision of the double floating point format is exceeded. Thus, one can be even bolder, if after 10 steps the iteration does not converge, the initial point was bad, be it that it leads to a periodic orbit or divergence to infinity. Most likely, close-by initial points will behave similarly, so make a non-trivial change in the initial point for the start of the next run of the iteration."LutzL[16]
Pseudocode
editThe following is an example of using the Newton's Method to help find a root of a function f
which has derivative fprime
.
The initial guess will be and the function will be so that .
Each new iterative of Newton's method will be denoted by x1
. We will check during the computation whether the denominator (yprime
) becomes too small (smaller than epsilon
), which would be the case if , since otherwise a large amount of error could be introduced.
%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
code
edit- c [17]
- matlab : newtonm function [18]
- Maxima CAS
- scilab [19]
- Mathematica[20]
- Animation engine for explanatory math videos in python
error analysis
editapplications of Newton's method in case of iterated complex quadratic polynomials
editrecurrence relations
editA recurrence relation is an equation ( map) that recursively defines a sequence of points (for example orbit). To compute the sequence one needs:
- starting point
- map (equation = function = recurrence relation )
definitions
editDefine the iterated function by:
Please note the other function:
which is used for :
Pick arbitrary names for the iterated function and it's derivatives These derivatives can be computed by recurrence relations.
description | math notation | arbitrary names | Initial conditions | Recurrence steps | common names |
---|---|---|---|---|---|
iterated function | z or f | ||||
first derivative with respect to z | d (from derivative) or p ( from prime) of f' | ||||
second derivative with respect to z | |||||
first derivative with respect to c | |||||
second derivative with respect to c |
Derivation of recurrence relations
editFirst basic rules:
- Iterated function is a function composition
- Applay chain rule for composite functions, so if then
Computing
editA = 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 */
}
or without complex numbers:
/* 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
so :
- bx = 2 (ax bx - ay by)
- by = 2 (ax by + ay bx)
- ax = cx - ay^2 + ax^2
- ay = cy + 2 ax ay
where:
- a = ax + ay*i
- b = bx + by*i
- c = cx + cy*i
center
editTo use it in a GUI program :
- click on the menu item : nucleus
- using a mouse mark rectangle region around center of hyperbolic component
A center ( or nucleus ) of period satisfies :
Applying Newton's method in one variable:
where initial estimate is arbitrary choosen ( abs(n) <2.0 )
stoppping rule :
arbitrary names | Initial conditions | Recurrence steps |
---|---|---|
/*
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;
}
boundary
editThe boundary of the component with center of period can be parameterized by internal angles.
The boundary point at angle (measured in turns) satisfies system of 2 equations :
Defining a function of two complex variables:
and applying Newton's method in two variables:
where
This can be expressed using the recurrence relations as:
where are evaluated at .
Example
edit"The process is for numerically finding one boundary point on one component at a time - you can't solve for the whole boundary of all the components at once. So pick and fix one particular nucleus n and one particular internal angle t " (Claude Heiland-Allen [21])
Lets try find point of boundary ( bond point )
of hyperbolic component of Mandelbrot set for :
- period p=3
- angle t=0
- center ( nucleus ) c = n = -0.12256116687665+0.74486176661974i
There is only one such point.
Initial estimates (first guess) will be center ( nucleus ) of this components. This center ( complex number) will be saved in a vector of complex numbers containing two copies of that nucleus:
The boundary point at angle (measured in turns) satisfies system of 2 equations :
so :
Then function g
Jacobian matrix is
where denominator d is :
Then find better aproximation of point c=b using iteration :
using as an starting point :
Maxima CAS code
editC++ code
edit/*
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;
}
size
editNewton's method (providing the initial estimate is sufficiently good and the function is well-behaved) provides successively more accurate approximations. How to tell when the approximation is good enough? One way might be to compare the center with points on its boundary, and continue to increase precision until this distance is accurate to enough bits.
Algorithm:
- given a center location estimate accurate to bits of precision
- compute a boundary location estimate accurate to the same precision, using center as starting value
- compute to find an estimate of the size of the component
- if it isn't zero, and if it is accurate enough (to a few 10s of bits, perhaps) then finish
- otherwise increase precision, refine center estimate to new precision, try again from the top
Measuring effective accuracy of the distance between two very close points might be done by comparing floating point exponents with floating point mantissa precision.
internal coordinate
editInternal coordinate of hyperbolic component of Mandelbrot set
It is possible to map hyperbolic component to closed unit disk centered at origin :
.
This relation is described by system of 2 equations :
where
- p is the period of the target hyperbolic component on parameter plane
- the desired internal angle of points c and b
- is an internal radius
- is a point of unit circle = internal coordinate
- c is a point of parameter plane
- z is a point of dynamic plane
- and
The algorithm by Claude Heiland-Allen for finding internal coordinate b from c :
- choose c
- check c ( bailout test on dynamical plane). When c is outside the Mandelbrot set, give up now ( or compute external coordinate );
- Start with period one :
- while
- Find periodic point such that using Newton's method in one complex variable;
- Find b by evaluating at
- If then return b
- otherwise continue with the next p :
To solve such system of equations one can use Newton method.
internal ray
edit/*
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);
}
external rays
editDynamic ray
edit// 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
Parameter ray
editequipotential line
editBoundaries of potential or escape time level sets ( equipotential curves) are polynomial lemniscates. Computing from symbolic equations gets exponentially more complicated as boundary gest closer to Mandelbrot or Julia set.
Is there a more efficient way of computing for the Mandelbrot lemniscates?[22]
Complex iteration of gives a degree polynomial in in work.
Suppose
Then
These can be calculated together in one inner loop (being careful not to overwrite old values that are still needed).
Now you can use Newton's root finding method to solve the implicit form by
Use the previous as initial guess for next . The increment in needs to be smaller than about for the algorithm to be stable. Note that (measured in turns) wraps around the unit circle times before gets back to its starting point.
This approach is inspired by An algorithm to draw external rays of the Mandelbrot set by Tomoki Kawahira
Periodic points
editSee also
editReferences
edit- ↑ NEWTON’S METHOD AND FRACTALS by AARON BURTON
- ↑ Newton, Chebyshev, and Halley Basins of Attraction; A Complete Geometric Approach by Bart D. Stewart
- ↑ "How to find all roots of complex polynomials by Newton's method", by Hubbard, Schliecher, and Sutherland
- ↑ Finding the all roots of a polynomial by using Newton-Raphson method.
- ↑ Everything You Always Wanted to Ask About Newton's Method But Were Afraid to Know by Ernst W. Mayer
- ↑ math.stackexchange question: how-to-solve-simultaneous-equations-using-newton-raphsons-method
- ↑ Iterative Methods for Sparse Linear Systems Luca Bergamaschi
- ↑ Newton's method for nonlinear systems by Mark S. Gockenbach 2003-01-23
- ↑ wikipedia : Jacobian_matrix_and_determinant
- ↑ Iterative Methods for Sparse Linear Systems Luca Bergamaschi
- ↑ Multiple Nonlinear Equations using the Newton-Raphson Method by Bruce A. Finlayson
- ↑ Stopping Criterion by Prof. Alan Hood
- ↑ math.stackexchange question : why-does-the-newton-raphson-method-not-converge-for-some-functions
- ↑ Numerical Methods for Scientists and Engineers by Richard Hamming, page 69
- ↑ StackExchange : A function for which the Newton-Raphson method slowly converges?
- ↑ stackoverflow question: how-to-detect-a-cycle-in-newton-raphson-method-in-c
- ↑ Numerical Recipies In C
- ↑ Numerical_Methods_in_Civil_Engineering : NonLinear Equations Matlab By Gilberto E. Urroz, September 2004
- ↑ Solving nonlinera equations with Scilab by G E Urroz
- ↑ Principles of Linear Algebra With Mathematica® The Newton–Raphson Method Kenneth Shiskowski and Karl Frinkle
- ↑ Home page of Claude Heiland-Allen
- ↑ math.stackexchange qestion: is-there-a-more-efficient-equation-for-the-mandelbrot-lemniscates
- NEWTON’S METHOD IN PRACTICE: FINDING ALL ROOTS OF POLYNOMIALS OF DEGREE ONE MILLION EFFICIENTLY by DIERK SCHLEICHER AND ROBIN STOLL
- NEWTON’S METHOD IN PRACTICE II: THE ITERATED REFINEMENT NEWTON METHOD AND NEAR-OPTIMAL COMPLEXITY FOR FINDING ALL ROOTS OF SOME POLYNOMIALS OF VERY LARGE DEGREES by DIERK SCHLEICHER AND ROBIN STOLL
- Finding polynomial roots by dynamical systems -- a case studyby Sergey Shemyakov, Roman Chernov, Dzmitry Rumiantsau, Dierk Schleicher, Simon Schmitt, Anton Shemyakov