# perturbation and aproximation

* how solution of the problem is altered for non-zero but small perturbation ( usually denoted as an espilon)
* problem can not be solved explicitely ( using symbolic methods), only numericla aproximation can be found

# Convergence

" ... methods of acceleration of convergence. Suppose you have a slowly convergent series, and want to know its sum (numerically). Just by summing x_1 + x_2 + ... + x_{1000} + ... x_{1000000} + ... + x_{1000000000} + ... you will get the required accuracy after 100 years. If you fit your x_n to c_2/n^2 + c_3/n^3 + (a few more terms), you will get the same accuracy of the sum in 1 second." Andrey (theoretical high energy physicist)
Above problem you can find in parabolic Julia sets.
" your sequence was generated as ... the output of a fixed point iteration $x_{i+1}=f(x_{i})$. ... The Levin transformation ... uses forward differences successively to remove error terms of an alternating series."
Examples:
* GSL library
* "checking for when the series approximation starts to deviate from the result of running it through the Levin u-transform available in the gnu scientific library. so far in my tests this is taking me right up to the event horizon of what will end up being the lowest iteration, even in locations that used to cause the other methods to bail on the series approximation way sooner, but it also does not seem to be going too far so as to cause glitchy artifacts. it seems to be just right, across the board, which i find to be rather exciting! " ( quaz0r about Mandelbrot )
* series calculator

# limitation of native double numbers format

* when 2 pixels separated by less than :$2^{-53}$
* when the size of a pixel drops below approx $10^{-308}$

# Precision

Precision in numerical compuiting
* Floating-point arithmetic Accuracy_problems in wikipedia
How to deal with lack of floating point precision:
* implement higher precision arithmetic than your hardware natively supports
** software emulation ( emulating a double with two floats, fixed points numbers , ...)
* use algorithms that are more numerically stable
* use representable numbers
how-many-floating-point-numbers-are-in-the-interval [0,1] ?
Libraries :
* Arb by Fredrik Johansson
See how arb library can be used for checking and adjusting precision :

## How many decimal digits are there in n bits ?

One digit of binary number needs one bit : there are 2 binary digits ( 0 and 1) and bits have 2 states. One digit of decimal number needs approximately 3.4 bits. There are 10 decimal numbers, 3 bits have 8 states which is not enough, 4 bits have 16 states which is too much.
| type | Total bits | Bits precision | Number of significant decimal digits |
| float | 32 | 24 | 6 |
| double | 64 | 53 | 15 |
| long double | 80 | 64 | 18 |
| mpfr_t | | 142 | 41 |
| mpfr_t | | 1000 | 294 |
See also: decimal-precision-of-binary-floating-point-numbers by RIck Regan

## Rules
### Rounding

* "compute with two more significant figures than your ultimate answer
* Round only after the last step in calculation. Never do further calculations with rounded numbers.
* for multiplying and dividing find the number of significant digits in each factor. The result will have the smaller number of significant digits.
* For powers and roots, the answer should have the same number of significant digits as the original"
* "When you add (or subtract), you keep as many decimal places as there are in the least accurate number.
* When you multiply (or divide), you keep as many significant digits as there are in the least accurate number."
See
* gcc -frounding-math option

# Examples of numerical computings

* Finding roots of equation

## Seed value

* fractalforums.org: julia-sets-true-shape-and-escape-time

## standard images ( without zoom)

* Roundoff Error
* Shadowing lemma : that every pseudo-orbit (which one can think of as a numerically computed trajectory with rounding errors on every step) stays uniformly close to some true trajectory

## External rays
### Parameter rays

The Wolf Jung test : The external parameter rays for angles (in turns)
* 1/7 (period 3)
* 321685687669320/2251799813685247 (period 51)
* 321685687669322/2251799813685247 ( period 51 )
Angles differ by about $10^{-15}$, but the landing points of the corresponding parameter rays are about 0.035 apart.

### Dynamic rays
#### angle precision

