## Taylor Series Fractals

The taylor series is a mechanism for approximating a function by summing polynomial terms.

First, some definitions:

 f(x) the function to approximate fn(x) the nth derivative of f (note that f0(x) = f(x)) n! n factorial, n! = n * (n - 1) * (n - 2) * ... * 2

The taylor series expansion of f(x) around a point (c) is then written:

```
n -> infinity
+----
\      fn(c) * (x - c) n
f(x) = +      -----------------
/             n!
+----
n = 0
```

As an example, consider the case where:

 f(x) => exp(x) c => 0

Since the derivative of exp(x) is exp(x) we know that fn(x) = exp(x) for all n.

Here are the first few terms of the series:

value of n term value of term
0 exp(0) * (x - 0)0 / 0! 1
1 exp(0) * (x - 0)1 / 1! x
2 exp(0) * (x - 0)2 / 2! x2 / 2
3 exp(0) * (x - 0)3 / 3! x3 / 6
4 exp(0) * (x - 0)4 / 4! x4 / 24

We can write the expansion of exp(x) around c = 0 more concisely as follows:

```        n -> infinity
+----
\      xn
f(x) = +      --
/      n!
+----
n = 0
```

Now lets look at how accurate the approximation is as we add more terms. In this case, lets use x = 3.

exp(3) n sum approximation difference
20.08553692318766774092 0 1 1 19.08553692318766774092
20.08553692318766774092 1 1 + x 4 16.08553692318766774092
20.08553692318766774092 2 1 + x + x2/2 8.5 11.58553692318766774092
20.08553692318766774092 3 1 + x + x2/2 + x3/6 13 7.08553692318766774092
20.08553692318766774092 4 1 + x + x2/2 + x3/6 + x4/24 16.375 3.71053692318766774092

As you can see, the approximation converges pretty rapidly on the actual value of exp(x).

Ok, so how do we apply this to generating a fractal? Firstly, all the statements about 'x' above apply equally well to 'z' where z is the complex number (x + iy).

In our case, assume 'current' is the point 'z' for which we wish to create an approximation of exp(z). Our algorithm (in pseudo code) looks like this:

``` set actualValue   = exp(current)
set maxIterations = 20
set tolerance     = .00001
set count         = 0
set zSum          = [0, 0]
set denominator   = 1

while count is less than maxIterations and
the magnitude of the difference between actualValue and zSum is
greater than tolerance
do
zSum  = zSum + currentcount / denominator

if count is greater than 0
denominator = denominator * (count + 1)
endif

count = count + 1
done
```

In this example, our expansion point (c) is zero, we can readily generalize our algorithm to handle other values of (c)

``` set actualValue   = exp(current)
set maxIterations = 20
set tolerance     = .00001
set count         = 0
set zSum          = [0, 0]
set denominator   = 1
set c             = [2, 3]
set cValue        = exp(c)

while count is less than maxIterations and
the magnitude of the difference between actualValue and zSum is
greater than tolerance
do
zSum  = zSum + cValue * (current - c)count / denominator

if count is greater than 0
denominator = denominator * (count + 1)
endif

count = count + 1
done
```

Finally, here is a full-on fractorama formula file demonstrating this technique

```fractal
{
mapping {
(-2.33672316384180778215,
3.41009999999999990905,
1.18621468926553697010,
8.00000000000000000000) => (150, 195)
}

formula
{
// The value we wish to approximate

zvalue = exp(current);

\$xavg = 0;
\$yavg = 0;
\$mavg = 0;

\$xmax = 0;
\$ymax = 0;
\$mmax = 0;

\$xtot = 0;
\$ytot = 0;
\$mtot = 0;

\$maxIterations = 15;
\$tolerance     = 1e-6;

zSum           = 0;
\$denominator   = 1;

// The point around which we'll create our expansion.  We use
// exp(current) here but we could just as well have used:
//
// point = [2, 3]
// point = sin(current)
// etc.
//
// There is no special significance to the fact that we're using
// exp(current) and we're attempting to approximate exp.

point = exp(current);

// Loop until we've been through the loop '\$maxIterations' times
// or we deem that zSum is close enough to zvalue

while(\$count < \$maxIterations && ssq(zvalue - zSum) > \$tolerance)
{
zSum += exp(point) * ((current - point) ^ \$count) / \$denominator;

if(\$count > 0) { \$denominator *= (\$count + 1); }

value = cot(zSum);

\$x = mag(value ^ .1);
\$y = mag(value ^ .2);
\$m = mag(value ^ .4);

\$xavg = (\$x + \$xavg * \$count) / (\$count + 1);
\$yavg = (\$y + \$yavg * \$count) / (\$count + 1);
\$mavg = (\$m + \$mavg * \$count) / (\$count + 1);

\$xdev = sqrt((\$x - \$xavg) * (\$x - \$xavg));
\$ydev = sqrt((\$y - \$yavg) * (\$y - \$yavg));
\$mdev = sqrt((\$m - \$mavg) * (\$m - \$mavg));

\$xmax = max(\$xmax, \$xdev);
\$ymax = max(\$ymax, \$ydev);
\$mmax = max(\$mmax, \$mdev);

\$xtot += \$xdev;
\$ytot += \$ydev;
\$mtot += \$mdev;
}

\$r = 0;
\$g = 0;
\$b = 0;

if(\$count > 0)
{
\$r = (\$xtot/\$count) / \$xmax;
\$g = (\$ytot/\$count) / \$ymax;
\$b = (\$mtot/\$count) / \$mmax;

\$r = get_sin_color(\$r, 255, 1);
\$g = get_sin_color(\$g, 255, 1);
\$b = get_sin_color(\$b, 255, 1);
}

set_color(\$r, \$g, \$b);
}
}
```

And here is the resulting image: 