A question about a fractal like iteratively defined function

436 Views Asked by At

I am trying to figure out what the following function $f:\Bbb{R}^3-\{\mathbf{0}\}\to\Bbb{R}$ defined below (in pseudocode) does:


function $f(\mathbf{v}\in \Bbb{R}^3-\{\mathbf{0}\})$
{
    $a=b=0;$      $(a,b\in \Bbb{R})$
    for ${(i=0; i<I; i=i+1)}$
    {
            $\mathbf{v}=\frac{|\mathbf{v}|}{||\mathbf{v}||^2}-\mathbf{A};$
            $a=a+|||\mathbf{v}||-b|;$
            $b=||\mathbf{v}||;$
    }
    return $ a;$
}

Where $I$ is some positive integer that determines the number of iterations the loop will run, and $\mathbf{A}$ is some vector in the set $\{(t,t,t)\in\Bbb{R}^3|t\in[-1,1]\}$.

(The notation $|\mathbf{v}|$ denotes component wise absolute value of the vector elements)


I do not know how this function looks like and what it does and I am only interested in vectors $\mathbf{v}\in[-1,1]^3-\{\mathbf{0}\}$.

Also what is $lim_{\mathbf{v}\to \mathbf{0}}f(v)$?


Here are some picks that were generated by my shader (maybe they can help):

When $I = 1$ and $\mathbf{A}$ is any vector in the set $\{(t,t,t)\in\Bbb{R}^3|t\in[-1,1]\}$ the images will be something like this enter image description here

As we increase the value of $I$ things get interesting and we get imaged like those:

When $I = 10$ and $\mathbf{A}=\mathbf{(0.5, 0.5, 0.5)}$ enter image description here

When $I = 10$ and $\mathbf{A}=\mathbf{(0.63, 0.63, 0.63)}$ enter image description here

When $I = 10$ and $\mathbf{A}=\mathbf{(0.74, 0.74, 0.74)}$ enter image description here

When $I = 10$ and $\mathbf{A}=\mathbf{(1, 1, 1)}$ enter image description here

Since the center of every image is solid white I can conjecture that $lim_{\mathbf{v}\to\mathbf{0}}f(\mathbf{v})=\infty$


the main part of my GLSL fragment shader that produces the images and implements the function $f$ defind above is

vec3 color = vec3(0.);
vec3 v = vec3(2 * (TexCoord.xy - .5), 0.); // TexCoord.xy in [0,1]^2
float a = 0.;
float b = 0.;

// The loop that I do not understand
for (uint i = 0U; i < I; i++)
{
    v = abs(v) / dot(v,v) - A;
    a += abs(length(v) - b);
    b = length(v);
}

a = a * a * a; // Adding some contrast

float brightness = 0.015;
color = vec3(.1, .01, .0001) * a * brightness;

FragColor = vec4(color, 1.);

Any help will be appreciated.