For rotational number ( internal angle) 1/34 ray for external angle :
$t_{a}={\frac {2^{33}}{2^{34}-1}}={\frac {8589934592}{17179869183}}=0.5000000000291038$
lands on the alfa fixed point :
$z_{a}=0.491486549841951+0.091874758908285i$
It is not a ray for angle :
$t_{b}={\frac {1}{2}}=0.5$
which land on the point :
$z_{b}=-0.508513450158049+0.091874758908285i$
Difference between external angles of the rays is :
$dt=t_{b}-t_{a}=-0.0000000000291038$
and between landing points of the rays points is :

## Escaping test

This test was introduced by by John Milnor. See also analysis by Mark Braverman and roundoff error by Robert P. Munafo
Comment by Mark McClure : " an escape time algorithm would take forever to generate that type of image, since the dynamics are so slow there. If you want resolution of 1/100, it would take roughly 2*10^8 iterates to move the point z0=0.01 to z=2 by iterating f(z)=z+z5 "
Here's an orbit computed with 1000 decimal digits, escaping after about 300000 iterations by gerrit

### Cases
#### Nonparabolic case

Lets take simple hyperbolic case where parameter c is :
Here repelling fixed point z_f is :

#### Parabolic case

Lets take simple parabolic case where parameter c is :
Here parabolic fixed point z_f is :

### Test

Lets take point z of exterior of Julia set but lying near fixed point :
where n is a positive integer
Check how many iterates i needs point z to reach target set ( = escape) :
where ER is Escape Radius
Show relationship between :
* n
* Last Iteration
* type of numbers used for computations ( float, double, long double, extended, arbitrary precision )

### Programs

See FractalForum for evaldraw script

### Results

Maximal n for which program does not fall into fixed point (columns = number types (precision of significand); rows = case) :
| | float (24) | double (53) | long double (64) | MPFR (80) | MPFR (100) |
| hyperbolic | 23 | 52 | 63 | 79 | 99 |
| parabolic | 12 | 26 | 32 | | |
The results for standard C types ( float, double, long double) and MPFR precision are the same for the same precision
Relation between : Last iteration, n ( in columns) and case ( in rows) :
| | 1 | 2 | 3 | 4 | 5 | 24 | 53 | 64 | 80 | 100 |
| hyperbolic | 1 | 2 | 3 | 4 | 5 | 24 | 53 | 64 | 80 | 100 |
| parabolic | 3 | 5 | 10 | 19 | 35 | 16 778 821 | | | | |
Relation between number of iterations and time of computation in hyperbolic case :
Parabolic case :

### Analysis

Maximal n in hyperbolic case $n=n_{h}$ is almost the same as the precision of the significand
In parabolic case maximal $n=n_{p}$ is
Last Iteration ( escape time = iteration fro which abs(zn) > ER ) is : in hyperbolic case equal to n :
in parabolic case equal to 2^n :
Time of computations is proportional to number of iterations. In hyperbolic case is is short. In parabolic case grows quickly as number of iterations.
Checking one point if escapes in parabolic case :
* for n = 34 take about one hour ( 5 595 seconds )
* for n = 40 take about one day
* for n = 45 take about one month
* for n = 50 take about one year

### Q&A
#### Why programs fails ?

Cancellation of significant digits and loss of significance (LOS).
The program fails because of limited precision of used number types. Addition of big (zp) and small number (distance) gives number which has more decimal digits then can be saved ( floating point type has only 7 decimal digits). Some of the most right digits are cancelled and iteration goes into an infinite loop.
For example : when using floating point in parabolic case lets take
so
It is greater than machine epsilon:
so this addition still works :
After multiplication it gives :
next step is addition. Because floating point format saves only 7 decimal digits it is truncated to :
Here relative error is to big and
is smaller then FLT_EPSILON/2.0 = 0.0000000596;
Solution : increase precision !
Explanation in polish

#### Why parabolic dynamics is so weak ?

Parabolic dynamics
* numerical problems
* Distance to fixed point for various types of local discrete dynamics near fixed point . See how parabolic dynamics is complicated
* Orbits near fixed point of fat Douady rabbit Julia set. See that to analyze one orbit one have to include period and petal

## Distance estimation

Example

