Benutzer:Georg-Johann/Mathematik

A Julia Set (in white) visualized using CCA

 

Visualising Julia sets

Julia sets and Fatou sets can be nice. Their computergraphical generation can be hard.

In the remainder of this page, we work out a method which can be used to visualise the Julia set of a complex function ƒ. The advantage is that you don't need to know attractors of the iteration

$z\mapsto f(z)$

The generated images will be smooth in the Fatou part.

Overview: Imaging methods

Note that there exist other approaches to color complex dynamical systems like

Inverse Iteration Method (IIM)

Compute the preimages of ƒ, i.e. compute the reverse orbit. Because the stability of the fixed points turns from attractive to repelling and vice versa, one just choses a complex number and looks where it goes under the invers iteration. The trouble is that the results will not uniformly distributet in $J$ and that you have to compute the inverse of ƒ.

Coloring by speed of attraction (CSA)

Taking a value from the lattice of points to color, perform iterations until the iterative is close to an attractor. Color the point according to the number of iterations needed to bring it close enough to the attractor.

This method is commonly used to visualize Julia sets of polynomials and Julia sets that are attached to Newton's famous method for finding the zeros of a function. Polynomials or degree > 1 always have infinity as a super-attractive fixed point. The rational function that occurs in Newton's method has always the function's roots as attractive fixed points. However, in both cases there may be other attractors, which – moreover – need not to consist of just one point.

Escape Time Algorithm (ETA)

If ∞ is an attractor, i.e. a fixed point of the process, then color the point according to the number of iterations – the time – it takes until one sees the point escapes towards ∞. If the point does not escape during the maximum number of iterations, the point is colored as belonging to the Julia set or to the basin of some other attractor. This method works for polynomials. The most prominent Julia sets are the ones for z→z²+c where c is an element of the Mandelbrot set or not far away from the mandelbrot set. If you see a picture of such a Julia set, it is likely that ETA had been used to get the picture.

Cauchy Convergence Algorithm (CCA)

In the remainder of this page, I will present a different approach whose idea as basically the same es that of Escape Time Algorithm. However, no basin of attraction must be known in advance and the different basins of attraction can be separated and be colored differently. The approch uses the notion of Cauchys convergence. Instead of observing the orbit of a point, this method observes how the distance of two nearby points z and z+ε evolves as these two values are iterated. If the difference tends to 0, then the point heads for an attractor. If the difference does not approach 0, then the point is close to (or a part of) the Julia set.

The Metric

Let $\varrho$ be a canonical projection of the compactified complex plane onto the Riemann sphere S₂:

$\varrho :{\overline {\mathbb {C} }}=\mathbb {C} \cup \{\infty \}\rightarrow S_{2}$

This gives us a metric d: As distance between two points in the complex plane we take their distance on the sphere, i.e. the length of the orthodrome. This means that the metric is bounded by π and even the distance to ∞ (which is now the north pole) is finite.

In order to compute the distance between two points z and w we rotate the sphere S₂ in such a way that

1. w maps to 0
2. z maps to the positive real axis

After these transformations the distance can be computed quite easily. The rotation can be accomplished by a squence of isometric Möbius transformations. All an all, we get

$t_{w}:z\mapsto {\begin{cases}|z|\;,&{\text{ if }}w=0\\{\textstyle {\frac {1}{|z|}}}\;,&{\text{ if }}w=\infty \\{\textstyle {\frac {1}{|w|}}}\;,&{\text{ if }}z=\infty \\&\\{\textstyle {\frac {1}{|w|}}}\cdot \left|{\frac {z{\overline {w}}-|w|^{2}}{z{\overline {w}}+1}}\right|\;,&{\text{ else}}\end{cases}}$

which is left as an exercise to the reader. The bar denotes complex conjugation. The metric is then

$d(z,w)=\,\!2\arctan(t_{w}(z))$

The nice feature that is introduced by d is that sequences that formerly diverged against infinity now converge towards infinity, i.e. towards the north pole of S₂.

Stability of Orbits

Recall the definition of the Julia set for a contraction mapping ƒ. The definition implies some facts on the stability of the iteration

$z\mapsto f(z)\mapsto f^{2}(z)\mapsto f^{3}(z)\mapsto \cdots$

where ƒⁿ denotes the n-th iterative of ƒ:

$f^{\,n}=f\circ f\circ \cdots \circ f$

The set

$\{f^{\,n}(z)\}_{n\in \mathbb {N} _{0}}$

is called Orbit of z (under ƒ).

The orbits of two points z and w behave similar − in some sense − if z and w lie close together and belong to the Fatou set F_ƒ which is the complement of the Julia set J_ƒ of ƒ. If z is an element of J_ƒ then z and w will behave quite different, even if w itself is an element of J_ƒ.

To get a notion of the stability of an orbit, we set

$\Sigma _{n}(z)=\Sigma _{n}(z,\varepsilon )=\sum _{k\leqslant n}d\,(f^{k}(z),f^{k}(z+\varepsilon ))$

for a small ε and with the metric d from above. This means that we take two points which are close together, and then we summarize their distances as ƒ makes them jump around on the Riemann sphere. Note that for any fixed ε the sum can diverge for large n even if z is a Fatou point.

However, we can use Σ_n to measure how close a point is to J_ƒ: the larger the sum is, the more instable is the iteration and the closer is the point to the Julia set of ƒ.

To destinguish points of (or close to) the Julia set from points in the Fatou set, we need a creterion. To get it, we compute all the Σ-values for the points that we want to color, i.e. for all points in a lattice Γ. After computing these values, we do a little bit of statistics to get the expected value E and the standard deviation σ for the set of Σ-values: Let

${\begin{aligned}E_{n}&=E_{n}(\Gamma ,\varepsilon )=E\left\{\log(\delta +\Sigma _{n}(z,\varepsilon ))\right\}_{\!z\,\in \,\Gamma }\\\sigma _{n}&={\sqrt {D_{n}}}\\\end{aligned}}$

Remind that

${\begin{aligned}EX&={\textstyle {\frac {1}{|X|}}}\sum \limits _{x\in X}x\\D\!X&=EX^{2}-(EX)^{2}\\\end{aligned}}$

It turns out that the values Σ are widely spread over several scales. Therefore, we do not use Σ directly. Instead, we use the logarithm of Σ. The value δ is just a small constant to avoid the logarithm's input to be zero.

Coloring

Now, we have all we need to color a point:

1. choose
   1. a lattice $\Gamma$ of points $z$ to color
   2. the number $n$ of iterations to perform
   3. a small $\varepsilon$
2. for all points, compute $\Sigma _{n}(z,\varepsilon )$
3. compute $E_{n}$ and $\sigma _{n}$
4. compute

$J(z)=\,\ln(\delta +\Sigma _{n}(z))-K$

for some constant $K$. Because $K$ will be used to separate points that belong to the Julia set ($J(z)>0$) from points that to not ($J<0$), reasonable values for $K$ are greater than $E_{n}$. Try settings like $K=E_{n}+2\sigma _{n}$ or $E_{n}+\sigma _{n}$ or the like. If $J(z)>0$ then we color $z$ as belonging to the Julia set. If $J(z)<0$ we can use that value to shade the Fatou set. If we know some attractor, we can check if ƒ^$n$$(z)$ is close to it and use that information, too.

To map values to the valid ranges for saturation and brightness we use the function $h$ from section helper function h.

Modifications

The computation of $E_{n}$ takes a lot of time. The visualisation process needs two passes:

1. compute $E_{n}$ from all $\Sigma _{n}(z)$
2. color the points using $E_{n}$ and recomputed $\Sigma _{n}(z)$

Alternatively

1. compute and store all $\Sigma _{n}(z)$
2. compute $E_{n}$ and color the points

An other approach looks like that:

Find the smallest $n$ so that $d_{k}(z,z+\varepsilon )$ is below a given bound for all $k\geq n$. If no such $n$ can be found, then assume $z$ to be an element of the Julia set. Otherwise, use $n$ to color the Fatou set.

This algorithm is a variant of the escape time algorithm (ETA). Note that in ETA the point does not really escape (at least if we are on the sphere), it just converges to ∞. This approach is similar. However, we don't need to know an attractive fixed point of ƒ.

Up to now, I didn't try the modified version. One disadvantage may be that the Fatou set will no more appear smooth colored. Then I am not sure if this modification is really an advantage, because the iteration must be done until a given maximum number of iterations is reached. Note that even if $d_{n}$ is under the bound for some $n$ the distance $d_{k}$ can rise again. I cannot say if this effect is crucial or can be neglected...

Gallery

• 
  $z\mapsto z^{2}$
• 
  $z\mapsto z^{2}+i$
• 
  $z\mapsto z^{2}-1$
• 
  ${}_{z\mapsto z^{2}-0.742+i/10}$
• 
  $z\mapsto N_{z^{3}-1}$
• 
  $z\mapsto N_{z^{3}-2z+2}$

$N$ denotes the Newton operator

$N:f\mapsto {\mathit {id}}-{\frac {f}{f'}}$

Using Critical Points Absorption (CPA)

The previous method yeilds neat, smooth colorings and requires least knowledge about the dynamics of the process. However, it is quite time consuming. Teh following approach is an extension of escape time algorithm (ETA) for polynomials.

Let ƒ be a polynomial of degree d > 1 over C. Such a polynomial has at most d critical points: infinity and the at most d–1 zeroes of ƒ′. It is well known that each attractor of the process z→ƒ(z) absorbs at least one critical point. Suppose z_K is a critical value, then ƒⁿ(z_K) comes arbitrarirly close to one of the attractive cycles of ƒ.

A process for a quadratic polynomial ƒ(z) = z² + c is the simplest case: The critical values are 0 and  ∞. As ∞ absorbs itself, only 0 is left, and we have the following cases. M denotes the Mandelbrot set.

$c\notin M$

Easy case. 0 is absorbed by ∞, and J(ƒ) consists just of Cantor dust. Escape time to ∞ can be used to color all points.

$c\in M\setminus \partial M$

If c is an element of the interior of M, then 0 will be absorbed by a (super) attractive cycle. Compute

$w=w(c,n)=f^{n}(0)$

for a sufficiently large n. As w (basically) only depends on c, it can be precomputed before starting the visualization of J(ƒ). Notice that w is the element of a cycle that might have more than one element, i.e. w is only unique modulo that cycle.

Coloring for a point z in C is then as easy as ETA:

• If z is absorbed by ∞, use escape time to color z.
• If ƒⁿ(z) comes close to w, i.e if |ƒⁿ(z) – w| < ε for the first time, use that n to color z.
• If ƒⁿ(z) neither comes close to ∞ nor comes close to w, then color z as part of J(ƒ). Only few z will fall into this category.

Visualising complex functions

Color space in HSV-projection

Suppose we want to visualise a complex valued function like

${\begin{aligned}f:\mathbb {C} &\;\to \mathbb {C} \\z&\;\mapsto w=f(z)\\\end{aligned}}$

In order to color $w$ we decompose it into its absolute value $|w|$ and its argument $\arg(w)$.

Then, we assign the color

$\mathrm {HSV} \left({\tfrac {1}{2\pi }}\arg(w),\;1-g_{a_{s},b_{s}}(|w|),\;g_{a_{v},b_{v}}(|w|)\right)$

to the point representing $z$. In this HSV color space, all values are in the range from 0 to  1. The first component (hue) of the HSV color depends only on the argument of $w$ and the second and third component (saturation and value) depend only on the absolute value of $w$. We use a transformation $g_{a,b}$ on $w$ to map it into the interval $[0,1]$. For $g$ see section helper functions.

Examples

The values indexed s (saturation) control the transition saturated colors→white (resp. gray scale), i.e. intermediate values → infinity. The values indexed h (value/valenz) control the transition black→bright, i.e. zero→nonzero. Parameter a controls where the transition takes place: a is just the radius of the circle dividing the two regions (dark/bright, saturated/gray, etc.) Parameter b controls how sharp the transition is: b small = soft, b large = sharp.

The following images all show the range [-10,10]×[-10,10]$i$ and use $a_{s}=5$ (radius of rainbow) and $a_{v}=1$ (radius of black disc).

In the image with swapped meanings of S and V zero is printed in white and infinity in black.

Gallery

• Some examples of colorings
• 
  $b_{s}=b_{v}=10$
• 
  $b_{s}=b_{v}=1$
• 
  $b_{s}=0.05,b_{v}=1/3$
• 
  $b_{s}=0.05,b_{v}=1/3$
  Roles of S and V interchanged

• Analysis of critical point of z→z²+c
• 
  n=1: starting values (color map)
• 
  n=2
• 
  n=3
• 
  n=4
• 
  n=17
• 
  n=18
• 
  n=19
• 
  n=20

Helper functions

Helper function t

some $t\,$

is a smooth, monotone sigmoidal transition from −1 to 1 that satisfies $t'(0)=1$ and $t(x)=-t(-x)$. There are many choices for it. Some of them are

${\begin{aligned}t:\mathbb {R} &\to (-1,1)\\x&\mapsto \tanh \,x\\&\mapsto {\tfrac {2}{\pi }}\arctan \left({\tfrac {\pi }{2}}x\right)\\&\mapsto {\tfrac {2}{\pi }}\operatorname {gd} \left({\tfrac {\pi }{2}}x\right)\\&\mapsto {\frac {x}{\sqrt {1+x^{2}}}}\\&\mapsto {\frac {x}{1+|x|}}\\\end{aligned}}$

with gd denoting the gudermannian function.

• All functions in a Desmos plot.

Helper function h

$h_{a}\,$

Helper function $h_{a}$ is almost the same like $t$ but it maps to $(0,1)$ and we can adjust the speed of transition by parameter $a$.

${\begin{aligned}h_{a}:\mathbb {R} &\to (0,1)\\x\;&\mapsto {\tfrac {1}{2}}+{\tfrac {1}{2}}t(2ax)\end{aligned}}$

with $a>0$. Then $h_{a}$ is symmetric to the point (0, 1/2) and

${\begin{alignedat}{2}h_{a}&(0)&&={\tfrac {1}{2}}\\h'_{a}&(0)&&=a\\h_{a}&(-\infty )&&=0^{+}\\h_{a}&(\infty )&&=1^{-}\\\end{alignedat}}$

Negative values are mapped to values between 0 and 1/2. Positive values are mapped to values between 1/2 and 1. The parameter $a$ controls how fast the transition will be.

If we want a falling function, we can use the symmetry

$h_{a}(x)\,=1-h_{a}(-x)=1-h_{-a}(x)$

i.e. we negate $x$ or $a$.

Helper function g

$g_{a,b}\,$

This function maps the positive numbers to the interval $[0,1)$.

${\begin{aligned}g_{a,b}:\mathbb {R} ^{\geqslant \,0}&\to (0,1)\\x\;&\mapsto {\tfrac {1}{2}}+{\tfrac {1}{2}}\,t(2b\cdot a\cdot u(x/a))\end{aligned}}$

for some function $u$ that is defined below. If $u$ is appropriately chosen then for $g$ the following holds

${\begin{array}{lll}g_{a,b}(0)&=&0\\g_{a,b}(\infty )&=&1\\g_{a,b}(a)&=&1/2\\g'_{a,b}(a)&=&b\\g'_{a,b}(x)&>&0{\text{ for }}x>0\end{array}}$

This means that $g_{a,b}$

• grows continuously from 0 to 1 as $x$ grows from $0$ to $\infty$
• we can control where $g$ crosses the middle between 0 and 1 by specifying parameter $a$
• we can control how fast $g$ passes the point $(a,1/2)$ by specifying parameter $b$

We are left with determining the finction $u$ on $\mathbb {R} ^{>0}$ with

$u\,$

$u$ must satisty

${\begin{array}{lll}u(0^{+})&=&-\infty \\u(\infty )&=&\infty \\u(1)&=&0\\u'(1)&=&1\end{array}}$

For $u$ we set

$u(x)={\tfrac {1}{2}}\left(x-{\tfrac {1}{x}}\right)$

Helper function w

$w_{a}\,$

This function maps the positive numbers to the interval $[0,1)$.

${\begin{aligned}w_{a}:\mathbb {R} ^{\geqslant \,0}&\to [0,1)\\x\;&\mapsto t(ax)\end{aligned}}$

By $a$ we can control its slope in the origin:

$w'(0)=\,a$

• Interactive Desmos plot of $w$

Circular arc through (0,0) and (1,1)

A circular arc through (0,0) and (1,1) that has given slope of $a$ in (0,0) and $-1/a$ in (1,1):

$\operatorname {arc} (x)={\begin{cases}x&;{\text{ if }}a=1\\y_{0}+\operatorname {sign} (\ln a){\sqrt {r^{2}-(x-x_{0})^{2}}}&;{\text{ if }}a\neq 1\\\end{cases}}$

where

${\begin{aligned}x_{0}&={\frac {a}{a-1}}\\y_{0}&=1-x_{0}\\r&={\sqrt {x_{0}^{2}+y_{0}^{2}}}\\\end{aligned}}$

The center of the circle is incident on the line $y=1-x$.

• Desmos plot

Bézier-Curves

Quadratic

Suppose we have a smooth function

${\begin{aligned}f:[t_{0},t_{1}]&\;\to \;\mathbb {R} ^{2}\\t&\;\mapsto \;(x(t),\,y(t))\\\end{aligned}}$

and like to draw an approximation of it using quadratic Béziers. Obviously, the two end points lie on ƒ and the control point lies on the crosspoint — for simplicity, we assume it exists — of the two tangents through the bounding points. In other words, we have to determine the intersection of the two lines

$g_{n}:\textstyle {\binom {x_{n}}{y_{n}}}+\lambda _{n}\cdot \textstyle {\binom {{\dot {x}}_{n}}{{\dot {y}}_{n}}}$

with canonical abbreviations like

${\dot {x}}_{n}={\tfrac {dx}{dt}}(t_{n})$

This leads to a determining equation for the λ's

${\binom {x_{0}-x_{1}}{y_{0}-y_{1}}}={\binom {-{\dot {x}}_{0}\quad {\dot {x}}_{1}}{-{\dot {y}}_{0}\quad {\dot {y}}_{1}}}{\binom {\lambda _{0}}{\lambda _{1}}}$

whose solution is

${\binom {\lambda _{0}}{\lambda _{1}}}={\frac {1}{{\dot {x}}_{1}{\dot {y}}_{0}-{\dot {x}}_{0}{\dot {y}}_{1}}}{\binom {{\dot {y}}_{1}\quad -{\dot {x}}_{1}}{{\dot {y}}_{0}\quad -{\dot {x}}_{0}}}{\binom {x_{0}-x_{1}}{y_{0}-y_{1}}}$

This gives us the intersection, i.e. the control point, by evaluating one of (or the two of) the g's at the end point(s).

Cubic

Desmos: Interactive, cubic Bézier

Suppose we have a smooth function

${\begin{aligned}f:[t_{0},t_{1}]&\;\to \;\mathbb {R} ^{2}\\t&\;\mapsto \;(x(t),\,y(t))\\\end{aligned}}$

and like to draw an approximation of it using cubic Béziers

${\begin{aligned}b(t)&=(1-t)^{3}P_{0}+3(1-t)^{2}tP_{1}+3(1-t)t^{2}P_{2}+t^{3}P_{3}\\&=P_{0}+3(P_{1}-P_{0})t+3(P_{2}-2P_{1}+P_{0})t^{2}+(P_{3}-3P_{2}+3P_{1}-P_{0})t^{3}\\&=P_{3}+3(P_{2}-P_{3})(1-t)+3(P_{1}-2P_{2}+P_{3})(1-t)^{2}+(P_{0}-3P_{1}+3P_{2}-P_{3})(1-t)^{3}\end{aligned}}$

From the two last representations we see immediately the value of the derivatives of b in the end points. We want the first and second derivatives in the end points to point in the same direction as the derivatives of ƒ do. Thus

This method applied to the standard parametrization (cos t, sin t) of a quarter of the unit circle (in black):
We get the coordinates (1, 1/2) and (1/2, 1) for the control points (red). The resulting bézier (orange) does not approximate the circle as good as possible.

${\begin{aligned}\alpha \cdot {\dot {f}}_{0}&=P_{1}-P_{0}\\\beta \cdot {\dot {f}}_{1}&=P_{3}-P_{2}\\\gamma \cdot {\ddot {f}}_{0}&=P_{0}-2P_{1}+P_{2}\\\delta \cdot {\ddot {f}}_{1}&=P_{1}-2P_{2}+P_{3}\\\end{aligned}}$

and together with ƒ₀=P₀ and ƒ₁=P₃ we get the linear system

${\begin{pmatrix}-2{\dot {x}}_{0}&-{\dot {x}}_{1}&-{\ddot {x}}_{0}&0\\-2{\dot {y}}_{0}&-{\dot {y}}_{1}&-{\ddot {y}}_{0}&0\\{\dot {x}}_{0}&2{\dot {x}}_{1}&0&-{\ddot {x}}_{1}\\{\dot {y}}_{0}&2{\dot {y}}_{1}&0&-{\ddot {y}}_{1}\end{pmatrix}}\cdot {\begin{pmatrix}\alpha \\\beta \\\gamma \\\delta \end{pmatrix}}={\begin{pmatrix}x_{0}-x_{1}\\y_{0}-y_{1}\\x_{1}-x_{0}\\y_{1}-y_{0}\end{pmatrix}}$

Note that in the special case of x(t) = t we get

$\alpha =\beta ={\tfrac {1}{3}}(t_{1}-t_{0})$

i.e. the x-coordinates of the control points P₁ and P₂ divide the interval [t₀,t₁] with ratios 1:2 resp. 2:1. This is quite astonishing because in order to get the control points we do not have to evaluate second derivatives of ƒ. This is due to properties of Bernstein polynomials.

To solve the above system we use standard technique like the adjugate.

Also note that if both second derivatives of x or both second derivatives of y happen to vanish the above system won't have full rank, i.e. the corresponding matrix won't be invertible. However, it's sufficient to determine α and β to get the control points and for vanishing second derivatives we get the less complicated systems

${\begin{pmatrix}-2{\dot {x}}_{0}&-{\dot {x}}_{1}\\{\dot {x}}_{0}&2{\dot {x}}_{1}\end{pmatrix}}\cdot {\binom {\alpha }{\beta }}={\binom {x_{0}-x_{1}}{x_{1}-x_{0}}}\qquad {\text{ if }}{\ddot {x}}_{0,1}=0$

resp.

${\begin{pmatrix}-2{\dot {y}}_{0}&-{\dot {y}}_{1}\\{\dot {y}}_{0}&2{\dot {y}}_{1}\end{pmatrix}}\cdot {\binom {\alpha }{\beta }}={\binom {y_{0}-y_{1}}{y_{1}-y_{0}}}\qquad {\text{ if }}{\ddot {y}}_{0,1}=0$

Solving this yields

${\binom {\alpha }{\beta }}={\frac {1}{3}}(x_{1}-x_{0}){\binom {{\dot {x}}_{0}}{{\dot {x}}_{1}}}\qquad {\text{ if }}{\ddot {x}}_{0,1}=0$

and ditto for y.

---

As you can see in the image above, there is some room for improvement and therefore we work out a second approach. Again, we set P₀=ƒ₀, P₃=ƒ₁ and constrain the two control points P₁ resp. P₂ to lie on the tangent through the end point next to it. This leaves us again with a two dimensional space to search in. Instead of imposing properties on second derivative, we now simply force the bézier to meet the curve ƒ in a third point Q, i.e. we set

$Q=f_{t}=f(t\cdot t_{0}+(1-t)\cdot t_{1})\quad {\stackrel {.}{=}}\quad b_{[P_{0},P_{1},P_{2},P_{3}]}(t)$

for some t in (0,1) and b denoting the bézier. The condition on the end and control points again reads as

This method applied to the standard parametrization (cos t, sin t) of a quarter of the unit circle (in black):
We get the coordinates (1, ω) and (ω, 1) for the control points (red) with
   ${\scriptstyle \omega ={\frac {4}{3}}({\sqrt {2}}-1)}$
The resulting bézier approximation (orange) is fine. There is no visual difference between the circle and the bézier. The additional point Q is indicated halfway on the bow.

${\begin{aligned}P_{0}&=f_{0}\\P_{1}&=f_{0}+\alpha \cdot {\dot {f}}_{0}\\P_{2}&=f_{1}-\beta \cdot {\dot {f}}_{1}\\P_{3}&=f_{1}\\\end{aligned}}$

Putting this together with the condition on Q we get

${\begin{aligned}&{\binom {x_{t}-x_{0}}{y_{t}-y_{0}}}+(2t^{3}-3t^{2}){\binom {x_{1}-x_{0}}{y_{1}-y_{0}}}\\&\qquad \qquad =3t(1-t){\binom {(1-t){\dot {x}}_{0}\quad -t{\dot {x}}_{1}}{(1-t){\dot {y}}_{0}\quad -t{\dot {y}}_{1}}}{\binom {\alpha }{\beta }}\end{aligned}}$

Note that we can use more than one point Q. Suppose we like to use n points Q_i to guide the bézier. Each Q_i will add two more lines in the above linear system, i.e. the system will look like

$w=M\!\cdot \!{\tbinom {\alpha }{\beta }}$

with a 2n-dimensional vector w and a 2n×2 matrix M. In general, this system is overdetermined and thus has no solution. Therefore, we solve the 2-dimensional system

$M^{\top }w=M^{\top }\!M\!\cdot \!{\tbinom {\alpha }{\beta }}$

instead which yields a least-square solution for α and β.

We make the above system a little bit more explicit for the case of one additional point at t = 1/2. The linear system then reduces to

${\binom {\alpha }{\beta }}={\frac {4}{3}}\,{\frac {1}{{\dot {x}}_{1}{\dot {y}}_{0}-{\dot {x}}_{0}{\dot {y}}_{1}}}{\binom {{\dot {y}}_{1}\quad -{\dot {x}}_{1}}{{\dot {y}}_{0}\quad -{\dot {x}}_{0}}}{\binom {x_{0}+x_{1}-2x_{1/2}}{y_{0}+y_{1}-2y_{1/2}}}$

Note that the same matrix already occured in the computation for quadratic bézier curves above, however the matrix now gets multiplied with a vector describing the curvature, whereas in the quadratic case it gets multiplied with a vector descriping the direction.

If the special case of a linear function x we get

$\alpha =\beta ={\frac {4}{3}}\cdot {\frac {y_{0}+y_{1}-2y_{1/2}}{{\dot {y}}_{1}-{\dot {y}}_{0}}}$

noname

Suppose we have a smooth function

${\begin{aligned}f:[t_{0},t_{3}]&\;\to \;\mathbb {R} ^{2}\\t&\;\mapsto \;(x(t),\,y(t))\\\end{aligned}}$

from the reals to the plane and like to draw an approximation of it using cubic Béziers.

Let u be a second function with fixed points at t₀ and t₃ which is smooth and monotone. Then the composition ƒ o u will yield exactly the same plot for any such u. Applying function u means that drawing the curve at different, arbitrarily increasing and descreasing speeds does not change the way the plot of the curve looks like. That is nice from the plotter's point of view. However, this generates some difficulties when we try to approximate the curve, which eventially turns out to be a playing field for calculus of variations that will lead to formulas and partial differential equations much too complex for practical considerations. So we look for some characteristics that are invariant under reparameterisations u.

The derivative of ƒ o u is the velocity

$v={\frac {\mathrm {d} }{\mathrm {d} t}}(f\circ u)={\dot {f}}\circ u\cdot {\dot {u}}$

This furmula is the rationale why the first derivative of the bézier curve shall point in the same direction as the first derivative of ƒ: the direction is independent of the parametrisation u. The speed v always points into the same direction, no matter how fast we drive along ƒ. This is no more true for the acceleration

$a={\dot {v}}={\frac {\mathrm {d} ^{2}}{\mathrm {d} t^{2}}}(f\circ u)={\dot {f}}\circ u\cdot {\ddot {u}}\,+\,{\dot {u}}^{2}\!\cdot \!{\ddot {f}}\circ u$

In the remainder of the essay we will only use properties of the curve at its end points. Thus, before we proceed, let's use the fact that u has fixed points in the t-valueas that yield the end points in order so simplify the formulas for v and $a$. The formulas then read

$v={\dot {f}}\cdot {\dot {u}}\quad {\text{ and }}\quad a={\dot {f}}\cdot {\ddot {u}}\,+\,{\dot {u}}^{2}\!\cdot \!{\ddot {f}}$

provided we evaluate these quandities at one of the end points.

We can think of $a$ as being composed of two orthogonal components: one in direction of motion that speeds up or slows down the pencil, and on perpendicular to v which leads to a change in direction. The projection of $a$ onto the direction normal to v is parallel to

${\begin{aligned}\langle a,v^{\bot }\rangle \cdot v^{\bot }&=\langle a,{\tbinom {~{\dot {y}}}{-{\dot {x}}}}\rangle \cdot {\tbinom {~{\dot {y}}}{-{\dot {x}}}}\\&={\binom {{\ddot {x}}\,{\dot {y}}^{2}\,-\,{\ddot {y}}\,{\dot {x}}\,{\dot {y}}}{{\ddot {y}}\,{\dot {x}}^{2}\,-\,{\ddot {x}}\,{\dot {y}}\,{\dot {x}}}}\cdot {\dot {u}}^{4}\\&\parallel {\binom {{\ddot {x}}\,{\dot {y}}^{2}\,-\,{\ddot {y}}\,{\dot {x}}\,{\dot {y}}}{{\ddot {y}}\,{\dot {x}}^{2}\,-\,{\ddot {x}}\,{\dot {y}}\,{\dot {x}}}}\end{aligned}}$

$L(a,b)=\int _{a}^{b}\!{\sqrt {{\dot {x}}^{2}+{\dot {y}}^{2}}}\,\mathrm {d} t$

$K(a,b)=\int _{a}^{b}{\frac {{\dot {x}}\,{\ddot {y}}\,-\,{\dot {y}}\,{\ddot {x}}}{{\dot {x}}^{2}+{\dot {y}}^{2}}}\,\mathrm {d} t$

Sphärengleiche Linear Transforms

Given a linear transfrom in n-dimensional euklidean Space:

${\begin{aligned}A:\quad \mathbb {R} ^{n}&\;\to \;\mathbb {R} ^{n}\\x&\;\mapsto \;A\cdot x\end{aligned}}$

We call two linear transforms A and B spärengleich if they map the unit sphere to the same set:

$A\sim B\quad \Leftrightarrow \quad A\cdot D_{n}=B\cdot D_{n}\qquad {\text{with}}\qquad D_{n}=\{x\in \mathbb {R} ^{n}\;/\;|x|\leqslant 1\}$

Obviously, this ~ is an equivalence relation and we look for a representant of each equivalence class. We observe that this relation preserves the spectral norm

$A\sim B\quad \Rightarrow \quad \|A\|_{2}=\|B\|_{2}$

and that orthogonal matrices don't change the equivalence class:

$Q^{\top }Q=\mathrm {id} \quad \Rightarrow \quad A\sim A\cdot Q$

for any A. The last line is immediately clear because orthogonal transformations map spheres to themselves. The preservation of spectral norm follows from the definition of spectral norm. Let

$A=U_{A}\cdot \Sigma _{A}\cdot V_{A}^{\top }$

be the singular value decomposition (svd) of A. Then we have

$A\sim U_{A}\cdot \Sigma _{A}$

and

$A\sim B\quad \Rightarrow \quad \Sigma _{A}=\Sigma _{B}$

However, the converse is not true.

Arcsin

arccos and its asymptote $\scriptstyle {\sqrt {2(1-x)}}$ in orange.

In order to approximate arcsin and arccos if the argument is close to 1 the following expansions might be useful:

${\begin{aligned}\arccos(1-x)&={\sqrt {2x}}\cdot a(x)\\\arcsin(1-x)&={\frac {\pi }{2}}-{\sqrt {2x}}\cdot a(x)\end{aligned}}$

with a rational power series

${\begin{aligned}a(x)&=1+{\frac {1}{12}}x+{\frac {3}{160}}x^{2}+{\frac {5}{896}}x^{3}+{\frac {35}{18432}}x^{4}+{\mathcal {O}}(x^{5})\\&=1+{\frac {1}{2}}\cdot {x \over 3\cdot 2}+{\frac {1\cdot 3}{2\cdot 4}}\cdot {x^{2} \over 5\cdot 2^{2}}+{\frac {1\cdot 3\cdot 5}{2\cdot 4\cdot 6}}\cdot {x^{3} \over 7\cdot 2^{3}}+\cdots \\&=\sum _{j=0}^{\infty }{\binom {2j}{j}}{\frac {x^{j}}{(2j+1)8^{j}}}\\&=\sum _{j=0}^{\infty }{\frac {(2j-1)!!}{(2j)!!}}{\frac {x^{j}}{(2j+1)\,2^{j}}}\end{aligned}}$

where !! denotes the double factorial. The radius of convergence of a is 2. We start by observing that

$\arccos(1-x)=2\arcsin {\bigl (}{\sqrt {x/2}}{\bigr )}$

which can easily verified by differentiaton. It follows that

${\begin{aligned}a(x)={\arccos(1-x) \over {\sqrt {2x}}}={\arcsin({\sqrt {x/2}}) \over {\sqrt {x/2}}}=\sum _{j=0}^{\infty }{\binom {2j}{j}}{\frac {x^{j}}{(2j+1)8^{j}}}\end{aligned}}$

Also note the following half-argument relations for –2 ≤ x ≤ 2:

${\begin{aligned}2\arcsin(x/2)&=\operatorname {sgn}(x)\arccos(1-x^{2}/2)\\2\arccos(x/2)&=\pi -\operatorname {sgn}(x)\arccos(1-x^{2}/2)\end{aligned}}$

N-th derivative

${\begin{aligned}{\frac {d^{n}}{dx^{n}}}\arcsin(x)&\;=\;\sum _{j=0}^{\left\lfloor {\frac {n-1}{2}}\right\rfloor }{\binom {n-1}{2j}}(2j-1)!!\,(2n-3-2j)!!\,{\frac {x^{n-1-2j}}{(1-x^{2})^{n-j-1/2}}}\\{\frac {d^{n}}{dx^{n}}}\arccos(x)&\;=\;-{\frac {d^{n}}{dx^{n}}}\arcsin(x)\end{aligned}}$

Again, !! denotes the double factorial. Note that for k < 0, k!! = 1.

Approximation of arcsin and arccos

The relative error of $x\cdot a(2x^{2})$ against $\arcsin(x)$ in $[0,1/2]$ stays below 6·10⁻¹⁸.

The relative error of ${\sqrt {2-2x}}\cdot a(1-x)$ against $\arccos(x)$ in $[1/2,1]$ stays below 6·10⁻¹⁸.

For 64-bit IEEE double computation (53-bit mantissa) of arcsin and arccos we use the following approach:

• If 1/2 ≤ |x| ≤ 1, compute

$\arccos(|x|)={\sqrt {2-2|x|}}\,a(1-|x|)$

• If |x| ≤ 1/2, compute

$\arcsin(x)=x\cdot a(2x^{2})$

We have to evaluate a(x) for x in [0, 1/2] and use the following rational MiniMax approximation of order [5/4] with a relative error below 6·10⁻¹⁸:

$a(x)\approx p(x)/q(x)$

with

p(x) = 
+ 0.99999999999999999442491073135027586203
- 1.0352340338921976278427312087167692142 x
+ 0.35290206232981519813422591897720574012 x^2
- 0.043334831706416857056123518013656946650 x^3
+ 0.0012557428614630796315205218507940285622 x^4
+ 0.0000084705471128435769021718764878041684288 x^5

q(x) =
+ 1
- 1.1185673672255329236623716486696411533 x
+ 0.42736600959872448854098334016758333519 x^2
- 0.063555884849631716599421483898013782858 x^3
+ 0.0028820878185134035637440105959294542908 x^4

(Note: A [4/5] rational MiniMax is even better, the relative error stays below 4.9·10⁻¹⁸.)

Then we use the symmetries

$\arccos(x)+\arcsin(x)=\pi /2$

$\arccos(x)+\arccos(-x)=\pi$

$\arcsin(x)+\arcsin(-x)=0$

of arcsin and arccos to get the desired result(s) over all of $[-1,1]$.

In order to achieve IEEE single precision, we can use

$a(x)\approx {\frac {45.210185257899-18.617417552712x+x^{2}}{45.210185141956-22.384922725383x+2.0175735681637x^{2}}}$

over $[0,0.501]$ with a relative error below 2.6·10⁻⁹.

The relative error for arccos stays below 2.6·10⁻⁹, whereas the one for arcsin stays below 5·10⁻⁹. See a Desmos plot of the relative errors.

Error Analysis

Using functional equations to get arcsin and arccos for all values in $[-1,1]$ makes the relative error more complicated than usual. Let $\delta _{0}$ denote the relative error of $a(x)$ at 0. Then we get the following envelopes for the relative errors (up to sign):

$\delta _{\arcsin }(x)=\delta _{0}\cdot {\begin{cases}1&;{\text{ if }}|x|<1/2\\{\displaystyle {\frac {\arccos(|x|)}{\arcsin x}}}&;{\text{ if }}|x|>1/2\\\end{cases}}$

$\delta _{\arccos }(x)=\delta _{0}\cdot {\begin{cases}{\displaystyle {\frac {\arccos(-x)}{\arccos x}}}&;{\text{ if }}x<-1/2\\{\displaystyle {\frac {\arcsin x}{\arccos x}}}&;{\text{ if }}|x|<1/2\\1&;{\text{ if }}x>1/2\\\end{cases}}$

Some special values

$x$	−1	−0.5−	−0.5+	0	0.5−	0.5+	1	$\max([-1,1])$
$\delta _{\arcsin }/\delta _{0}$	0	2	1	1	1	2	0	2
$\delta _{\arccos }/\delta _{0}$	0	0.5	0.25	0	0.5	1	1	1

• Desmos plot of relative errors with envelopes

Moduli Space of plane Triangles

Moduli space of similar, plane triangles is again a triangle (kind of).

• Interactive Desmos plot

Dobble Spot It

Dobble is a card game to have a nice time with fast pattern recognition. Each card shows 8 different icons, and any two cards have exactly one icon in common. The task is to find this common icon as fast as possible.

This essay is about how such a deck of cards can be constructed, and it supplies some mathematical background. As Dobble is famous for that mathematical background, we won't get too much into the theory; many articles found on the net address that background. A deck of card consists of 55 cards, each showing 8 different icons out of a set of 57 icons. So how do we have to arrange these icons on the cards so that any two cards picked at random have exactly one icon in common?

The property

any two cards have exactly one icon in common

reminds of a theorem in plane geometry:

any two different lines in the plane meet in exactly one point

Well, almost. In Euclidean geometry there are lines in the plane that don't meet, namely parallel lines. Before we go into details, let's summarize the geometric objects and how to associate to them a deck of cards, and how properties of the game arise from properties in geometry.

The two different Ways to identify Objects of Dobble with Objects in Geometry

Projective Plane of Order Q	Deck of Cards Cards = Lines,  Icons = Points	Deck of Cards Cards = Points,  Icons = Lines
Q2+Q+1 Lines	Q2+Q+1 Cards	Q2+Q+1 Icons
Q2+Q+1 Points	Q2+Q+1 Icons	Q2+Q+1 Cards
Each line passes through Q+1 points	Each card shows Q+1 icons	Each icon is shown on Q+1 cards
Each point is incident on Q+1 lines	Each icon is shown on Q+1 cards	Each card shows Q+1 icons
Any two different lines meet in exactly one point	Any two different cards have exactly one icon in common	Any two different icons are shown together on exactly one card
Any two different points uniquely determine one line	Any two different icons are shown together on exactly one card	Any two different cards have exactly one icon in common

Properties of Dobble

Dual properties which do only apply to a complete deck of Dobble with 57 cards. As mentioned below however, Dobble in incomplete as it comes with 55 cards only. Therefore, the dual properties do not hold for the Dobble you can buy.

In order to overcome the problem with parallel lines in Euclidean geometry, we switch to projective geometry which doen't come with that shortcoming. Whereas points in the Euclidean plane can be regarded as pairs of numbers (x,y), points in the projective plane are triples (x : y : z) such that not all three are equal to zero. In addition, we consider two triples P and P' as the same if they are multiples of each other, i.e. if there is a number λ ≠ 0 such that

$(x:y:z)=(\lambda x':\lambda y':\lambda z')$

A line g in the projective plane is given by a triple (g_x : g_y : g_z) and a point P = (x : y : z) is incident on the line provided

$g_{x}x+g_{y}y+g_{z}z=0$

This is similar to the Euclidean case where a point is incident on a line if

$g_{x}x+g_{y}y+g_{z}=0$

but in projective space the additional symmetry in the formulae for the point-on-a-line relation has its counterpart in the additional symmetry of any two distinct lines meet in exactly one point. Notice that if z ≠ 0 we can divide the projective condition by z, and the outcome is basically the condition for Euclidean space. However, in the projective case there are also points with z = 0 which don't have a counterpart in Euclidean geometry. These points are sometimes called the horizon or points at infinity, but we don't need such fuzzy wording or any distinction of different classes of points to construct a deck of Dobble.

The second difference to ordinary geometry is that a deck of cards consists only of finitely many items: a finite number of icons, a finite number of cards, and a finite number of icons per card. This is handled by considering geometry over a finite field instead of geometry over the Reals. A field is an entity which features addition and multiplication, both commutative and connected by the distributive law. The addition of any element can be undone, and the element that doesn't change the result of an addition is called zero and denoted as 0. The multiplication with any element except 0 can be undone, and the element that doesn't change the result of a multiplication is called one and denoted as 1.

So let's switch to a finite field F with Q elements. The first observation is that there are only finitely many points in the plane. There are Q³ triples of the form (x : y : z) with x, y and z in F. As not all three coordinates shall be zero, we are left with Q³−1 non-zero triples. Triples which are a multiple of each other are regarded as the same point, and because there are Q−1 non-zero values in F which can serve as the factor λ from above, we get the total number of points of the projective plane over the finite field F of order Q as

${\frac {Q^{3}-1}{Q-1}}=Q^{2}+Q+1$

This is also the number of lines in that plane because the lines are also represented as triples. In order to see how many points are incident on each line, let's enumerate the points as

$(0:0:1),\;(0:1:z){\text{ and }}(1:y:z)$

If the x-coordinate is non-zero, we can divide all coordinates by x which gets us the points of the form (1 : y : z). If x is zero and y is non-zero, then divide by y to get the points of the form (0 : 1 : z). If both x and y are zero, that point can be represented as (0 : 0 : 1) because z must be non-zero then.

In order to compute which points are incident on which line, we could just do brute force and iterate over all lines and all points and test whether or not the relation from above is satisfied. But we can do better by working out explicit formulae. To that end we use a different relation to determine whether a point is incident on a line:

$g_{x}z+g_{y}y+g_{z}x=0$

This is just a rearrangement of points which does not change the global structure. The advantage is that in our enumeration of points and lines, the first coordinates, i.e. x and g_x, are always 0 or 1 which will simplify the computation.

Lines and Points incident on it

Line g = (gx : gy : gz)	Constraint(s) on coordinates of P	P = (x : y : z)	Case
$(0:0:1)$	$x=0$	$(0:0:1),\;(0:1:z)$	(1)
$(0:1:g_{z})$	$y+g_{z}x=0$	$(0:0:1),\;(1:-g_{z}:z)$	(2)
$(1:g_{y}:g_{z})$	$z+g_{y}y+g_{z}x=0$	$(0:1:-g_{y}),\;(1:y:-g_{z}-g_{y}y)$	(3)

In either case, there are Q+1 points on each line, and it's easy to verify that any two distinct lines have exactly one point in common. Due to the duality of points and lines, each point is incident on Q+1 different lines.

The case Q = 2: Fano Plane

Image 1: The seven 3-icon cards represented as tri-color lines. The seven icons are represented as 7 colors:

(0:0:1)
(0:1:0)
(0:1:1)
(1:0:0)
(1:0:1)
(1:1:0)
(1:1:1)

Image 2: Seven 3-icon cards at the intersections of the seven lines. Each line represents an icon with the same color coding like in Image 1.

Let's work out the simplest case of Q = 2, the Fano plane. The field F is the Galois Field GF(2), the field with the 2 elements 0 and 1. We expect 2²+2+1 = 7 lines and hence also 7 points, each line passing through 2+1 = 3 points, and each point incident on 3 lines.

Fano Plane

Case	Line	Points	Cards
(1)	(0:0:1)	(0:0:1), (0:1:0), (0:1:1)
(2)	(0:1:0)	(0:0:1), (1:0:0), (1:0:1)
	(0:1:1)	(0:0:1), (1:1:0), (1:1:1)
(3)	(1:0:0)	(0:1:0), (1:0:0), (1:1:0)
	(1:0:1)	(0:1:0), (1:0:1), (1:1:1)
	(1:1:0)	(0:1:1), (1:0:0), (1:1:1)
	(1:1:1)	(0:1:1), (1:0:1), (1:1:0)

The case Q = 7: Dobble

For Q = 7 we get Dobble: 57 icons on 57 cards, each card displaying 8 icons. But wait — a deck of Dobble consists only of 55 cards, not of 7²+7+1 = 57. Why that? Nobody knows! Two cards are "missing" and not contained in the game. These two missing cards would show the following combinations of icons:

• snowman, exclamation mark, dog, eye, light bulb, ladybug, skull, hammer
• snowman, question mark, gingerbread man, maple leaf, cactus, daisy, ice cube, dino

Due to these two missing cards, some of the statements won't hold for Dobble: For example, not all icons are present 8-fold in the entire game. In particular, the snowman is only present 6 times as it is the icon common to the two missing cards. And not for any combination of two icons there is a card showing that combination. For example, there is no card showing both dog and eye because this combination belongs to one of the missing cards.

Beyond Dobble

A central object in the construction from above is the field F which only exists if Q is the power of a prime number p, i.e. Q =pⁿ for some prime number p and a natural number n ≥ 1. What happens if Q is not the power of a prime? We can use a rather axiomatic approach and define the projective plane of order Q to be an entity consisting of Q²+Q+1 points, Q²+Q+1 lines, each line passing through Q+1 points, each point incident on Q+1 lines, any two points uniquely determining a line, and any two lines having exactly one point in common.

For Q = 1, the construction from above still works provided we take the "1" in the formulae literally and set the free variable (z in cases (1) and (2), y in case (3)) to  0. That yields a game with 3 cards, each showing 2 icons out of a set of 3 icons. This works even though there is no field with one element.

For Q > 1, the Bruck–Ryser theorem adds some constraints on Q: If a projective plane of order Q exists and Q is 1 or 2 mod 4, then Q must be the sum of two squares. Hence, if Q is 1 or 2 mod 4 but not the sum of two squares, e.g. Q = 6, then no projective plane of order Q exists. However, there are infinitely many numbers remaining to which the theorem does not apply, the first one being 10 — the only case where an answer is known so far. The result for 10 has been achieved by heavy computation. The next greater value which is not the power of a prime and where the Bruck–Ryser theorem does not apply is Q = 12 with 13 icons per card and 157 cards. Taking a pure combinatorial approach we get

${\binom {157}{13}}={\frac {157!}{(157-13)!\cdot 13!}}=3\,393\,796\,168\,826\,188\,475\approx 3.3\cdot 10^{18}$

different 13-icon subsets (possible cards) which can be built out of the 157 icons, and for a complete game we have to pick 157 from these 3.3·10¹⁸ cards in such a way that all the axioms are satisfied.

Linear Recurrence

Suppose the linear recurrence

$x_{n}=a_{1}x_{n-1}+a_{2}x_{n-2}+\cdots +a_{k}x_{n-k}=\sum _{j=1}^{k}a_{j}x_{n-j}$

for $n>k\geqslant 1$ where $x_{1}$, ..., $x_{k}$ are given numbers.

We want to determine an explicit representation of $x_{n}$. To that end, write the recurrence as:

$\underbrace {\left({\begin{array}{l}x_{n\;\;\;}\\x_{n-1}\\\;\;\vdots \\x_{n-k+2}\\x_{n-k+1}\\\end{array}}\right)} _{\displaystyle {=:y_{n}}}=\underbrace {\begin{pmatrix}a_{1}&a_{2}&\cdots &a_{k-1}&a_{k}\\1&0&\cdots &0&0\\0&1&\cdots &0&0\\\vdots &&\ddots &&\vdots \\0&0&\cdots &1&0\\\end{pmatrix}} _{\displaystyle {=:A\in K^{k\times k}}}\cdot \underbrace {\left({\begin{array}{l}x_{n-1}\\x_{n-2}\\\;\;\vdots \\x_{n-k+1}\\x_{n-k}\\\end{array}}\right)} _{\displaystyle {=:y_{n-1}}}$

so that it takes the form

$y_{n}=Ay_{n-1}=A^{n-k}y_{k}$

Now suppose $A$ has $k$ different eigenvectors $v_{j}$ and we know all of them, including the corresponding eigenvalues $\lambda _{j}$. Then we can write:

$y_{k}=\sum _{j=1}^{k}\beta _{j}v_{j}=V{\begin{pmatrix}\beta _{k}\\\vdots \\\beta _{1}\end{pmatrix}}={\begin{pmatrix}v_{k}&\cdots &v_{1}\end{pmatrix}}\cdot {\begin{pmatrix}\beta _{k}\\\vdots \\\beta _{1}\end{pmatrix}}$

where the $\beta _{j}$ are scalars in the algebraic closure of $K$ and $V$ is a matrix with the eigenvectors of $A$ as columns. Hence:

$y_{n}=A^{n-k}y_{k}=A^{n-k}{\Big (}\sum _{j=1}^{k}\beta _{j}v_{j}{\Big )}=\sum _{j=1}^{k}\beta _{j}A^{n-k}v_{j}=\sum _{j=1}^{k}\beta _{j}\lambda _{j}^{n-k}v_{j}\qquad (1)$

which leaves is with the computation of the $\beta _{j}$, the $v_{j}$ and the $\lambda _{j}$. Once we determined the eigenvectors, we get the $\beta _{j}$ by means of:

${\begin{pmatrix}\beta _{k}\\\vdots \\\beta _{1}\end{pmatrix}}=V^{-1}y_{k}$

Expanding the determinant of $A-\lambda E$ by expanding after it's top row, we find that all eigenvalues satisfy the characteristic equation

$\lambda ^{k}=\sum _{j=1}^{k}a_{j}\lambda ^{k-j}=a_{1}\lambda ^{k-1}+a_{2}\lambda ^{k-2}+\cdots +a_{k-1}\lambda +a_{k}$

From this we easily see that the eigenvectors of $A$ are:

$v_{j}=\left({\begin{array}{l}\lambda _{j}^{k-1}\\\;\;\vdots \\\lambda _{j}^{2}\\\lambda _{j}\\1\\\end{array}}\right)$

Due to (1), in order to get $x_{n}$ we take the top component of $y_{n}$ to get:

$x_{n}=\sum _{j=1}^{k}\beta _{j}\lambda _{j}^{n-k}\lambda _{j}^{k-1}=\sum _{j=1}^{k}\beta _{j}\lambda _{j}^{n-1}\qquad (2)$

Thus we are finished: Depending on the $a_{j}$, the eigenvalues can be computed explicitly or by numerical methods. From the eigenvalues we get the matrix $V$ which we use to compute the coefficients $\beta _{j}$ from the starting values $x_{1}$ ... $x_{k}$ so that we have determines all unknowns in (2).

Example 1

Let $a_{1}=a_{2}=1$ so that we get the Fibonacci sequence

$x_{n}=x_{n-1}+x_{n-2}$

with the characteristic equation

$\lambda ^{2}=\lambda +1$

Finding solutions of ln cosh x = αx + β

The graph of

$f(x)=\ln \cosh x$

is convex and loosely resembles a hyperbola

$h_{c}(x)={\sqrt {c^{2}+x^{2}}}-c$

In particular, $f(x)\approx |x|-\ln 2$ for large $|x|$. The convexity of $f$ allows for easy specification of the number $S(\alpha ,\beta )$ of solutions of

$f(x)~{\stackrel {!}{=}}~\alpha x+\beta$		(1)

If the line $\alpha x+\beta$ is tangent to $f$, then there is a uniqe solution of multiplicity two. If $|\alpha |>1$ then (1) has a unique solution, and for $|\alpha |\leqslant 1$ the number of solutions depend on whether the line runs below a tangent to $f$ or above:

$S(\alpha ,\beta )={\begin{cases}1,&{\text{if }}|\alpha |>1\\1,&{\text{if }}|\alpha |=1{\text{ and }}\beta >-\ln 2\\0,&{\text{if }}|\alpha |=1{\text{ and }}\beta \leqslant -\ln 2\\0,&{\text{if }}|\alpha |<1{\text{ and }}\beta <g(\alpha )\\1_{2},&{\text{if }}|\alpha |<1{\text{ and }}\beta =g(\alpha )\\2,&{\text{if }}|\alpha |<1{\text{ and }}\beta >g(\alpha )\\\end{cases}}$

with

$g(x)=-{\frac {1}{2}}\ln \left(1-x^{2}\right)-x\operatorname {artanh} x$

and where $1_{2}$ stand for "one solution of multiplicity two".

The procedure to approximate solutions of (1) is then as follows:

1. Compute $S=S(\alpha ,\beta )\in \{0,1,2\}$, the number of solutions of (1). If $S=0$, then terminate with "no solution".
2. Compute the solutions of

   $h_{c}(x)~{\stackrel {!}{=}}~\alpha x+\beta$		(2)

   for $c=1$. The hyperbola $h_{1}$ has the property that (2) has at least as many solutions like (1) has. If $S=1$ but (2) has two solutions $x_{1}$ and $x_{2}$, then use $(x_{1}+x_{2})/2$. The solutions of (2) are the solutions of

   $Ax^{2}+Bx+C=0\quad {\text{ that satisfy }}\quad \alpha x+\beta \geqslant -c$

   where

   ${\begin{aligned}A&=1-\alpha ^{2}\\B&=-2\alpha (\beta +c)\\C&=c^{2}-(b+c)^{2}\end{aligned}}$

   This is a quadratic equation in $x$, except when $A=0$ where it is linear with solution $x=-C/B$. In the quadratic case, the number of solutions are 0, 1 or 2, depending on whether the discriminant $D=B^{2}-4AC$ satisfies $D<0$, $D=0$ or $D>0$, respectively:

   $x_{1,2}={\frac {-B\pm {\sqrt {D}}}{2A}}$

3. We determined $S\in \{1,2\}$ approximate solutions in step 2. For each initial approximation $x$, perform Newton-Raphson iterations for (1) until the desired accuracy is reached:

   $x\mapsto x-{\frac {f(x)-\alpha x-\beta }{f'(x)-\alpha }}=x-{\frac {\ln(\cosh x)-\alpha x-\beta }{\tanh(x)-\alpha }}$

Finding solutions of a^x + b^x = c

For $a,b,c\in {\mathbb {R}}^{+}$ we want to determine approximate solutions $x\in {\mathbb {R}}$ of

$a^{x}+b^{x}=c$		(1)

Before handling the general case, sort out the special case $a=b$:

• If $a=b\neq 1$, then there is a single solution $x=\log _{a}(c/2)$.

• If $a=b=1$, $c=2$ we can take any $x\in {\mathbb {R}}$, and if $a=b=1$, $c\neq 2$ then there are no solutions.

In the general case of $a\neq b$, transform $a^{x}+b^{x}=c$ to the equivalent problem

$\ln \cosh x=\alpha x+\beta$		(2)

handled in the previous section, where

${\begin{aligned}\alpha &={\frac {\ln a+\ln b}{\ln a-\ln b}}\\\beta &=\ln(c/2)\end{aligned}}$

The solutions of (1) are then given by

${\frac {2x}{\ln b-\ln a}}$

where $x$ ranges over all solutions of (2).

Finding solutions of a^x − b^x = c

For $a,b\in {\mathbb {R}}^{+}$, $c\in {\mathbb {R}}$ we want to determine approximate solutions $x\in {\mathbb {R}}$ of

$a^{x}-b^{x}=c$		(1)

If $c=0$, then the solution is $x=0$, except when $a=b$ which allows any $x$.

If $c<0$, then we solve the equivalent problem $b^{x}-a^{x}=-c$ with $-c>0$. Thus, without loss of generality, we may assume $c>0$ in the remainder.

Before treating the general case, sort out the special cases $a=b$ and $b=1$:

• If $a=b$, then there is no restriction on $x$ if $c=0$, and if $c\neq 0$ there is no solution.

• If $a\neq b=1$, then the solution is $x=\log _{a}(1+c)$.

In order to solve the general case $b\neq 1$, transform $a^{x}-b^{x}=c$ to the equivalent problem

$\ln \cosh x=\alpha x+\beta$		(2)

handled in a previous section, where

${\begin{aligned}\alpha &=1-2\log _{b}a\\\beta &=(\log _{b}a-1)\ln c\end{aligned}}$

The solutions of (1) are then given by

${\frac {\ln c-2x}{\ln b}}$

where $x$ ranges over all solutions of (2).

Finding solutions of a^x + b^x = c^x

For $a,b,c\in {\mathbb {R}}^{+}$ we want to determine approximate solutions $x\in {\mathbb {R}}$ of

$a^{x}+b^{x}=c^{x}$		(1)

• If $a=c$ or $b=c$, then there is no solution.

• If $a=b\neq c$, then the solution is $x=\ln 2/(\ln a-\ln b)$.

In order to solve the general case, transform $a^{x}+b^{x}=c^{x}$ to the equivalent problem

$\ln \cosh x=\alpha x+\beta$		(2)

handled in a previous section, where

${\begin{aligned}\alpha &={\frac {2\ln c-\ln a-\ln b}{\ln a-\ln b}}\\\beta &=-\ln 2\end{aligned}}$

The solutions of (1) are then given by

${\frac {2x}{\ln a-\ln b}}$

where $x$ ranges over all solutions of (2).

Finding solutions of x^r = ax + b

For $r,a,b\in {\mathbb {R}}$ we want to determine approximate solutions $x\in {\mathbb {R}}$ of

$x^{r}=ax+b$		(1)

Split the problem into several special cases. Once a special case has been treated, the assumption is that it won't occur in the cases below.

r = 0

• $a=0$: All possible $x$ are solutions if $b=1$. If $b\neq 1$, there are no solutions.
• $a\neq 0$: The solution is $x=(1-b)/a$.

r = 1

• $a=1$: All possible $x$ are solutions if $b=0$. If $b\neq 0$, there are no solutions.
• $a\neq 1$: The solution is $x=b/(1-a)$.

r < 0

Solve the equivalent problem

$y^{1-r}=by+a$		(2)

where $1-r>1$, a case handled below. Solutions of (1) are then of the form $x=1/y$, where $y$ iterates over all non-zero solutions of (2).

In the remainnder, we will have to consider three different kinds of exponents $r$:

• Odd exponents: The exponent has a representation $r=p/q$ where $p$ and $q$ are odd integers.
• Even exponents: The exponent has a representation $r=p/q$ where $p$ is an even integer and $q$ is an odd integer.
• Real exponents: The exponent neither fits the odd case nore the even case.

a = 0

The solutions of (1) are the solutions of

$x^{r}=b$		(3)

which are:

$x={\begin{cases}b^{1/r},&{\text{ if }}r{\text{ is odd}}\\b^{1/r},&{\text{ if }}b\geqslant 0{\text{ and }}r{\text{ is real}}\\\pm b^{1/r},&{\text{ if }}b\geqslant 0{\text{ and }}r{\text{ is even}}\\{\text{none}},&{\text{ if }}b<0{\text{ and }}r{\text{ is even or real}}\\\end{cases}}$		(4)

0 < r < 1

Set $x=y^{1/r}$ and solve the problem

$y^{1/r}={\frac {1}{a}}y-{\frac {b}{a}}$		(5)

with $1/r>1$. Solutions of (1) are then $x=y^{1/r}$ where $y$ iterates over all solutions of (5). If $r$ is not odd, then the solutions must satisfy $x\geqslant 0$.

In addition, for even exponents there are solutions $x=-y^{1/r}$ that satisfy $x<0$, where $y$ ranges over solutions of

$y^{1/r}=-{\frac {1}{a}}y-{\frac {b}{a}}$		(6)

r > 1

All other cases have been reduced to this one (or have been handled as special cases). Noice that we have $a\neq 0$. In order to simplify this presentation, assume that we only have odd and even exponents. The real case is mapped to the even case by assuming

$f(x)=x^{r}$		(7)

is an even function. After we computed all solutions, we ditch negative solutions if the came from the real case.