Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

Vector form of Julia fractal

Posted by peeterjoot on December 27, 2010

[Click here for a PDF of this post with nicer formatting]

Motivation.

As outlined in [1], 2-D and N-D Julia fractals can be computed using the geometric product, instead of complex numbers. Explore a couple of details related to that here.

Guts

Fractal patterns like the mandelbrot and julia sets are typically using iterative computations in the complex plane. For the Julia set, our iteration has the form

\begin{aligned}Z \rightarrow Z^p + C\end{aligned} \hspace{\stretch{1}}(2.1)

where p is an integer constant, and Z, and C are complex numbers. For p=2 I believe we obtain the Mandelbrot set. Given the isomorphism between complex numbers and vectors using the geometric product, we can use write

\begin{aligned}Z &= \mathbf{x} \hat{\mathbf{n}} \\ C &= \mathbf{c} \hat{\mathbf{n}},\end{aligned} \hspace{\stretch{1}}(2.2)

and reexpress the Julia iterator as

\begin{aligned}\mathbf{x} \rightarrow (\mathbf{x} \hat{\mathbf{n}})^p \hat{\mathbf{n}} + \mathbf{c}\end{aligned} \hspace{\stretch{1}}(2.4)

It’s not obvious that the RHS of this equation is a vector and not a multivector, especially when the vector \mathbf{x} lies in \mathbb{R}^{3} or higher dimensional space. To get a feel for this, let’s start by write this out in components for \hat{\mathbf{n}} = \mathbf{e}_1 and p=2. We obtain for the product term

\begin{aligned}(\mathbf{x} \hat{\mathbf{n}})^p \hat{\mathbf{n}} &= \mathbf{x} \hat{\mathbf{n}} \mathbf{x} \hat{\mathbf{n}} \hat{\mathbf{n}} \\ &= \mathbf{x} \hat{\mathbf{n}} \mathbf{x} \\ &= (x_1 \mathbf{e}_1 + x_2 \mathbf{e}_2  )\mathbf{e}_1(x_1 \mathbf{e}_1 + x_2 \mathbf{e}_2  ) \\ &= (x_1 + x_2 \mathbf{e}_2 \mathbf{e}_1 )(x_1 \mathbf{e}_1 + x_2 \mathbf{e}_2  ) \\ &= (x_1^2 - x_2^2 ) \mathbf{e}_1 + 2 x_1 x_2 \mathbf{e}_2\end{aligned}

Looking at the same square in coordinate representation for the \mathbb{R}^{n} case (using summation notation unless otherwise specified), we have

\begin{aligned}\mathbf{x} \hat{\mathbf{n}} \mathbf{x} &= x_k \mathbf{e}_k \mathbf{e}_1x_m \mathbf{e}_m  \\ &= \left(x_1 + \sum_{k>1} x_k \mathbf{e}_k \mathbf{e}_1\right)x_m \mathbf{e}_m  \\ &= x_1 x_m \mathbf{e}_m +\sum_{k>1} x_k x_m \mathbf{e}_k \mathbf{e}_1 \mathbf{e}_m \\ &= x_1 x_m \mathbf{e}_m +\sum_{k>1} x_k x_1 \mathbf{e}_k +\sum_{k>1,m>1} x_k x_m \mathbf{e}_k \mathbf{e}_1 \mathbf{e}_m \\ &= \left(x_1^2 -\sum_{k>1} x_k^2\right) \mathbf{e}_1+2 \sum_{k>1} x_1 x_k \mathbf{e}_k +\sum_{1 < k < m, 1 < m < k} x_k x_m \mathbf{e}_k \mathbf{e}_1 \mathbf{e}_m \\ \end{aligned}

This last term is zero since \mathbf{e}_k \mathbf{e}_1 \mathbf{e}_m = -\mathbf{e}_m \mathbf{e}_1 \mathbf{e}_k, and we are left with