## Zoom

* manipulating numbers below $10^{-308}$ in double format shall be avoided.
What precision do I need for zoom ?
* Pixel density
* Pixel spacing is a distance between the centers of each two-dimensional pixel
One can check it using such program ( automatic math precision) :

## Tent map

Numerical error in computing orbit of the tent map with parameter m = 2. Aftere 50 iterations of double numbers all orbits fall into zero.
Result:

# See also

* A few notes on computation by Mark McClure
* The Numerical Tours of Data Sciences, by Gabriel Peyré, gather Matlab, Python, Julia and R experiments to explore modern mathematical data sciences
* Numerical Errors in Logistic Map Calculation
* Keisan Online Calculator
* Finite arithmetic and error analysis in python by Gonzalo Galiano Casas and Esperanza García Gonzalo
* math.stackexchange question: is-there-an-equation-for-the-number-of-iterations-required-to-escape-the-mandelb ?

# References

* ↑ mathoverflow question : /rounding-errors-in-images-of-julia-sets
* ↑ ChaosAndFractals by Mark McClure. 2.11 A few notes on computation
* ↑ limit of series at gmane.org discussion
* ↑ math.stackexchange questions : levins-u-transformation
* ↑ math.stackexchange questions : accelerating-convergence-of-a-sequence
* ↑ gsl-1.0 : Series Acceleration
* ↑ gsl manual : Series-Acceleration
* ↑ Fractal Forum : series-acceleration
* ↑ wolframalpha series calculator
* ↑ Fractalshades documentation
* ↑ What you never wanted to know about floating point but will be forced to find out by volkerschatz
* ↑ stackoverflow : Dealing with lack of floating point precision in OpenCL particle system
* ↑ Heavy computing with GLSL – Part 5: Emulated quadruple precision by Henry Thasler
* ↑ how-many-floating-point-numbers-are-in-the-interval [0,1 ? by Daniel Lemire]
* ↑ Arb by Fredrik Johansson
* ↑ Computing π with Chudnovsky and GMP by Beej Jorgensen
* ↑ [:w:Extended precision: Extended precision in wikipedia]
* ↑ [:w:GNU MPFR|GNU MPFR in wikipedia]
* ↑ Significant Digits and Rounding Copyright © 2003–2014 by Stan Brown
* ↑ mathforum : Rules for Significant Figures and Decimal Places
* ↑ Roundoff Error Robert P. Munafo, 1996 Dec 3.
* ↑ wikipedia : Shadowing_lemma
* ↑ Wolf Jung's test for precision of drawing parameter rays
* ↑ One can find it using program Mandel by Wolf Jung using Ray to point menu position ( or Y key)
* ↑ Dynamics in one complex variable: introductory lectures version of 9-5-91 Appendix G by John W. Milnor
* ↑ Parabolic Julia Sets are Polynomial Time Computable Mark Braverman
* ↑ Roundoff Error by Robert P. Munafo, 1996 Dec 3.
* ↑ math.stackexchange questions : what-is-the-shape-of-parabolic-critical-orbit
* ↑ Parabolic Julia Sets are Polynomial Time Computable by Mark Braverman
* ↑ fractalforums : numerical-problem-with-bailout-test
* ↑ wikipedia :Floating point , EEE_754
* ↑ wikipedia : Significant figures
* ↑ wikipedia : Loss of significance
* ↑ IEEE Arithmetic from Numerical Computation Guide by Oracle
* ↑ Machine epsilon
* ↑ Dyskusja po polsku na pl.comp.os.linux.programowanie
* ↑ Precision and mandel zoom using DEM/M
* ↑ Fractalshades-doc : math
* ↑ reenigne blog : arbitrary precision mandelbrot sets
* ↑ hpdz : Bignum by Michael Condron
* ↑ fractint : Arbitrary Precision and Deep Zooming
* ↑ chaospro documentation : parmparm
* ↑ Pixel density
* ↑ Automatic Math-Precision Robert P. Munafo, 1993 Apr 15.
* Do Numerical Orbits of Chaotic Dynamical Processes Represent True Orbits?