# Fractals/Mathematics/Vector field

Vector field[1]

Here mainly numerical methods for time independent 2D vector fields are described.

# Dictionary

• vector function is a function that gives vector as an output
• field : space (plane, sphere, ... )
• field line is a line that is everywhere tangent to a given vector field
• scalar/ vector / tensor:
• Scalars are real numbers used in linear algebra. Scalar is a tensor of zero order
• Vector is a tensor of first order. Vector is an extension of scalar
• tensor is an extension of vector

## Vector

Forms of 2D vector:[2]

• [z1] ( only one complex number when first point is known , for example z0 is origin
• [z0, z1] = two complex numbers
• 4 scalars ( real numbers)
• [x, y, dx , dy]
• [x0, y0, x1, y1]
• [x, y, angle, magnitude]
• 2 scalars : [x1, y1] for second complex number when first point is known , for example z0 is origin

The numerical gradient of a function:

• "is a way to estimate the values of the partial derivatives in each dimension using the known values of the function at certain points."[3]

Gradient function G of function f at point (x0,y0)

  G(f,x0,y0) = (x1,y1)


Input

• function f
• point (x0,y0) where gradient is computed

Output:

• vector from (x0,yo) to (x1,y1) = gradient

See:

Computation:[4]

• "The gradient is calculated as: (f(x + h) - f(x - h)) / (2*h) where h is a small number, typically 1e-5 f(x) will be called for each input elements with +h and -h pertubation. For gradient checking, recommend using float64 types to assure numerical precision."[5]
• in matlab[6][7]
• in R[8]
• python[9]

## equation

ODE means Ordinary Differential Equation, where "ordinary" means with derivative respect to only one variable (like ${\displaystyle {\frac {d}{dt}}}$ ), as opposed to an equation with partial derivatives (like ${\displaystyle {\frac {\partial }{\partial x}}}$ , ${\displaystyle {\frac {\partial }{\partial y}}}$ , ...) called PDE. (matteo.basei)


# Field types

Criteria for classification

Force types A gravitational force fields

• the field lines are the solution of ${\displaystyle g'=F}$
• the trajectory of a test mass is the solution of ${\displaystyle mg''=F}$

where

• g is the standard gravity
• m is a mass
• F is a force field

An electric field

• The field lines are the paths that a point positive charge would follow as it is forced to move within the field. Field lines due to stationary charges have several important properties, including always originating from positive charges and terminating at negative charges, they enter all good conductors at right angles, and they never cross or close in on themselves.  The field lines are a representative concept; the field actually permeates all the intervening space between the lines. More or fewer lines may be drawn depending on the precision to which it is desired to represent the field. The study of electric fields created by stationary charges is called electrostatics.

# Algorithm

• plane (parameter plane or dynamic plane)
• scalar function
• vector function
• create scalar field using scalar function ( potential)
• create vector field from scalar field using vector function ( gradient of the potential)
• compute:
• filed lines ( stream lines )
• contour lines ( equipotential lines )
• map whole field using
• Line Integral Convolution (LIC)

//https://editor.p5js.org/ndeji69/sketches/EA17R4HHa
// p5 js Tutorial】Swirling Pattern using Gradient for Generative Art by Nekodigi

function curl(x, y){
var EPSILON = 0.001;//sampling interval
//Find rate of change in X direction
var n1 = noise(x + EPSILON, y);
var n2 = noise(x - EPSILON, y);
//Average to find approximate derivative
var cx = (n1 - n2)/(2 * EPSILON);

//Find rate of change in Y direction
n1 = noise(x, y + EPSILON);
n2 = noise(x, y - EPSILON);

//Average to find approximate derivative
var cy = (n1 - n2)/(2 * EPSILON);

//return new createVector(cx, cy);//gradient toward higher position
return new createVector(cy, -cx);//rotate 90deg
}

function draw() {
tint(255, 4);
image(noiseImg, 0, 0);//fill with transparent noise image
//fill(0, 4);
//rect(0, 0, width, height);

strokeWeight(4);//particle size
stroke(255);

for(var i=0; i<particles.length; i++){
var p = particles[i];//pick a particle
point(p.pos.x, p.pos.y);
}
}


## separatrix

The separatrix is clearly visible by numerically solving for trajectories backwards in time. Since when solving for the trajectories forwards in time, trajectories diverge from the separatrix, when solving backwards in time, trajectories converge to the separatrix.

## Field line computing

Problem Statement:

• Field line tracing ( not curve sketching[11]}
• drawing contour maps ( in computer graphic) = Numerical continuation ( in math)
• compute an integral curve from a seed point through a vector field without any analysis of its structure on the uniform grid ( raster scan or pixels)

Methods ( solvers) available for the field-lines[12][13]

  None of these 4 methods generate an exact answer, but they are (from left to right) increasingly more accurate. They also take (from left to right) more and more time to finish as they require more samples for each iteration.
You won't be able to create reliably closed curves using iterative sampling methods as small errors at any step may be amplified in successive steps. There is also no guarantee that the field-line ends up in the exact coordinate where it started.
The Grasshopper metaball solver on the other hand uses a marching squares algorithm which is capable of finding closed loops because it is a grid-cell approach and sampling inaccuracy in one area doesn't carry over to another.
However the solving of iso-curves is a very different process    from the solving of particle trajectories through fields. ...
Typically field lines shoot to infinity rather than form closed loops. That is one reason why I chose the RK methods here, because marching-cubes is very bad at dealing with things that tend to infinity.[16]


#### Construction

Given a vector field ${\displaystyle \mathbf {F} (\mathbf {x} )}$  and a starting point ${\displaystyle \mathbf {x} _{\text{0}}}$  a field line can be constructed iteratively by finding the field vector at that point ${\displaystyle \mathbf {F} (\mathbf {x} _{\text{0}})}$ . The unit tangent vector at that point is: ${\displaystyle \mathbf {F} (\mathbf {x} _{\text{0}})/|\mathbf {F} (\mathbf {x} _{\text{0}})|}$ . By moving a short distance ${\displaystyle ds}$  along the field direction a new point on the line can be found

${\displaystyle \mathbf {x} _{\text{1}}=\mathbf {x} _{\text{0}}+{\mathbf {F} (\mathbf {x} _{\text{0}}) \over |\mathbf {F} (\mathbf {x} _{\text{0}})|}ds}$

Then the field at that point ${\displaystyle \mathbf {F} (\mathbf {x} _{\text{1}})}$  is found and moving a further distance ${\displaystyle ds}$  in that direction the next point of the field line is found

${\displaystyle \mathbf {x} _{\text{2}}=\mathbf {x} _{\text{1}}+{\mathbf {F} (\mathbf {x} _{\text{1}}) \over |\mathbf {F} (\mathbf {x} _{\text{1}})|}ds}$

By repeating this and connecting the points,the field line can be extended as far as desired. This is only an approximation to the actual field line, since each straight segment isn't actually tangent to the field along its length, just at its starting point. But by using a small enough value for ${\displaystyle ds}$ , taking a greater number of shorter steps, the field line can be approximated as closely as desired. The field line can be extended in the opposite direction from ${\displaystyle \mathbf {x} _{\text{0}}}$  by taking each step in the opposite direction by using a negative step ${\displaystyle -ds}$ .

### rk4 numerical integration method

Fourth-order Runge-Kutta (RK4) in case of 2D time independent vector field

${\displaystyle F}$  is a vector function that for each point p

p = (x, y)

in a domain assigns a vector v

${\displaystyle v=F(p)=F(x,y)=(F_{1}(x,y),F_{2}(x,y))}$

where each of the functions ${\displaystyle F_{i}}$  is a scalar function:

${\displaystyle F_{i}:\mathbb {R} ^{2}\to \mathbb {R} }$

A field line is a line that is everywhere tangent to a given vector field.

Let r(s) be a field line given by a system of ordinary differential equations, which written on vector form is:

${\displaystyle {\frac {dr}{ds}}=F(r(s))}$

where:

• s representing the arc length along the field line, like for example continous iteration count
• ${\displaystyle r(0)=r_{0}}$  is a seed point

#### 2 variables

Given a seed point ${\displaystyle r_{0}}$  on the field line, the update rule ( RK4) to find the next point ${\displaystyle r_{i}}$ along the field line is[17]

${\displaystyle r_{i+1}=r_{i}+{\frac {h(k_{1}+2k_{2}+2k_{3}+k_{4})}{6}}}$

where:

• h is the step size along field line = ds
• k are the intermediate vectors:

${\displaystyle {\begin{array}{lcl}k_{1}=F(r_{i})\\k_{2}=F(r_{i}+{\frac {h}{2}}k_{1})\\k_{3}=F(r_{i}+{\frac {h}{2}}k_{2})\\k_{4}=F(r_{i}+hk_{3})\end{array}}}$

#### only x

Here

${\displaystyle y'={\frac {dy}{dx}}=f(x)}$


Given a seed point ${\displaystyle r_{0}}$  on the field line, the update rule ( RK4) to find the next point ${\displaystyle r_{i}}$ along the field line is[18]

${\displaystyle x_{i+1}=x_{i}+dx=x_{i}+h}$
${\displaystyle y_{i+1}=y_{i}+{\frac {h(k_{1}+2k_{2}+2k_{3}+k_{4})}{6}}}$


where:

• h is the step size along field line = dx
• k are the intermediate vectors:

${\displaystyle {\begin{array}{lcl}k_{1}=F(x_{i})\\k_{2}=F(x_{i}+{\frac {h}{2}}k_{1})\\k_{3}=F(x_{i}+{\frac {h}{2}}k_{2})\\k_{4}=F(x_{i}+hk_{3})\end{array}}}$

Examples:

• ${\displaystyle {\frac {dy}{dx}}=x}$
• ${\displaystyle y^{'}+x^{2}+x=0}$ [19]
• ${\displaystyle {\frac {dy}{dx}}=x^{2}-x-2}$

## Visualisation of vector field

Plot types (Visualization Techniques for Flow Data) : [20]

• Glyphs = Icons or signs for visualizing vector fields
• simplest glyph = Line segment (hedgehog plots)
• arrow plot = quiver plot = Hedgehogs (global arrow plots)
• Characteristic Lines [21]
• streamlines = curve everywhere tangential to the instantaneous vector (velocity) field (time independent vector field). For time independent vector field streaklines = Path lines = streak lines [22]
• texture (line integral convolution = LIC)[23]
• Topological skeleton [24]
• fixed point extraction ( Jacobian)
  "path lines, streak lines, and stream lines are identical for stationary flows" Leif Kobbelt[25]


### quiver plot

Definition

• "A quiver plot displays velocity vectors as arrows with components (u,v) at the points (x,y)"[26]

# Example fields

## fBM

fBM stands for Fractional Brownian Motion

## SDF

SDF = Signed Distance Function[27]

• by dimension ( 1D, 2D, 3D, ...)
• by color
• by distance function ( Euclid distance,
• algorithm [28]
• the efficient fast marching method,
• ray marching [29]
• like Marching Parabolas, a linear-time CPU-amenable algorithm.
• Min Erosion, a simple-to-implement GPU-amenable algorithm
• fast sweeping method
• the more general level-set method.
• visualisation ( gray gradient, LSM,
• simple predefined figures or arbitrary shape

It is not

### circle

// The MIT License
// Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
// GSLS

// Signed distance to a disk

// List of some other 2D distances: https://www.shadertoy.com/playlist/MXdSRf
//
// and iquilezles.org/articles/distfunctions2d

float sdCircle( in vec2 p, in float r )
{
return length(p)-r;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
vec2 p = (2.0*fragCoord-iResolution.xy)/iResolution.y;
vec2 m = (2.0*iMouse.xy-iResolution.xy)/iResolution.y;

float d = sdCircle(p,0.5);

// coloring
vec3 col = (d>0.0) ? vec3(0.9,0.6,0.3) : vec3(0.65,0.85,1.0); // exterior / interior
col *= 1.0 - exp(-6.0*abs(d)); // adding a black outline ( gray gradient) to the circle
col *= 0.8 + 0.2*cos(150.0*d); //  adding waves
col = mix( col, vec3(1.0), 1.0-smoothstep(0.0,0.01,abs(d)) ); // note: adding white border to the circle

fragColor = vec4(col,1.0);
}


### hg_sdf

hg_sdf: A glsl library for building signed distance functions

### blender

• b3dsdf = A toolkit of 2D/3D distance functions, sdf/vector ops and various utility shader nodegroups (159+) for Blender 2.83+.

# Examples

## Videos

### rboyce1000

• Bifurcation of Quartic Polynomial Vector Fields by @rboyce1000
The coloured curves are the separatrices (i.e. real flow lines reaching infinity) of the complex ODE dz/dt = p(z) = z^4 + O(z^2), where the four roots of p(z) are pictured as the black dots: one fixed at the origin, and the remaining three forming the vertices of an equilateral triangle centered at the origin and rotating.
Bifurcation occurs at certain critical angles of the rotation, where separatrices instantaneously merge to form homoclinic orbits. Following bifurcation, the so-called 'sectorial pairing' is permuted.
There are a total of 5 possible sectorial pairings for the quartic polynomial vector fields (enumerated by the 3rd Catalan number).
Three out of the five possibilities can be seen in 2/3 video, while the remaining two can be seen in part 1/3
In 3/3 example, at a bifurcation we have that either: only two of the four roots are centers (the other two remaining attached to separatrices), or NONE of the roots are centers (a phenomenon which does not occur for the quadratic or cubic polynomial vector fields).
This video is inspired by the work of A. Douady and P. Sentenac. ( rboyce1000)


Algorithm

• create polynomial with desired properities
• f(z) = z*g(z) with root at origin
• g(z) is a 3-rd root of unity = ${\displaystyle z^{3}-1}$

f(z) = z(z^3 - 1)

One can check it with Maxima CAS

 z:x+y*%i;
(%o1)                              %i y + x
(%i2) p:z*(z^3-1);
3
(%o2)                    (%i y + x) ((%i y + x)  - 1)
(%i3) display2d:false;

(%o3) false
(%i4) r:realpart(p);

(%o4) x*((-3*x*y^2)+x^3-1)-y*(3*x^2*y-y^3)
(%i5) m:imagpart(p);

(%o5) x*(3*x^2*y-y^3)+y*((-3*x*y^2)+x^3-1)
(%i6) plotdf([r,m],[x,y]);

(%o6) "/tmp/maxout28945.xmaxima"

(%i8) s:solve([p],[x,y]);

(%o8) [[x = %r1,y = (2*%i*%r1+%i+sqrt(3))/2],
[x = %r2,y = (2*%i*%r2+%i-sqrt(3))/2],[x = %r3,y = %i*%r3],
[x = %r4,y = %i*%r4-%i]]



to rotate it around origin let's change 1 with :${\displaystyle \quad e^{2\pi it}}$  ( multiplier of the fixed point) where t is a proper fraction in turns

${\displaystyle f_{t}(z)=z(z^{3}-e^{2\pi it})}$

${\displaystyle f_{t}(z)=z(z^{3}-1.5^{3}e^{2\pi it})}$