(I've encountered this function while trying to achieve some interesting fractal like visual effects while doing graphics programming)

EDIT

In all pictures, the parameter $Z$ denotes the $Z$ coordinate of the vector $\mathbf{v}$

$I=90$, $A=0$ and $Z=0$ enter image description here

$I=90$, $A=0.5$ and $Z=0$ enter image description here $I=90$, $A=1$ and $Z=0$ enter image description here $I=100$, $A=0$ and $Z=0$ enter image description here $I=100$, $A=0.5$ and $Z=0$ enter image description here $I=100$, $A=1$ and $Z=0$ enter image description here

$I=100$, $A= 0.042373 $ and $Z=0$ enter image description here $I=100$, $A= 0.906780 $ and $Z=0$ enter image description here $I=100$, $A= 0.932203 $ and $Z=0$ enter image description here

$I=100$, $A= 0.046610 $ and $Z= 0.894068 $ (cool dotted effect) enter image description here

$I=16$, $A= 0.194915 $ and $Z=0$ enter image description here $I=29$, $A= 0.194915 $ and $Z=0$ enter image description here

$I=500$, $A=0$ and $Z=0$ enter image description here $I=500$, $A=0.5$ and $Z=0$ (nearly all white) enter image description here $I=500$, $A=1$ and $Z=0$ enter image description here

$I=1000$, $A=0$ and $Z=0$ enter image description here

$I=1000$, $A=0.5$ and $Z=0$ (nearly all white) enter image description here

$I=1000$, $A=1$ and $Z=0$ enter image description here


Here are some images that view the first iteration of the loop on the xy plane with $A=(0.5,0.5,0)$:

(I will try to generate images for greater iterations)

Above the xy plane view enter image description here

Below the xy plane view enter image description here

2

There are 2 best solutions below

11
On BEST ANSWER

I feel the need to take this to a more mathematical realm. So here's a commentary to long for comments.

Set-Up

We have a system of recurrence relations defined by,

$$(1) \quad \vec v_{n+1}={{|\vec v_n|} \over {||\vec v_n||^2}}-\vec A$$ $$(2) \quad a_{n+1}=a_n+| \ ||\vec v_n||-b_n|$$ $$(3) \quad b_{n+1}=||\vec v_n||$$

First let's reduce this to a system of two relations. $(3)$ is redundant, we can simply use $||v_{n+1}||$. We now have,

$$(1) \quad \vec v_{n+1}={{|\vec v_n|} \over {||\vec v_n||^2}}-\vec A$$ $$(4) \quad a_{n+1}=a_n+| \ ||\vec v_n||-||\vec v_{n-1}|| \ |$$

$(1)$ happens to be independent of $(4)$ so I'd suggest that we understand $(1)$ before trying to figure out $(4)$, which is just the picture generator.

Interpretations

It may be hard to interpret what $(1)$ and $(4)$ are actually doing. Here's a briefer. $(4)$ is easy. It measures how volatile the path taken from $v_0$ to $v_n$ was. A higher value means the path was longer and more volatile. A lower value means it either didn't have a value explosion, or that it converged to some value. The interpretation of $(1)$ is more subtle. I still don't know exactly what it does, but it seems to normalize the vector to a unit vector, then it does that again, then it subtracts a constant and we start over. Because the vector is normalized, most of the interesting behavior happens in the unit circle around the origin.

This system is Chaotic

We can start by constructing a toy model. We'll do this by collapsing $(1)$ to only a recurrence in $1$-d.

$$(5) \quad v_{n+1}={1 \over {|v_n|}}-c$$

Which should retain some of the same characteristics. You seem interested in the case where $v \rightarrow 0$. Using the toy model $(5)$ we see that,

$$v_0=0$$ $$v_1=\infty$$ $$v_2=-c$$ $$v_3={1 \over c}-c$$ $$v_4={1 \over {|{1 \over c}-c|}}-c$$ And so on...

The general behavior of the iterations is very similar to continued fractions. The only difference is the absolute value sign. This changes the behavior dramatically. For instance here's the behavior of $v_0=0.7$ and $c=0.5$ under iteration.

Figure 1

The behavior is chaotic. So basically, this toy model shows you why your getting incredibly complex behavior in the $3$-d model. It's the inclusion of the absolute value that's doing it! Compare this to the version without it here.

Figure 2

Note that it has an explicit formula for the evolution.

Singularity Analysis

Let's move to the $2$-d case, equivalent to setting all of the $z$ components to $0$. I'll skip the algebra, but here's the evolution of the vector system.

$$(6) \quad x_{n+1}={{|x_n|} \over {x_n^2+y_n^2}}-v_x$$ $$(7) \quad y_{n+1}={{|y_n|} \over {x_n^2+y_n^2}}-v_y$$

Once again, we'll ignore the detailed analysis of the path indicator function. (I hope to show why that is later). By inspection we can see that the function will become infinite whenever a vector passes through infinity. By inspection we can see that this happens if the initial components satisfy,

$$x_0^2+y_0^2=0$$

So it really wouldn't be difficult to prove that for the $3$-d case with $z_0$, the same equation must be satisfied. So this answers your question. $\lim_{\vec v \to 0} a_n=\infty$. So the path indicator function becomes infinite. That's why the origin is white in your pictures. The path indicator function also becomes infinite at iteration $n=2$ if,

$$v_x^2 \cdot (x_0^2+y_0^2)+v_y^2 \cdot (x_0^2+y_0^2)+1=2 \cdot v_x \cdot |x_0|+2 \cdot v_y \cdot |y_0|$$

If one could set up a suitable renormalization scheme, it might be possible to find all the zeros within some suitable error tolerance.

Show and Tell

Here's my attempt at providing pictures. If you wondering about the quality, it was made on a DS. We have a logarithmic path indicator function. It's identical to yours, except I take a log after we're done calculating it. The range of the graphic is the square from $-1 \le x \le 1$ and $-1 \le y \le 1$. Use the origin to base off what the colors mean, as far the path indicator function is concerned.

Below is what results when the constant vector is $v_x=v_y=0.1$.

Figure 3

Using $v_x=v_y=0.5$, we get qualitatively different behavior.

Figure 4

Let's add some asymmetry with $v_x=0.25, \ v_y=0$. The line at $y=0$ corresponds to the $1$-D case we discussed above.

Figure 5

Fixed and Singularity Points

So the pictures I posted have very intricate and complex patterns. However, there are two pieces of these pictures that can be readily explained. These are the fixed and singularity points. These values are pretty much what they sound like. Fixed points are points where the path indicator function eventually becomes fixed in value. Singularity points are points where the path indicator function becomes infinite after some iteration. Finding these points becomes very difficult so I'll be only giving an intro to this topic. In the one dimensional case, the recurrence $(5)$ has fixed points at,

$$(8) \pm \cfrac{-c+\sqrt{c^2+4}}{2}$$

If your being keen, you'll notice that the negative solution is false. However, if your very keen you'll notice that the negative solution will evolve to the real fixed point in one iteration. Basically a fixed point of $-2$ is false, but after one iteration, it will become $2$ which is a fixed point.

To find the singularity points, we have to find values such that the iteration becomes infinite. In our case, we can do this by finding values such that the iteration is $0$. If we find these values, we know that the very next value will be infinite. Looking at $(5)$ again, we find that for an initial point $v_0$ to evolve to an infinite point after $1$ iteration, $v_0$ must equal $0$. For, the case of two iterations, we solve,

$$0={1 \over {|v_0|}}-c$$

Which implies that $v_0=\pm {1 \over c}$. We can then move on to finding the values of $v_0$ that give infinite values after three iterations and so forth.

Let's use the third picture for an example. We have $v_x=0.25, \ v_y=0$ so we can our above analysis quite nicely. Using $(8)$ we find the fixed points to be

$$\pm {{\sqrt{65}-1} \over 8}=\pm 0.8828...$$

That corresponds to the dark areas on the $y=0$ line in the third picture. Similarly we can find the first three singularity points to be, $0$, $\pm 4$, and $\pm 4/17=0.235...$. By the way, formulas exist to find the next $6$ singularity points.

There are similar, albeit more complicated, formulas for the two and three dimensional cases. For instance, I provide a formula that can be utilized for finding singularity points in another section.

Potpourri

The final interesting thing to note is that the iteration of $v$, in the two dimensional case, can be turned into a complex iteration equivalent to,

$$z_{n+1}=\cfrac{\operatorname{Re}(z_n)+\operatorname{Im}(z_n)}{|z_n|^2}$$

Where we have the real and imaginary part operators involved.

I may or may not add something else...

6
On

So first out some theoretical pondering

When comes to theoretical investigation, since we have considered 2d vectors so far, I will make the simplification to move to $z = r+wi \in \mathbb{C}^1$ where , the operation $|{\bf v}|$ takes $z = r + wi$ to $|r| + |w|i$, i.e.

  1. Moves complex number to the first quadrant by a series of reflections in the real and imaginary axes
  2. Normalizes it's length.
  3. Subtracts a prescribed vector which is also always in the first quadrant.

Also, we can modify $\frac{|v|}{\|v\|^2} \rightarrow \frac{|v|+\epsilon}{\|v\|^2+\epsilon}$, to avoid possible issue of having the denominator too close to 0 to cause numerical instability or division by 0.

The rest of the operations do not involve the information of which quadrant we were in any way. Therefore we should be able to perform the inverse operations on $\bf A$ instead of forward operations on $\bf v$ each iteration. Then what $a$ does is basically to sum the fluctuations of the length of $\bf v$ over the iterations. The absolute value of this number is the distance in the complex plane to the point $\bf A$. So we measure how fast we swing around $\bf A$ in some sense. The dramatic points light up and the ones who remain at the same distance are dark.


Here is my octave implementation. Code can be found at the bottom. It sure does some strange fractal-like thing anyway. But before the code, here's some animations.

I choose the colors as gray scale ( each channel the same ): $\|a\|^{\frac{1}{4}}$ and the plots are on the square $x,y \in \left[-\frac{1}{2},\frac{1}{2}\right]^2$

enter image description here

Above is for ${\bf A} = \frac{1}{3}[1,1,1]^T$, 25 iterations.

enter image description here

Above is for ${\bf A} = \frac{1}{2}[1,1,1]^T$, 25 iterations.

enter image description here

Above is for ${\bf A} = [1,1,1]^T$, 25 iterations.

And here's the Octave code:

A = [1;1;1]/8;

lRes = 512;

%[xs,ys,zs] = >>ndgrid(linspace(-0.5,0.5,lRes),linspace(-0.5,0.5,lRes),linspace(eps,eps,1));

[xs,ys,zs] = ndgrid(linspace(-1,1,lRes),linspace(-1,1,lRes),linspace(eps,eps,1));

v = [vec(xs)';vec(ys)';vec(zs)']';

lIters = 30;

a=0; b=0;

for i_=1:lIters

v = abs(v)./(sum(v.^2,2)) - A';

a = a + abs(sum(v.^2,2).^0.5-b);

b = sum(v.^2,2).^0.5;

lImage = reshape(a(:,:,1),[lRes,lRes]).^0.25;

imagesc(lImage); colormap gray; drawnow;

end