\begin{aligned}\mathbf{x} \hat{\mathbf{n}} \mathbf{x} =\left(x_1^2 -\sum_{k>1} x_k^2\right) \mathbf{e}_1+2 \sum_{k>1} x_1 x_k \mathbf{e}_k,\end{aligned} \hspace{\stretch{1}}(2.5)

a vector, even for non-planar vectors. How about for an arbitrary orientation of the unit vector in \mathbb{R}^{n}? For that we get

\begin{aligned}\mathbf{x} \hat{\mathbf{n}} \mathbf{x} &=(\mathbf{x} \cdot \hat{\mathbf{n}} \hat{\mathbf{n}} + \mathbf{x} \wedge \hat{\mathbf{n}} \hat{\mathbf{n}} ) \hat{\mathbf{n}} \mathbf{x}  \\ &=(\mathbf{x} \cdot \hat{\mathbf{n}} + \mathbf{x} \wedge \hat{\mathbf{n}} ) (\mathbf{x} \cdot \hat{\mathbf{n}} \hat{\mathbf{n}} + \mathbf{x} \wedge \hat{\mathbf{n}} \hat{\mathbf{n}} )   \\ &=((\mathbf{x} \cdot \hat{\mathbf{n}})^2 + (\mathbf{x} \wedge \hat{\mathbf{n}})^2) \hat{\mathbf{n}}+ 2 (\mathbf{x} \cdot \hat{\mathbf{n}}) (\mathbf{x} \wedge \hat{\mathbf{n}}) \hat{\mathbf{n}}\end{aligned}

We can read 2.5 off of this result by inspection for the \hat{\mathbf{n}} = \mathbf{e}_1 case.

It is now straightforward to show that the product (\mathbf{x} \hat{\mathbf{n}})^p \hat{\mathbf{n}} is a vector for integer p \ge 2. We’ve covered the p=2 case, justifying an assumption that this product has the following form

\begin{aligned}(\mathbf{x} \hat{\mathbf{n}})^{p-1} \hat{\mathbf{n}} = a \hat{\mathbf{n}} + b (\mathbf{x} \wedge \hat{\mathbf{n}}) \hat{\mathbf{n}},\end{aligned} \hspace{\stretch{1}}(2.6)

for scalars a and b. The induction test becomes

\begin{aligned}(\mathbf{x} \hat{\mathbf{n}})^{p} \hat{\mathbf{n}} &= (\mathbf{x} \hat{\mathbf{n}})^{p-1} (\mathbf{x} \hat{\mathbf{n}}) \hat{\mathbf{n}} \\ &= (\mathbf{x} \hat{\mathbf{n}})^{p-1} \mathbf{x} \\ &= (a + b (\mathbf{x} \wedge \hat{\mathbf{n}}) ) ((\mathbf{x} \cdot \hat{\mathbf{n}} )\hat{\mathbf{n}} + (\mathbf{x} \wedge \hat{\mathbf{n}}) \hat{\mathbf{n}}) \\ &= ( a(\mathbf{x} \cdot \hat{\mathbf{n}} )^2 - b (\mathbf{x} \wedge \hat{\mathbf{n}})^2 ) \hat{\mathbf{n}}+ ( a + b(\mathbf{x} \cdot \hat{\mathbf{n}} ) ) (\mathbf{x} \wedge \hat{\mathbf{n}}) \hat{\mathbf{n}}.\end{aligned}

Again we have a vector split nicely into projective and rejective components, so for any integer power of p our iterator 2.4 employing the geometric product is a mapping from vectors to vectors.

There is a striking image in the text of such a Julia set for such a 3D iterator, and an exersize left for the adventurous reader to attempt to code that based on the 2D p=2 sample code they provide.

References

[1] L. Dorst, D. Fontijne, and S. Mann. Geometric Algebra for Computer Science. Morgan Kaufmann, San Francisco, 2007.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

 
%d bloggers like this: