View source for Piecewise linear continuation

'''Simplicial continuation''', or '''piecewise linear continuation''' (Allgower and Georg),<ref name="one">Eugene L. Allgower, K. Georg, "Introduction to Numerical Continuation Methods", ''SIAM Classics in Applied Mathematics'' 45, 2003.</ref><ref name="three">E. L. Allgower, K. Georg, "Simplicial and Continuation Methods for Approximating Fixed Points and Solutions to Systems of Equations", ''SIAM Review'', Volume 22, 28-85, 1980.</ref> is a one-parameter [[numerical continuation|continuation method]] which is well suited to small to medium embedding spaces. The algorithm has been generalized to compute higher-dimensional manifolds by (Allgower and Gnutzman)<ref name="two">Eugene L. Allgower, Stefan Gnutzmann, "An Algorithm for Piecewise Linear Approximation of Implicitly Defined Two-Dimensional Surfaces", ''SIAM Journal on Numerical Analysis'', Volume 24, Number 2, 452-469, 1987.</ref> and (Allgower and Schmidt).<ref name="four">Eugene L. Allgower, Phillip H. Schmidt, "An Algorithm for Piecewise-Linear Approximation of an Implicitly Defined Manifold", ''SIAM Journal on Numerical Analysis'', Volume 22, Number 2, 322-346, April 1985.</ref>

The algorithm for drawing contours is a simplicial continuation algorithm, and since it is easy to visualize, it serves as a good introduction to the algorithm.

== Contour plotting ==

The contour plotting problem is to find the zeros (contours) of <math> f(x,y)=0\,</math> (<math> f(\cdot)\,</math> a smooth scalar valued function) in the square <math>0\leq x \leq 1, 0\leq y \leq 1\,</math>,
<center>
[[Image:Contours.gif|An example of contours]] [[Image:ContoursA.gif|Contours, three-dimensional view]]
</center>
The square is divided into small triangles, usually by introducing points at the corners of a regular square mesh <math>ih_x\leq x\leq (i+1)h_x\,</math>, <math>jh_y\leq y \leq (j+1)h_y\,</math>, making a table of the values of <math>f(x_i,y_j)\,</math> at each corner <math>(i,j)\,</math>, and then dividing each square into two triangles. The value of <math>f(x_i,y_j)\,</math> at the corners of the triangle defines a unique Piecewise Linear interpolant <math>lf(x,y)\,</math> to <math>f(\cdot)\,</math> over each triangle. One way of writing this interpolant on the triangle with corners
<math>(x_0,y_0),~(x_1,y_1),~(x_2,y_2)\,</math> is as the set of equations

: <math> (x,y) = (x_0,y_0)+(x_1-x_0,y_1-y_0)s+(x_2-x_0,y_2-y_0)t\,</math>
: <math> 0\leq s\,</math>
: <math> 0\leq t\,</math>
: <math> s+t \leq 1\,</math>
: <math> lf(x,y) = f(x_0,y_0)+(f(x_1,y_1)-f(x_0,y_0))s+(f(x_2,y_2)-f(x_0,y_0))t\,</math>

The first four equations can be solved for <math>(s,t)\,</math> (this maps the original triangle to a right unit triangle), then the remaining equation gives the interpolated value of <math>f(\cdot)\,</math>. Over the whole mesh of triangles, this piecewise linear interpolant is continuous.
<center>
[[Image:LinearInterpolant.gif|An example of a triangulation and marked vertices]] [[Image:LinearInterpolantA.gif|Linear interpolant, three-dimensional view]]
</center>
The contour of the interpolant on an individual triangle is a line segment (it is an interval on the intersection of two planes). The equation for the line can be found, however the points where the line crosses the edges of the triangle are the endpoints of the line segment.
<center>
[[Image:Contour.gif|The unique linear interpolant on a simplex and its zero set]][[Image:TriangleContour.gif|The contour of the linear interpolant over a triangle]]
</center>
The contour of the piecewise linear interpolant is a set of curves made up of these line segments. Any point on the edge connecting <math>(x_0,y_0)\,</math> and <math>(x_1,y_1)\,</math> can be written as

: <math>(x,y) = (x_0,y_0) + t (x_1-x_0,y_1-y_0),\,</math>

with <math>t\,</math> in <math>(0,1)\,</math>, and the linear interpolant over the edge is

: <math> f \sim f_0 + t (f_1-f_0)\,</math>

So setting <math> f = 0\,</math>

: <math>t = -f_0/(f_1-f_0)\,</math> and <math> (x,y) = (x_0,y_0)-f_0*(x_1-x_0,y_1-y_0)/(f_1-f_0)\,</math>

Since this only depends on values on the edge, every triangle which shares this edge will produce the same point, so the contour will be continuous. Each triangle can be tested independently, and if all are checked the entire set of contour curves can be found.

== Piecewise linear continuation ==

Piecewise linear continuation is similar to contour plotting (Dobkin, Levy, Thurston and Wilks),<ref name="five">[[David P. Dobkin]], Silvio V. F. Levy, [[William Thurston|William P. Thurston]] and Allan R. Wilks, "Contour Tracing by Piecewise Linear Approximations", ''ACM Transactions on Graphics'', 9(4) 389-423, 1990.</ref> but in higher dimensions. The algorithm is based on the following results:

=== Lemma 1 ===

<center>
{| class="wikitable"
|-
| If F(x) maps <math>\mathbb{R}^n</math> into <math>\mathbb{R}^{n-1}</math>, there is a unique linear interpolant on an '(n-1)'-dimensional [[simplex]] which agrees with the function values at the vertices of the simplex.
|}
</center>

An '(n-1)'-dimensional simplex has n vertices, and the function F assigns an 'n'-vector to each. The simplex is [[convex set|convex]], and any point within the simplex is a [[convex combination]] of the vertices. That is:

If x is in the interior of an (n-1)-dimensional simplex with n vertices <math> v_i </math>, then there are positive scalars <math>0<\alpha_i</math> such that

<center>
: <math> \mathbf{x} = \sum_i \alpha_i \mathbf{v}_i </math>
: <math> \sum_i \alpha_i = 1.\,</math>
</center>

If the vertices of the simplex are [[Linear independence|linearly independent]] the non-negative scalars <math>\alpha</math> are unique for each point x, and are called the [[Barycentric coordinates (mathematics)|barycentric coordinates]] of x. They determine the value of the unique [[interpolation|interpolant]] by the formula:
<center>
: <math>LF = \sum_i \alpha_i F(\mathbf{v}_i)</math>
</center>

=== Lemma 2 ===

<center>
{| class="wikitable"
|-
| An (n-1)-dimensional simplex can be tested to determine if it contains the origin.
|}
</center>

There are basically two tests. The one which was first used labels the vertices of the simplex with a vector of signs (+/-) of the coordinates of the vertex. For example the vertex (.5,-.2,1.) would be labelled (+,-,+). A simplex is called ''completely labelled'' if there is a vertex whose label begins with a string of "+" signs of length 0,1,2,3,4,...n. A completely labelled simplex contains a neighborhood of the origin. This may be surprising, but what underlies this result is that for each coordinate of a completely labelled simplex there is a vector with "+" and another with a "-". Put another way, the smallest cube with edges parallel to the coordinate axes and which covers the simplex has pairs of faces on opposite sides of 0. (i.e. a "+" and a "-" for each coordinate).

The second approach is called ''vector labelling''. It is based on the barycentric coordinates of the vertices of the simplex. The first step is to find the barycentric coordinates of the origin, and then the test that the simplex contains the origin is simply that all the barycentric coordinates are positive and the sum is less than 1.

=== Lemma 3 ===

<center>
{| class="wikitable"
|-
| There is a triangulation (the Coxeter-Freudenthal-Kuhn triangulation [1]) which is invariant under the pivot operation
<center><math>P_{v_i} (v_0,v_1,...,v_n) := (v_0,v_1,...,v_{i-1},\tilde v_i,v_{i+1},...,v_n)</math></center>
where
<center>
<math>
\tilde v_i = \left\{ \begin{array}{lcl}
v_1+v_n-v_0 & &i=0 \\
v_{i+1}+v_{i-1}-v_i&\qquad\qquad&0\\
v_{n-1}+v_0-v_n & &i=n \\ \end{array}\right.
</math>
</center>
|}
</center>

<center>
[[Image:Simplical3dOne.gif|The first step of three-dimensional simplicial continuation]] [[Image:Simplical3dTwo.gif|The second step of three-dimensional simplicial continuation]]
</center>

<center>
[[Image:Simplicial.gif]]
</center>

== References ==
{{reflist}}

[[Category:Numerical analysis]]