View source for Calculus/Discrete vector calculus

This chapter will present an analog to vector calculus where space now consists of discrete "lumps". The purpose of this chapter is to provide an intuitive basis for vector calculus. The portrayal of vector calculus in this chapter will enable the generalization of vector calculus to '''non-Euclidean''' geometries.

== The essence of Calculus ==

The key approach to calculus is to start with a model where each variable, using <math>x</math> as an example, takes on discrete quantities separated by small intervals of width <math>\Delta x</math>. <math>\Delta x</math> is a measure of the model's "coarseness" (when <math>\Delta x</math> is increasing) or "fineness" (when <math>\Delta x</math> is decreasing). When handling quantities at the scale of <math>\Delta x</math>, an approximation for which the error '''[[w:Relative change and difference|relative]]''' to <math>\Delta x</math> vanishes as <math>\Delta x</math> becomes infinitely small, becomes exact when <math>x</math> is treated as a continuous variable.

As an example, let <math>\Delta x</math> be small, and restrict <math>x</math> to the quantities <math>\{..., -3\Delta x, -2\Delta x, -\Delta x, 0, +\Delta x, +2\Delta x, +3\Delta x, ...\}</math>. Let <math>x_i = i\Delta x</math> for any integer <math>i</math>.

The difference <math>\Delta y(i) = x_{i+1}^2 - x_i^2 = 2x_i\Delta x + \Delta x^2</math> approaches 0 as <math>\Delta x \to 0^+</math>, but it is not accurate to approximate <math>\Delta y(i)</math> with 0, as the error relative to <math>\Delta x</math> does not vanish as <math>\Delta x \to 0^+</math>: in truth, <math>\frac{\Delta y(i) - 0}{\Delta x} = 2x_i + \Delta x \not\to 0</math> as <math>\Delta x \to 0</math>.

It is however, accurate to approximate <math>\Delta y(i)</math> with <math>\Delta z(i) = 2x_i\Delta x</math> since the error relative to <math>\Delta x</math> does vanish as <math>\Delta x \to 0^+</math>: it is the case that <math>\frac{\Delta y(i) - \Delta z(i)}{\Delta x} = \Delta x \to 0</math> as <math>\Delta x \to 0</math>.

When we articulate continuous space into small but non-zero volumes for the purpose of approximating vector calculus, it is important to note that the '''only''' approximations that will be made, will be those that entail relative errors that vanish as the articulation of space becomes more and more fine.

In order to demonstrate the necessity of having the error relative to the differential <math>\Delta x</math> approach 0 as <math>\Delta x \to 0</math> as opposed to the absolute error, the following example will be used:

Consider the closed interval <math>[a, b]</math> divided into <math>N</math> intervals: <math>[x_0, x_1], [x_1, x_2], ..., [x_{N-1}, x_N]</math> where <math>a = x_0 < x_1 < ... < x_N = b</math>. From each <math>i = 1, 2, \dots, N</math> choose a representative point <math>x_i^* \in [x_{i-1}, x_i]</math>, and let <math>\Delta x_i = x_i - x_{i-1}</math>. Let <math>f : \R^2 \to \R</math> and <math>g : \R^2 \to \R</math> be two functions whose first parameter is a representative point, and the second parameter is the interval width. Both <math>f</math> and <math>g</math> will be assumed to approach 0 as the interval width becomes infinitely small. Consider the two sums <math>F = \sum_{i=1}^N f(x_i^*,\Delta x_i)</math> and <math>G = \sum_{i=1}^N g(x_i^*,\Delta x_i)</math>. It will be assumed that <math>F</math> and <math>G</math> converge to <math>F_\infty</math> and <math>G_\infty</math> respectively as <math>N \to +\infty</math> and <math>\max(\Delta x_i) \to 0^+</math>.

The question of interest are circumstances under which <math>F_\infty = G_\infty</math>. If other words, as the interval <math>[a, b]</math> becomes more and more articulated into smaller and smaller intervals, what is required so that approximating <math>F</math> with <math>G</math> becomes exact?

If <math>\max_{x \in [a,b]}|f(x, \Delta x) - g(x, \Delta x)| \to 0</math> as <math>\Delta x \to 0^+</math>, this is not sufficient to guarantee that <math>F_\infty = G_\infty</math>. In contrast, if <math>\max_{x \in [a,b]}\frac{|f(x, \Delta x) - g(x, \Delta x)|}{\Delta x} \to 0</math> as <math>\Delta x \to 0^+</math>, this will force <math>F_\infty = G_\infty</math>. To prove this, the following derivation will be used:

<math>|F - G| = \left|\sum_{i=1}^N f(x_i^*,\Delta x_i) - \sum_{i=1}^N g(x_i^*,\Delta x_i)\right| = \left|\sum_{i=1}^N (f(x_i^*,\Delta x_i) - g(x_i^*,\Delta x_i))\right|</math>

<math> \leq \sum_{i=1}^N |f(x_i^*,\Delta x_i) - g(x_i^*,\Delta x_i)| = \sum_{i=1}^N \Delta x_i\frac{|f(x_i^*,\Delta x_i) - g(x_i^*,\Delta x_i)|}{\Delta x_i}</math>

<math> \leq \sum_{i=1}^N \Delta x_i\max_j\left(\frac{|f(x_j^*,\Delta x_j) - g(x_j^*,\Delta x_j)|}{\Delta x_j}\right) </math>
<math> = \left(\sum_{i=1}^N \Delta x_i\right)\max_j\left(\frac{|f(x_j^*,\Delta x_j) - g(x_j^*,\Delta x_j)|}{\Delta x_j}\right) </math>
<math> = (b-a)\max_j\left(\frac{|f(x_j^*,\Delta x_j) - g(x_j^*,\Delta x_j)|}{\Delta x_j}\right) </math>

For each <math>j = 1, 2, \dots, N</math>, it is assumed that <math>\frac{|f(x_j^*,\Delta x_j) - g(x_j^*,\Delta x_j)|}{\Delta x_j} \to 0</math> as <math>\Delta x_j \to 0^+</math>. Therefore <math>\max_j\left(\frac{|f(x_j^*,\Delta x_j) - g(x_j^*,\Delta x_j)|}{\Delta x_j}\right) \to 0</math> and hence <math>|F - G| \to 0</math> as <math>N \to +\infty</math> and <math>\max(\Delta x_i) \to 0^+</math>. Therefore <math>F_\infty = G_\infty</math>.

[[File:Riemann sums.svg|thumb|center|upright=4.0|The left image shows the sums <math>F</math> and <math>G</math>. The white area of each strip depicts <math>f</math>, and the red area of each strip depicts <math>g-f</math>. The total area highlighted in red depicts the difference between <math>G</math> and <math>F</math>. As <math>N</math> increases, the center image shows how the total red area does not decrease if it is not the case that for each strip that the "ratio of red area to the strip width", which is the height of the red area, fails to approach 0 for all strips. This is despite the fact that the red area (absolute difference) in each strip is approaching 0. The right image shows how the total red area approaches 0 if the "red area to strip width ratio" for each strip is approaching 0.]]

== Lumped space ==

=== Directed graphs ===

The lumped approximation of space will require the use of a "directed graph".

A "directed graph" <math>\mathcal{G} = \langle \mathcal{N}, \mathcal{E} \rangle</math> consists of a set <math>\mathcal{N}</math> of nodes/vertices and a set <math>\mathcal{E}</math> of directed edges. For each edge <math>e \in \mathcal{E}</math>, <math>e</math> has an origin node <math>t_\bullet(e) \in \mathcal{N}</math> that is the beginning of edge <math>e</math>, and a terminal node <math>t^\bullet(e) \in \mathcal{N}</math> that is the end of edge <math>e</math>. It will also be allowed that edges can form loops (<math>t_\bullet(e) = t^\bullet(e)</math>) and that multiple edges can have both the same beginning and end (<math>t_\bullet(e_1) = t_\bullet(e_2)</math> and <math>t^\bullet(e_1) = t^\bullet(e_2)</math>), though this will be irrelevant for the cases that will be considered.

An arbitrary function <math>f: \mathcal{N} \to \R</math> over the set of nodes will be referred to as a "'''node based function'''", and is effectively a scalar field. The notation <math>\langle n \in \mathcal{N} : f(n) \rangle</math> will be used if the input parameter is ambiguous. A node based function can also be denoted by exhaustively listing the assignment to each node: <math>\{n_1 \mapsto f(n_1), n_2 \mapsto f(n_2), \dots \}</math>. If a node based function is denoted by a complex expression <math>F</math>, then the evaluation of <math>F</math> at a single node <math>n</math> will be denoted using the pipe-subscript notation: <math>F|_n</math>.

An arbitrary function <math>f: \mathcal{E} \to \R</math> over the set of edges will be referred to as an "'''edge based function'''", and is effectively a scalar or vector field. The notation <math>\langle e \in \mathcal{E} : f(e) \rangle</math> will be used if the input parameter is ambiguous. An edge based function can also be denoted by exhaustively listing the assignment to each edge: <math>\{e_1 \mapsto f(e_1), e_2 \mapsto f(e_2), \dots \}</math>. If an edge based function is denoted by a complex expression <math>F</math>, then the evaluation of <math>F</math> at a single edge <math>e</math> will be denoted using the pipe-subscript notation: <math>F|_e</math>. When a real number assigned to an edge denotes a "vector", a positive value denotes an orientation along the edge's preferred direction, and a negative value denotes an orientation against the edge's preferred direction.

[[File:Example directed graph.svg|thumb|An example directed graph.]]

An example directed graph is shown to the right. The node set is <math>\mathcal{N} = \{n_1, n_2, n_3, n_4, n_5\}</math> and the edge set is <math>\mathcal{E} = \{e_1, e_2, e_3, e_4, e_5, e_6\}</math>. The origin and terminal nodes of each edge are: <math>(t_\bullet(e_1), t^\bullet(e_1)) = (n_1, n_2)</math>; <math>(t_\bullet(e_2), t^\bullet(e_2)) = (n_1, n_3)</math>; <math>(t_\bullet(e_3), t^\bullet(e_3)) = (n_2, n_4)</math>; <math>(t_\bullet(e_4), t^\bullet(e_4)) = (n_4, n_3)</math>; <math>(t_\bullet(e_5), t^\bullet(e_5)) = (n_4, n_5)</math>; and <math>(t_\bullet(e_6), t^\bullet(e_6)) = (n_5, n_1)</math>.

=== Interpreting nodes and edges ===

[[File:Directed graph discrete volumes.svg|thumb|Two irregular volumes with an irregular boundary are characterized as two points linked by a curve. Each point is a node in a directed graph, and the curve denotes the directed edge that links the two points. The "length" of the directed edge is the distance between the two points, and the "cross-sectional area" of the directed edge is the area of the projection of the boundary surface onto a plane that is perpendicular to the displacement between the two points.]]

Each node <math>n \in \mathcal{N}</math> will correspond to both a point <math>\omega_0(n)</math> in space, and a volume <math>\omega_3(n)</math> that contains <math>\omega_0(n)</math>. <math>\tau_3(n)</math> will denote the volume of <math>\omega_3(n)</math>, also referred to as the volume of node <math>n</math>. The volume of each node must be positive: <math>\tau_3(n) > 0</math> for all <math>n \in \mathcal{N}</math>. In the example directed graph, an example choice of volumes for each node is <math>\tau_3 = \{n_1 \mapsto 3, n_2 \mapsto 2, n_3 \mapsto 2, n_4 \mapsto 1, n_5 \mapsto 4\}</math>.

Each edge <math>e \in \mathcal{E}</math> will correspond to both a curve <math>\omega_1(e)</math> originating from <math>\omega_0(t_\bullet(e))</math> and terminating at <math>\omega_0(t^\bullet(e))</math>, and a surface <math>\omega_2(e)</math> that separates <math>\omega_3(t_\bullet(e))</math> from <math>\omega_3(t^\bullet(e))</math> oriented from <math>\omega_3(t_\bullet(e))</math> to <math>\omega_3(t^\bullet(e))</math>.

<math>\tau_1(e)</math> will denote the distance <math>|\omega_0(t^\bullet(e))-\omega_0(t_\bullet(e))|</math> from <math>\omega_0(t_\bullet(e))</math> to <math>\omega_0(t^\bullet(e))</math>, which will be referred to as the length of edge <math>e</math>. The length of each edge must be positive: <math>\tau_1(e) > 0</math> for all <math>e \in \mathcal{E}</math>.
<math>\tau_2(e)</math> will denote the area of <math>\omega_2(e)</math> projected onto a plane that is perpendicular to the displacement <math>\omega_0(t^\bullet(e))-\omega_0(t_\bullet(e))</math>, which will be referred to as the cross-sectional area, thickness, or width of edge <math>e</math>. The cross-sectional area of each edge must be positive: <math>\tau_2(e) > 0</math> for all <math>e \in \mathcal{E}</math>.
In the example directed graph, an example choice of lengths for each edge is <math>\tau_1 = \{e_1 \mapsto 1, e_2 \mapsto 1, e_3 \mapsto 2, e_4 \mapsto 1, e_5 \mapsto 3, e_6 \mapsto 4\}</math>, and an example choice of areas for each edge is <math>\tau_2 = \{e_1 \mapsto 1, e_2 \mapsto 2, e_3 \mapsto 1, e_4 \mapsto 2, e_5 \mapsto 3, e_6 \mapsto 3\}</math>.

<div style="font-weight:bold;">In summary an arbitrary directed graph <math>\mathcal{G}</math> will be characterized by: 
*the set of nodes <math>\mathcal{N}</math>
*the set of edges <math>\mathcal{E}</math> 
*the start <math>t_\bullet(e)</math> and end <math>t^\bullet(e)</math> of each edge <math>e \in \mathcal{E}</math>
*the volume <math>\tau_3(n)</math> of each node <math>n \in \mathcal{N}</math>
*the length <math>\tau_1(e)</math> and area <math>\tau_2(e)</math> of each edge <math>e \in \mathcal{E}</math>
</div>

The point <math>\omega_0(n)</math> and space <math>\omega_3(n)</math> associated with each node <math>n \in \mathcal{N}</math>, and the curve <math>\omega_1(e)</math> and surface <math>\omega_2(e)</math> associated with each edge <math>e \in \mathcal{E}</math>, are not given consideration after the approximate directed graph model has been established. The parameters in the above list are all that is needed to utilize the directed graph model.

=== Points and multi-points ===

A "point" in space is quantified by a single node <math>n_0 \in \mathcal{N}</math>, which corresponds to the space point <math>\omega_0(n_0)</math>.

Point <math>n_0</math> can be denoted by the node based function <math>\delta_0(n;n_0) : \mathcal{N} \to \R</math> where for each node <math>n \in \mathcal{N}</math> it is the case that <math>\delta_0(n;n_0) = \left\{\begin{array}{cc} 1/\tau_3(n) & (n = n_0) \\ 0 & (n \neq n_0) \\ \end{array}\right.</math>. The term <math>\frac{1}{\tau_3(n)}</math> arises from the fact that the point is "spread out" over the volume of <math>n</math>. The notation <math>\delta_0(n_0)</math> denotes <math>\delta_0(n;n_0)</math> with the parameter <math>n</math> not shown.

A superposition of points is a "multi-point". A multi-point is effectively a set of points with varying weights, which can be any real number. Let the set <math>W = \{(n_1,w_1),(n_2,w_2),...,(n_k,w_k)\}</math> of node/weight pairs denote a multi-point. For each <math>i = 1, 2, ..., k</math>: <math>w_i</math> is the weight associated with node <math>n_i</math>. The node based function that denotes multi-point <math>W</math> is: <math>\delta_0(n;W) = \sum_{i=1}^k w_i\delta_0(n;n_i) = \sum_{i=1}^k \left\{\begin{array}{cc} w_i/\tau_3(n) & (n = n_i) \\ 0 & (n \neq n_i) \\ \end{array}\right.</math>. Any node based function that returns real values can denote a multi-point.

'''Any node based function that denotes "density" is best interpreted as a multi-point.'''

Below is shown two examples of a multi-point <math>W</math> and the corresponding node based function <math>\delta_0(n;W)</math>. On the left, there is a collection of points with a weight of +1 denoted by red dots, and points with a weight of -1 denoted by blue dots. For each volume <math>\omega_3(n)</math>, the weights of all of the contained points are added together and the total weight is divided by <math>\tau_3(n)</math> to spread it out over the volume <math>\omega_3(n)</math>. On the right, shades of red denote positive values for <math>\delta_0(n;W)</math>, while shades of blue denote negative values for <math>\delta_0(n;W)</math>. White means that <math>\delta_0(n;W) = 0</math> for the current node <math>n</math>.

{{Gallery
|height=225
|width=450
|File:Multi point function 1.svg|Multi point example on the left, and function on the right.
|File:Multi point function 2.svg|Multi point example on the left, and function on the right.
}}

=== Volumes and multi-volumes ===

[[File:Example directed graph with volumes.svg|thumb|The example directed graph with the example multi-volume.]]

The discrete analog to volumes using a directed graph is that a volume <math>U</math> is a collection of nodes: <math>U \subseteq \mathcal{N}</math>. The combined volume region denoted by <math>U</math> is <math>\omega_3(U) = \bigcup_{n \in U} \omega_3(n)</math> and the total volume amount is <math>\tau_3(U) = \sum_{n \in U} \tau_3(n)</math>.

Volume <math>U</math> can be denoted by a node based membership function: <math>\delta_3(n;U) : \mathcal{N} \to \R</math> where for each node <math>n \in \mathcal{N}</math>, it is the case that <math>\delta_3(n;U) = \left\{\begin{array}{cc} 1 & (n \in U) \\ 0 & (n \notin U) \\ \end{array}\right.</math>. The notation <math>\delta_3(U)</math> denotes <math>\delta_3(n;U)</math> with the parameter <math>n</math> not shown.

A superposition of volumes forms a "multi-volume". Any node based function that returns real values can be the membership function of a multi-volume. For example, given the node based function <math>f = \{n_1 \mapsto -2, n_2 \mapsto 0, n_3 \mapsto 4, n_4 \mapsto 0.5, n_5 \mapsto 0.5\}</math> can be expressed as the following superpositions:

<math>f(n) = -2\delta_3(n;\{n_1\}) + 4\delta_3(n;\{n_3\}) + 0.5\delta_3(n;\{n_4,n_5\})</math> is the superposition of volumes <math>\{n_1\}, \{n_3\}, \{n_4, n_5\}</math> with the respective weights <math>-2, 4, 0.5</math>.

<math>f(n) = -2\delta_3(n;\{n_1\}) + 0.5\delta_3(n;\{n_3,n_4,n_5\}) + 3.5\delta_3(n;\{n_3\})</math> is the superposition of volumes <math>\{n_1\}, \{n_3,n_4,n_5\}, \{n_3\}</math> with the respective weights <math>-2, 0.5, 3.5</math>.

It is important to note that given an arbitrary node based function, decomposing the multi-volume denoted by the function into a superposition of volumes is '''not unique''', as shown above where a node based function is decomposed into multiple superpositions of volumes.

'''Any node based function that denotes a "potential" is best interpreted as a multi-volume.'''

=== Paths and multi-paths ===

[[File:Example directed graph with paths.svg|thumb|The example directed graph with the example multi-path.]]

The discrete analog to paths using a directed graph is that a path <math>P</math> is a series of edges and orientations: <math>P = \langle (e_1,s_1), (e_2,s_2), \dots, (e_k,s_k) \rangle</math> where <math>e_i \in \mathcal{E}</math> is the <math>i^\text{th}</math> edge and <math>s_i \in \{+1, -1\}</math> is an indicator that is <math>+1</math> if <math>e_i</math> is traversed in the preferred direction, and is <math>-1</math> if <math>e_i</math> is traversed in the opposite direction. The path must be unbroken: for each <math>i \in \{1, 2, \dots, k-1\}</math> it must be that <math>\left\{\begin{array}{cc} t^\bullet(e_i) & (s_i = +1) \\ t_\bullet(e_i) & (s_i = -1) \\ \end{array}\right. = \left\{\begin{array}{cc} t_\bullet(e_{i+1}) & (s_{i+1} = +1) \\ t^\bullet(e_{i+1}) & (s_{i+1} = -1) \\ \end{array}\right.</math>. The combined path denoted by <math>P</math> is <math>\omega_1(P) = \bigcup_{i = 1}^k s_i\omega_1(e_i)</math> (path <math>-\omega_1(e_i)</math> is path <math>\omega_1(e_i)</math> the orientation reversed) and the total length is <math>\tau_1(P) = \sum_{i=1}^k \tau_1(e_i)</math>.

A path <math>P</math> can be denoted by the edge based function <math>\delta_1(e;P) : \mathcal{E} \to \R</math> where for each edge <math>e \in \mathcal{E}</math>, it is the case that <math>\delta_1(e;P) = \sum_{i=1}^k s_i\delta_1(e;e_i)</math> where <math>\delta_1(e;e_i) = \left\{\begin{array}{cc} 1/\tau_2(e) & (e = e_i) \\ 0 & (e \neq e_i) \\ \end{array}\right.</math>. The term <math>\frac{1}{\tau_2(e)}</math> arises from the fact that the path is "spread out" over the thickness of <math>e</math>. The notation <math>\delta_1(P)</math> denotes <math>\delta_1(e;P)</math> with the parameter <math>e</math> not shown.

A superposition of paths forms a "multi-path". Any edge based function that returns real values can denote a multi-path. For example, assuming that <math>\tau_2(e) = 1</math> for all <math>e \in \mathcal{E}</math>, the edge based function <math>f = \{e_1 \mapsto 1, e_2 \mapsto -1, e_3 \mapsto 1, e_4 \mapsto 2, e_5 \mapsto 0, e_6 \mapsto 0\}</math> can be expressed as the following superpositions:

<math>f(e) = \delta_1(e;\langle (e_1,+1) \rangle) - \delta_1(e;\langle (e_2,+1) \rangle) + \delta_1(e;\langle (e_3,+1) \rangle) + 2\delta_1(e;\langle (e_4,+1) \rangle)</math> is the superposition of paths <math>\langle (e_1,+1) \rangle, \langle (e_2,+1) \rangle, \langle (e_3,+1) \rangle</math>, and <math>\langle (e_4,+1) \rangle</math> with the respective weights <math>1, -1, 1, 2</math>.

<math>f(e) = \delta_1(e;\langle (e_2,-1), (e_1,+1), (e_3,+1), (e_4,+1) \rangle) + \delta_1(e;\langle (e_4,+1) \rangle)</math> is the superposition of paths <math>\langle (e_2,-1), (e_1,+1), (e_3,+1), (e_4,+1) \rangle</math>, and <math>\langle (e_4,+1) \rangle</math> with the respective weights <math>1, 1</math>.

It is important to note that given an arbitrary edge based function, decomposing the multi-path denoted by the function into a superposition of paths is '''not unique''', as shown above where an edge based function is decomposed into multiple superpositions of paths.

'''Any edge based function that denotes "flow density" (flow per unit cross-sectional area) is best interpreted as a multi-path.'''

Below is shown two examples of a multi-path <math>P</math> and the corresponding edge based function <math>\delta_1(e;P)</math>. On the left a superposition of multiple paths is shown. For each surface <math>\omega_2(e)</math>, the total "flow" through <math>\omega_2(e)</math> is computed and is then divided by <math>\tau_2(e)</math> to spread it out over the surface <math>\omega_2(e)</math>. On the right, the flow enters a surface through the green side and leaves through the red side. The density of the flow is depicted by the strength of the shading. If the green to red orientation coincides with the orientation of <math>e</math>, then <math>\delta_1(e;P)</math> is positive. If the green to red orientation is opposite the orientation of <math>e</math>, then <math>\delta_1(e;P)</math> is negative.

{{Gallery
|height=250
|width=500
|File:Multi path function 1.svg|Multi path example on the left, and edge based function on the right.
|File:Multi path function 2.svg|Multi path example on the left, and edge based function on the right.
}}

=== Surfaces and multi-surfaces ===

[[File:Example directed graph with surfaces.svg|thumb|The example directed graph with the example multi-surface.]]

The discrete analog to oriented surfaces using a directed graph is that an oriented surface <math>S</math> is a collection of edges and orientations: <math>S = \{(e_1,s_1),(e_2,s_2),\dots,(e_k,s_k)\}</math> where all of the edges <math>e_1, e_2, \dots, e_k \in \mathcal{E}</math> are distinct and <math>s_i \in \{+1, -1\}</math> is an indicator that describes the surface's orientation. Each edge <math>e_i</math> is an edge that passes through the surface, and <math>s_i</math> is <math>+1</math> if <math>e_i</math> passes through the surface in the preferred direction, and <math>s_i</math> is <math>-1</math> if <math>e_i</math> passes through the surface in the opposite direction. The combined surface denoted by <math>S</math> is <math>\omega_2(S) = \bigcup_{i = 1}^k s_i\omega_2(e_i)</math> (surface <math>-\omega_2(e_i)</math> is surface <math>\omega_2(e_i)</math> with the orientation reversed) and the total area is <math>\tau_2(S) = \sum_{i=1}^k \tau_2(e_i)</math>.

A surface <math>S</math> can be denoted by the edge based function <math>\delta_2(e;S) : \mathcal{E} \to \R</math> where for each edge <math>e \in \mathcal{E}</math>, it is the case that <math>\delta_2(e;S) = \left\{\begin{array}{cc} +1/\tau_1(e) & ((e,+1) \in S) \\ -1/\tau_1(e) & ((e,-1) \in S) \\ 0 & (\text{else}) \\ \end{array}\right.</math>. The term <math>\frac{1}{\tau_1(e)}</math> arises from the fact that the surface is "spread out" over the length of <math>e</math>. The notation <math>\delta_2(S)</math> denotes <math>\delta_2(e;S)</math> with the parameter <math>e</math> not shown.

A superposition of surfaces forms a "multi-surface". Any edge based function that returns real values can denote a multi-surface. For example, assuming that <math>\tau_1(e) = 1</math> for all <math>e \in \mathcal{E}</math>, the edge based function <math>f = \{e_1 \mapsto 2, e_2 \mapsto 2, e_3 \mapsto 0, e_4 \mapsto 0, e_5 \mapsto 0, e_6 \mapsto -1\}</math> can be expressed as the following superpositions:

<math>f(e) = 2\delta_2(e;\{(e_1,+1)\}) + 2\delta_2(e;\{(e_2,+1)\}) - \delta_2(e;\{(e_6,+1)\})</math> is the superposition of surfaces <math>\{(e_1,+1)\}, \{(e_2,+1)\},</math> and <math>\{(e_6,+1)\}</math> with the respective weights <math>2, 2, -1</math>.

<math>f(e) = \delta_2(e;\{(e_1,+1),(e_2,+1),(e_6,-1)\}) + \delta_2(e;\{(e_1,+1),(e_2,+1)\})</math> is the superposition of surfaces <math>\{(e_1,+1),(e_2,+1),(e_6,-1)\}</math> and <math>\{(e_1,+1),(e_2,+1)\}</math> with the respective weights <math>1</math> and <math>1</math>.

It is important to note that given an arbitrary edge based function, decomposing the multi-surface denoted by the function into a superposition of surfaces is '''not unique''', as shown above where an edge based function is decomposed into multiple superpositions of surfaces.

'''Any edge based function that denotes the "rate of gain" (rate of gain per unit length) is best interpreted as a multi-surface.'''

Below is shown two examples of a multi-surface <math>S</math> and the corresponding edge based function <math>\delta_2(e;S)</math>. On the left a superposition of multiple surfaces is shown. For each curve <math>\omega_1(e)</math>, the net number of times edge <math>e</math> penetrates <math>S</math> in the preferred direction, referred to as the "gain" across <math>\omega_1(e)</math>, is computed and is then divided by <math>\tau_1(e)</math> to spread the gain out over the curve <math>\omega_1(e)</math>. On the right, the net gain is positive if an edge is traversed from the green end to the red end. The amount of gain is depicted by the strength of the shading. If the green to red orientation coincides with the orientation of <math>e</math>, then <math>\delta_2(e;S)</math> is positive. If the green to red orientation is opposite the orientation of <math>e</math>, then <math>\delta_2(e;S)</math> is negative.

{{Gallery
|height=250
|width=500
|File:Multi surface function 1.svg|Multi surface example on the left, and edge based function on the right.
|File:Multi surface function 2.svg|Multi surface example on the left, and edge based function on the right.
}}

=== Intersections ===

==== Multi-point multi-volume intersections ====

[[File:Point volume intersections.svg|thumb|A multi-volume and a multi-point is shown. The red volumes and points have a weight of +1, and the cyan volumes and points have a weight of -1. The intersection contributed by a point is its weight times the weight assigned to the volume that contains the point. The total intersection weight in this image is: N = (+1)x(+1+1+1-1-1) + 0x(+1+1-1-1-1) + (-1)x(+1-1) = +1.]]

Given an arbitrary point <math>n_0 \in \mathcal{N}</math> and volume <math>U \subseteq \mathcal{N}</math>, let <math>\Sigma_\mathcal{N}(n_0,U) = \left\{\begin{array}{cc} 1 & (n_0 \in U) \\ 0 & (n_0 \notin U) \end{array}\right.</math> denote an indicator function that returns 1 if point <math>n_0</math> is contained by volume <math>U</math> and returns 0 if otherwise.

Node based function <math>\delta_0(n;n_0) = \left\{\begin{array}{cc} 1/\tau_3(n) & (n = n_0) \\ 0 & (n \neq n_0) \end{array}\right.</math> denotes the point <math>n_0</math>, and node based function <math>\delta_3(n;U) = \left\{\begin{array}{cc} 1 & (n \in U) \\ 0 & (n \notin U) \end{array}\right.</math> denotes the volume <math>U</math>.

The indicator function <math>\Sigma_\mathcal{N}(n_0,U)</math> can be computed via the following sum: 
<math>\Sigma_\mathcal{N}(n_0,U) = \sum_{n \in \mathcal{N}} \left\{\begin{array}{cc} 1 & (n = n_0) \\ 0 & (n \neq n_0) \end{array}\right\}\cdot\left\{\begin{array}{cc} 1 & (n \in U) \\ 0 & (n \notin U) \end{array}\right\} = \sum_{n \in \mathcal{N}} \delta_0(n;n_0)\delta_3(n;U)\tau_3(n)</math>. 
The term <math>\tau_3(n)</math> cancels out the <math>\frac{1}{\tau_3(n)}</math> in <math>\delta_0(n;n_0)</math>.

<math>\Sigma_\mathcal{N}(n_0,U)</math> counts the number of times point <math>n_0</math> occurs in volume <math>U</math> which can only be 0 or 1 for simple points and volumes. This count can be generalized to sets of points. Let <math>W \subseteq \mathcal{N}</math> denote a collection of points. The collection of points can be denoted by the node based function <math>\delta_0(n;W) = \left\{\begin{array}{cc} 1/\tau_3(n) & (n \in W) \\ 0 & (n \notin W) \end{array}\right.</math>. The number of points from <math>W</math> that are contained by volume <math>U</math> can again be computed by the sum: 
<math>\Sigma_\mathcal{N}(W,U) = \sum_{n \in \mathcal{N}} \left\{\begin{array}{cc} 1 & (n \in W) \\ 0 & (n \notin W) \end{array}\right\}\cdot\left\{\begin{array}{cc} 1 & (n \in U) \\ 0 & (n \notin U) \end{array}\right\} = \sum_{n \in \mathcal{N}} \delta_0(n;W)\delta_3(n;U)\tau_3(n)</math>.

In the general case, given a node based function <math>f</math> that denotes a multi-point, and a node based function <math>g</math> that denotes a multi-volume, then the total weight of all the occurrences of <math>f</math> being inside of <math>g</math> is <math>\Sigma_\mathcal{N}(f,g) = \sum_{n \in \mathcal{N}} f(n)g(n)\tau_3(n)</math>. The quantity <math>\Sigma_\mathcal{N}(f,g)</math> is called the '''net intersection''' or the '''total intersection weight''' of <math>f</math> with <math>g</math>.

It should also be noted that <math>\Sigma_\mathcal{N}(f,g)</math> is commutative: <math>\Sigma_\mathcal{N}(f,g) = \Sigma_\mathcal{N}(g,f)</math>

==== Multi-path multi-surface intersections ====

[[File:Path surface intersections.svg|thumb|A path is shown intersecting a surface multiple times. The boundary of the surface has a counter-clockwise orientation. The net number of times that the surface is crossed in the preferred direction is +1-1+1+1+1 = 3.]]

Given an arbitrary path <math>P = \langle (e_{P,1},s_{P,1}), (e_{P,2},s_{P,2}), ..., (e_{P,k_P},s_{P,k_P}) \rangle</math> and surface <math>S = \{ (e_{S,1},s_{S,1}), (e_{S,2},s_{S,2}), ..., (e_{S,k_S},s_{S,k_S}) \}</math>, let <math>\Sigma_\mathcal{E}(P,S)</math> denote the number of times path <math>P</math> passes through surface <math>S</math> in the preferred direction minus the number of times path <math>P</math> passes through surface <math>S</math> in the opposite direction. For an arbitrary <math>i = 1, 2, ..., k_P</math>, and an arbitrary <math>j = 1, 2, ..., k_S</math>, if <math>e_{P,i} = e_{S,j}</math> then the <math>i^\text{th}</math> link of path <math>P</math> coincides with the <math>j^\text{th}</math> segment of surface <math>S</math>. If <math>s_{P,i}s_{S,j} = +1</math>, then <math>s_{P,i} = s_{S,j}</math> and <math>P</math> passes through <math>S</math> in the preferred direction. If <math>s_{P,i}s_{S,j} = -1</math>, then <math>s_{P,i} \neq s_{S,j}</math> and <math>P</math> passes through <math>S</math> in the opposite direction.

<math>\Sigma_\mathcal{E}(P,S)</math> can be computed via the following sum: 
<math>\Sigma_\mathcal{E}(P,S) = \sum_{i=1}^{k_P}\sum_{j=1}^{k_S} \left\{\begin{array}{cc} s_{P,i}s_{S,j} & (e_{P,i} = e_{S,j}) \\ 0 & (e_{P,i} \neq e_{S,j}) \end{array}\right. = \sum_{e \in \mathcal{E}} \delta_1(e;P)\delta_2(e;S)\tau_2(e)\tau_1(e)</math>. The term <math>\tau_2(e)</math> cancels out the <math>\frac{1}{\tau_2(e)}</math> from <math>\delta_1(e;P)</math>, and the term <math>\tau_1(e)</math> cancels out the <math>\frac{1}{\tau_1(e)}</math> from <math>\delta_2(e;S)</math>.

In the general case, given an edge based function <math>f</math> that denotes a multi-path, and an edge based function <math>g</math> that denotes a multi-surface, then the total weight of all the intersections of <math>f</math> with <math>g</math> is <math>\Sigma_\mathcal{E}(f,g) = \sum_{e \in \mathcal{E}} f(e)g(e)\tau_2(e)\tau_1(e)</math>. The quantity <math>\Sigma_\mathcal{E}(f,g)</math> is called the '''net intersection''' or the '''total intersection weight''' of <math>f</math> with <math>g</math>.

It should also be noted that <math>\Sigma_\mathcal{E}(f,g)</math> is commutative: <math>\Sigma_\mathcal{E}(f,g) = \Sigma_\mathcal{E}(g,f)</math>

== Lumped coordinate system models ==

=== Model of Cartesian coordinates ===

[[File:Lumped cartesian coordinate model.svg|thumb|A directed graph that models the Cartesian coordinate system. One dimension is suppressed for simplicity.]]

The discrete model used for Cartesian coordinates will consist of an infinite 3 dimensional grid of nodes that spans <math>\R^3</math>. To establish the model, the following is needed:
* A spacing <math>\Delta x</math> for the x-coordinates.
* A spacing <math>\Delta y</math> for the y-coordinates.
* A spacing <math>\Delta z</math> for the z-coordinates.

For each triplet of integer indices <math>(i,j,k) \in \Z^3</math>, there exists a node <math>n_{(i,j,k)}</math> and 3 edges <math>e_{x,(i,j,k)}</math>, <math>e_{y,(i,j,k)}</math>, and <math>e_{z,(i,j,k)}</math>.

The origin and terminal nodes of <math>e_{x,(i,j,k)}</math> are <math>t_\bullet(e_{x,(i,j,k)}) = n_{(i,j,k)}</math> and <math>t^\bullet(e_{x,(i,j,k)}) = n_{(i+1,j,k)}</math>

The origin and terminal nodes of <math>e_{y,(i,j,k)}</math> are <math>t_\bullet(e_{y,(i,j,k)}) = n_{(i,j,k)}</math> and <math>t^\bullet(e_{y,(i,j,k)}) = n_{(i,j+1,k)}</math>

The origin and terminal nodes of <math>e_{z,(i,j,k)}</math> are <math>t_\bullet(e_{z,(i,j,k)}) = n_{(i,j,k)}</math> and <math>t^\bullet(e_{z,(i,j,k)}) = n_{(i,j,k+1)}</math>

The Cartesian coordinate associated with node <math>n_{(i,j,k)}</math> is <math>\omega_0(n_{(i,j,k)}) = (i\Delta x, j\Delta y, k\Delta z)</math>, and the space associated with node <math>n_{(i,j,k)}</math> is <math>\omega_3(n_{(i,j,k)}) = \left[i\Delta x - \frac{\Delta x}{2}, i\Delta x + \frac{\Delta x}{2}\right] \times \left[j\Delta y - \frac{\Delta y}{2}, j\Delta y + \frac{\Delta y}{2}\right] \times \left[k\Delta z - \frac{\Delta z}{2}, k\Delta z + \frac{\Delta z}{2}\right]</math>.

The volume of node <math>n_{(i,j,k)}</math> is <math>\tau_3(n_{(i,j,k)}) = \Delta x \Delta y \Delta z</math>.

The length of edge <math>e_{x,(i,j,k)}</math> is <math>\tau_1(e_{x,(i,j,k)}) = \Delta x</math>. 
The length of edge <math>e_{y,(i,j,k)}</math> is <math>\tau_1(e_{y,(i,j,k)}) = \Delta y</math>.
The length of edge <math>e_{z,(i,j,k)}</math> is <math>\tau_1(e_{z,(i,j,k)}) = \Delta z</math>.

The area of edge <math>e_{x,(i,j,k)}</math> is <math>\tau_2(e_{x,(i,j,k)}) = \Delta y\Delta z</math>. 
The area of edge <math>e_{y,(i,j,k)}</math> is <math>\tau_2(e_{y,(i,j,k)}) = \Delta z\Delta x</math>.
The area of edge <math>e_{z,(i,j,k)}</math> is <math>\tau_2(e_{z,(i,j,k)}) = \Delta x\Delta y</math>.

Given a scalar field <math>f(x,y,z)</math> over Cartesian coordinates, <math>f</math> will be approximated by the node based function <math>f_\text{approx} : \mathcal{N} \to \R</math>. Function <math>f_\text{approx}</math> is defined by <math>f_\text{approx}(n_{(i,j,k)}) = f(i\Delta x, j\Delta y, k\Delta z)</math> for each <math>(i,j,k) \in \Z^3</math>.

Given a vector field <math>\mathbf{F}(x,y,z) = F_x(x,y,z)\mathbf{i} + F_y(x,y,z)\mathbf{j} + F_z(x,y,z)\mathbf{k}</math> over Cartesian coordinates, <math>\mathbf{F}</math> will be approximated by the edge based function <math>F_\text{approx} : \mathcal{E} \to \R</math>. Function <math>F_\text{approx}</math> is defined by <math>F_\text{approx}(e_{x,(i,j,k)}) = F_x(i\Delta x, j\Delta y, k\Delta z)</math>; <math>F_\text{approx}(e_{y,(i,j,k)}) = F_y(i\Delta x, j\Delta y, k\Delta z)</math>; and <math>F_\text{approx}(e_{z,(i,j,k)}) = F_z(i\Delta x, j\Delta y, k\Delta z)</math>.

=== Model of Cylindrical coordinates ===

[[File:Lumped cylindrical coordinate model.svg|thumb|A directed graph that models the cylindrical coordinate system. For simplicity, the vertical/z dimension is not shown.]]

The discrete model used for Cylindrical coordinates will consist of an infinite 3 dimensional cylindrical grid of nodes that spans <math>\R^3</math>. To establish the model, the following is needed:

*A spacing <math>\Delta \rho</math> for the coordinate <math>\rho</math>.
*A large positive integer <math>N_\phi</math> that will yield a spacing <math>\Delta \phi = \frac{2\pi}{N_\phi}</math> for the coordinate <math>\phi</math>.
*A spacing <math>\Delta z</math> for the coordinate <math>z</math>.

Each node <math>n \in \mathcal{N}</math> will correspond to a unique triplet <math>(i,j,k)</math> of integers, and will be denoted by <math>n_{(i,j,k)}</math>. The cylindrical coordinate denoted by <math>n_{(i,j,k)}</math> will be <math>\omega_0(n_{(i,j,k)}) = (i\Delta\rho, j\Delta\phi, k\Delta z)</math>.

Clearly, due to the nature of cylindrical coordinates (such as the fact that coordinate <math>\phi</math> is cyclical), not every triplet <math>(i,j,k)</math> will be allowed. The following three domains will be defined: <math>I_0</math>, <math>I_1</math>, and <math>I_2</math>.

<math>(i,j,k) \in I_0 \iff (i \geq 1 \;\text{and}\; 0 \leq j < N_\phi)</math>

<math>(i,j,k) \in I_1 \iff (i \geq 0 \;\text{and}\; (i = 0 \implies j = 0) \;\text{and}\; 0 \leq j < N_\phi)</math>

<math>(i,j,k) \in I_2 \iff (i \geq 0 \;\text{and}\; 0 \leq j < N_\phi)</math>

There exists a node <math>n_{(i,j,k)}</math> for each <math>(i,j,k) \in I_1</math>.

There exists an edge <math>e_{\rho,(i,j,k)}</math> for each <math>(i,j,k) \in I_2</math> such that <math>t_\bullet(e_{\rho,(i,j,k)}) = \left\{\begin{array}{cc} n_{(i,j,k)} & (i > 0) \\ n_{(i,0,k)} & (i = 0) \\ \end{array}\right.</math> and <math>t^\bullet(e_{\rho,(i,j,k)}) = n_{(i+1,j,k)}</math>

There exists an edge <math>e_{\phi,(i,j,k)}</math> for each <math>(i,j,k) \in I_0</math> such that <math>t_\bullet(e_{\phi,(i,j,k)}) = n_{(i,j,k)}</math> and <math>t^\bullet(e_{\phi,(i,j,k)}) = \left\{\begin{array}{cc} n_{(i,j+1,k)} & (j < N_\phi-1) \\ n_{(i,0,k)} & (j = N_\phi-1) \\ \end{array}\right.</math>

There exists an edge <math>e_{z,(i,j,k)}</math> for each <math>(i,j,k) \in I_1</math> such that <math>t_\bullet(e_{z,(i,j,k)}) = n_{(i,j,k)}</math> and <math>t^\bullet(e_{z,(i,j,k)}) = n_{(i,j,k+1)}</math>

Next, the volume of each node and the lengths and area of each edge will be calculated. To simplify notation, let <math>l_\phi(i_\text{real}) = (i_\text{real}\Delta \rho)\Delta\phi</math> for any nonnegative real number <math>i_\text{real}</math>.

For each <math>(i,j,k) \in I_1</math> the cylindrical coordinate that corresponds to <math>n_{(i,j,k)}</math> is <math>\omega_0(n_{(i,j,k)}) = (i\Delta\rho, j\Delta\phi, k\Delta z)</math>, and the space (described using cylindrical coordinates) associated with <math>n_{(i,j,k)}</math> is <math>\omega_3(n_{(i,j,k)}) = \left\{\begin{array}{cc} [(i-\frac{1}{2})\Delta\rho, (i+\frac{1}{2})\Delta\rho] \times [(j-\frac{1}{2})\Delta\phi, (j+\frac{1}{2})\Delta\phi] \times  [(k-\frac{1}{2})\Delta z, (k+\frac{1}{2})\Delta z] & (i > 0) \\ 
{[0, \frac{\Delta\rho}{2}]} \times [0, 2\pi] \times [(k-\frac{1}{2})\Delta z, (k+\frac{1}{2})\Delta z] & (i = 0) \\ \end{array}\right.</math>

For each <math>(i,j,k) \in I_1</math> the volume of node <math>n_{(i,j,k)}</math> is <math>\tau_3(n_{(i,j,k)}) = \left\{\begin{array}{cc} \Delta\rho \cdot l_\phi(i) \cdot \Delta z & (i > 0) \\ \pi(\frac{\Delta\rho}{2})^2\Delta z & (i = 0) \\ \end{array}\right.</math>

For each <math>(i,j,k) \in I_2</math> the length and area of edge <math>e_{\rho,(i,j,k)}</math> is <math>\tau_1(e_{\rho,(i,j,k)}) = \Delta\rho</math> and <math>\tau_2(e_{\rho,(i,j,k)}) \approx l_\phi(i+\frac{1}{2}) \cdot \Delta z</math> respectively.

For each <math>(i,j,k) \in I_0</math> the length and area of edge <math>e_{\phi,(i,j,k)}</math> is <math>\tau_1(e_{\phi,(i,j,k)}) \approx l_\phi(i)</math> and <math>\tau_2(e_{\phi,(i,j,k)}) = \Delta\rho \cdot \Delta z</math> respectively.

For each <math>(i,j,k) \in I_1</math> the length and area of edge <math>e_{z,(i,j,k)}</math> is <math>\tau_1(e_{z,(i,j,k)}) = \Delta z</math> and <math>\tau_2(e_{z,(i,j,k)}) = \left\{\begin{array}{cc} \Delta\rho \cdot l_\phi(i) & (i > 0) \\ \pi(\frac{\Delta\rho}{2})^2 & (i = 0) \\ \end{array}\right.</math> respectively.

Given a scalar field <math>f(\rho,\phi,z)</math> over cylindrical coordinates, <math>f</math> will be approximated by the node based function <math>f_\text{approx} : \mathcal{N} \to \R</math>. Function <math>f_\text{approx}</math> is defined by <math>f_\text{approx}(n_{(i,j,k)}) = f(i\Delta\rho, j\Delta\phi, k\Delta z)</math> for each <math>(i,j,k) \in I_1</math>.

Given a vector field <math>\mathbf{F}(\rho,\phi,z) = F_\rho(\rho,\phi,z)\hat{\mathbf{\rho}} + F_\phi(\rho,\phi,z)\hat{\mathbf{\phi}} + F_z(\rho,\phi,z)\hat{\mathbf{z}}</math> over cylindrical coordinates, <math>\mathbf{F}</math> will be approximated by the edge based function <math>F_\text{approx} : \mathcal{E} \to \R</math>. Function <math>F_\text{approx}</math> is defined by:

<math>\forall (i,j,k) \in I_2 : F_\text{approx}(e_{\rho,(i,j,k)}) = \left\{\begin{array}{cc} F_\rho(i\Delta\rho, j\Delta\phi, k\Delta z) & (i > 0) \\ \lim_{\rho \to 0^+} F_\rho(\rho, j\Delta\phi, k\Delta z) & (i = 0) \\ \end{array}\right.</math>

<math>\forall (i,j,k) \in I_0 : F_\text{approx}(e_{\phi,(i,j,k)}) = F_\phi(i\Delta\rho, j\Delta\phi, k\Delta z)</math>

<math>\forall (i,j,k) \in I_1 : F_\text{approx}(e_{z,(i,j,k)}) = F_z(i\Delta\rho, j\Delta\phi, k\Delta z)</math>

=== Model of Spherical coordinates ===

[[File:Lumped spherical coordinate model.svg|thumb|A directed graph that models the spherical coordinate system. The image is a cross-section that contains the vertical/z axis, which is "up" in the image.]]

The discrete model used for Spherical coordinates will consist of an infinite 3 dimensional spherical grid of nodes that spans <math>\R^3</math>. To establish the model, the following is needed:

*A spacing <math>\Delta r</math> for the coordinate <math>r</math>.
*A large positive integer <math>N_\theta</math> that will yield a spacing <math>\Delta \theta = \frac{\pi}{N_\theta}</math> for the coordinate <math>\theta</math>.
*A large positive integer <math>N_\phi</math> that will yield a spacing <math>\Delta \phi = \frac{2\pi}{N_\phi}</math> for the coordinate <math>\phi</math>.

Each node <math>n \in \mathcal{N}</math> will correspond to a unique triplet <math>(i,j,k)</math> of integers, and will be denoted by <math>n_{(i,j,k)}</math>. The spherical coordinate denoted by <math>n_{(i,j,k)}</math> will be <math>\omega_0(n_{(i,j,k)}) = (i\Delta r, j\Delta\theta, k\Delta\phi)</math>.

Clearly, due to the nature of spherical coordinates (such as the fact that coordinate <math>\phi</math> is cyclical), not every triplet <math>(i,j,k)</math> will be allowed. The following three domains will be defined: <math>I_n</math>, <math>I_r</math>, <math>I_\theta</math>, and <math>I_\phi</math>.

<math>(i,j,k) \in I_n \iff (i \geq 0 \;\text{and}\; 0 \leq j \leq N_\theta \;\text{and}\; 0 \leq k < N_\phi \;\text{and}\; (i = 0 \implies j = 0) \;\text{and}\; ((j = 0 \;\text{or}\; j = N_\theta) \implies k = 0))</math>

<math>(i,j,k) \in I_r \iff (i \geq 0 \;\text{and}\; 0 \leq j \leq N_\theta \;\text{and}\; 0 \leq k < N_\phi \;\text{and}\; ((j = 0 \;\text{or}\; j = N_\theta) \implies k = 0))</math>

<math>(i,j,k) \in I_\theta \iff (i \geq 1 \;\text{and}\; 0 \leq j < N_\theta \;\text{and}\; 0 \leq k < N_\phi)</math>

<math>(i,j,k) \in I_\phi \iff (i \geq 1 \;\text{and}\; 0 < j < N_\theta \;\text{and}\; 0 \leq k < N_\phi)</math>

There exists a node <math>n_{(i,j,k)}</math> for each <math>(i,j,k) \in I_n</math>.

There exists an edge <math>e_{r,(i,j,k)}</math> for each <math>(i,j,k) \in I_r</math> such that <math>t_\bullet(e_{r,(i,j,k)}) = \left\{\begin{array}{cc} n_{(i,j,k)} & (i > 0) \\ n_{(i,0,0)} & (i = 0) \\ \end{array}\right.</math> and <math>t^\bullet(e_{r,(i,j,k)}) = n_{(i+1,j,k)}</math>

There exists an edge <math>e_{\theta,(i,j,k)}</math> for each <math>(i,j,k) \in I_\theta</math> such that <math>t_\bullet(e_{\theta,(i,j,k)}) = \left\{\begin{array}{cc} n_{(i,j,k)} & (j > 0) \\ n_{(i,j,0)} & (j = 0) \\ \end{array}\right.</math> and <math>t^\bullet(e_{\theta,(i,j,k)}) = \left\{\begin{array}{cc} n_{(i,j+1,k)} & (j < N_\theta-1) \\ n_{(i,j+1,0)} & (j = N_\theta-1) \\ \end{array}\right.</math>

There exists an edge <math>e_{\phi,(i,j,k)}</math> for each <math>(i,j,k) \in I_\phi</math> such that <math>t_\bullet(e_{\phi,(i,j,k)}) = n_{(i,j,k)}</math> and <math>t^\bullet(e_{\phi,(i,j,k)}) = \left\{\begin{array}{cc} n_{(i,j,k+1)} & (k < N_\phi-1) \\ n_{(i,j,0)} & (k = N_\phi-1) \\ \end{array}\right.</math>

Next, the volume of each node and the lengths and area of each edge will be calculated. To simplify notation, let <math>l_\theta(i_\text{real}) = (i_\text{real}\Delta r)\Delta\theta</math> for any nonnegative real number <math>i_\text{real}</math>, and let <math>l_\phi(i_\text{real},j_\text{real}) = (i_\text{real}\Delta r)\sin(j_\text{real}\Delta\theta)\Delta\phi</math> for any nonnegative real number <math>i_\text{real}</math> and any real number <math>j_\text{real} \in [0, N_\theta]</math>.

For each <math>(i,j,k) \in I_n</math> the spherical coordinate that corresponds to <math>n_{(i,j,k)}</math> is <math>\omega_0(n_{(i,j,k)}) = (i\Delta r, j\Delta\theta, k\Delta\phi)</math>, and the space (described using spherical coordinates) associated with <math>n_{(i,j,k)}</math> is <math>\omega_3(n_{(i,j,k)}) = \left\{\begin{array}{cc} [(i-\frac{1}{2})\Delta r, (i+\frac{1}{2})\Delta r] \times [(j-\frac{1}{2})\Delta\theta, (j+\frac{1}{2})\Delta\theta] \times  [(k-\frac{1}{2})\Delta \phi, (k+\frac{1}{2})\Delta \phi] & (i > 0 \;\text{and}\; 0 < j < N_\theta) \\ 
{[(i-\frac{1}{2})\Delta r, (i+\frac{1}{2})\Delta r]} \times [0, \frac{\Delta\theta}{2}] \times [0, 2\pi] & (i > 0 \;\text{and}\; j = 0) \\
{[(i-\frac{1}{2})\Delta r, (i+\frac{1}{2})\Delta r]} \times [\pi-\frac{\Delta\theta}{2},\pi] \times [0, 2\pi] & (i > 0 \;\text{and}\; j = N_\theta) \\
{[0, \frac{\Delta r}{2}]} \times [0, \pi] \times [0, 2\pi] & (i = 0) \\
\end{array}\right.</math>

For each <math>(i,j,k) \in I_n</math> the volume of node <math>n_{(i,j,k)}</math> is <math>\tau_3(n_{(i,j,k)}) \approx \left\{\begin{array}{cc} \Delta r \cdot l_\theta(i) \cdot l_\phi(i,j) & (i > 0 \;\text{and}\; 0 < j < N_\theta) \\ 
\Delta r \cdot \pi(\frac{1}{2}l_\theta(i))^2 & (i > 0 \;\text{and}\; (j = 0 \;\text{or}\; j = N_\theta)) \\ 
\frac{4}{3}\pi(\frac{\Delta r}{2})^3 & (i = 0) \\ \end{array}\right.</math>

For each <math>(i,j,k) \in I_r</math> the length and area of edge <math>e_{r,(i,j,k)}</math> is <math>\tau_1(e_{r,(i,j,k)}) = \Delta r</math> and 
<math>\tau_2(e_{r,(i,j,k)}) \approx \left\{\begin{array}{cc} l_\theta(i+\frac{1}{2}) \cdot l_\phi(i+\frac{1}{2},j) & (0 < j < N_\theta) \\ \pi(\frac{1}{2}l_\theta(i+\frac{1}{2}))^2 & (j = 0 \;\text{or}\; j = N_\theta)
 \\ \end{array}\right.</math> respectively.

For each <math>(i,j,k) \in I_\theta</math> the length and area of edge <math>e_{\theta,(i,j,k)}</math> is <math>\tau_1(e_{\theta,(i,j,k)}) \approx l_\theta(i)</math> and <math>\tau_2(e_{\theta,(i,j,k)}) \approx \Delta r \cdot l_\phi(i,j+\frac{1}{2})</math> respectively.

For each <math>(i,j,k) \in I_\phi</math> the length and area of edge <math>e_{\phi,(i,j,k)}</math> is <math>\tau_1(e_{\phi,(i,j,k)}) \approx l_\phi(i,j)</math> and <math>\tau_2(e_{\phi,(i,j,k)}) = \Delta r \cdot l_\theta(i)</math> respectively.

Given a scalar field <math>f(r,\theta,\phi)</math> over spherical coordinates, <math>f</math> will be approximated by the node based function <math>f_\text{approx} : \mathcal{N} \to \R</math>. Function <math>f_\text{approx}</math> is defined by <math>f_\text{approx}(n_{(i,j,k)}) = f(i\Delta r, j\Delta\theta, k\Delta\phi)</math> for each <math>(i,j,k) \in I_n</math>.

Given a vector field <math>\mathbf{F}(r,\theta,\phi) = F_r(r,\theta,\phi)\hat{\mathbf{r}} + F_\theta(r,\theta,\phi)\hat{\mathbf{\theta}} + F_\phi(r,\theta,\phi)\hat{\mathbf{\phi}}</math> over spherical coordinates, <math>\mathbf{F}</math> will be approximated by the edge based function <math>F_\text{approx} : \mathcal{E} \to \R</math>. Function <math>F_\text{approx}</math> is defined by:

<math>\forall (i,j,k) \in I_r : F_\text{approx}(e_{r,(i,j,k)}) = \left\{\begin{array}{cc} F_r(i\Delta r, j\Delta\theta, k\Delta\phi) & (i > 0) \\ \lim_{r \to 0^+} F_r(r, j\Delta\theta, k\Delta\phi) & (i = 0) \\ \end{array}\right.</math>

<math>\forall (i,j,k) \in I_\theta : F_\text{approx}(e_{\theta,(i,j,k)}) = \left\{\begin{array}{cc} F_\theta(i\Delta r, j\Delta\theta, k\Delta\phi) & (j > 0) \\ \lim_{\theta \to 0^+} F_\theta(i\Delta r, \theta, k\Delta\phi) & (j = 0) \\ \end{array}\right.</math>

<math>\forall (i,j,k) \in I_\phi : F_\text{approx}(e_{\phi,(i,j,k)}) = F_\phi(i\Delta r, j\Delta\theta, k\Delta\phi)</math>

=== Model of Curvilinear coordinates ===

A curvilinear coordinate system is a coordinate system such as cylindrical or spherical coordinates where the position may not follow a straight line as one of the coordinates is changed. To model an arbitrary curvilinear coordinate system in 3 dimensional space, let the coordinates be denoted by <math>\xi_1, \xi_2, \xi_3</math>.

To avoid unnecessary repetition, the following notation will be used:

* Given any number <math>k \in \R</math>, and subset <math>C \subseteq \{1,2,3\}</math>, then <math>[k]_C</math> denotes the triple <math>[k]_C = \left(\left\{\begin{array}{cc} k & (1 \in C) \\ 0 & (1 \notin C) \end{array}\right.,\left\{\begin{array}{cc} k & (2 \in C) \\ 0 & (2 \notin C) \end{array}\right.,\left\{\begin{array}{cc} k & (3 \in C) \\ 0 & (3 \notin C) \end{array}\right.\right)</math>
* Given any two triples of numbers <math>(x_1,x_2,x_3), (y_1,y_2,y_3) \in \R^3</math>, then <math>(x_1,x_2,x_3) + (y_1,y_2,y_3) = (x_1+y_1,x_2+y_2,x_3+y_3)</math>
* Given any two triples of numbers <math>(x_1,x_2,x_3), (y_1,y_2,y_3) \in \R^3</math>, then <math>(x_1,x_2,x_3)(y_1,y_2,y_3) = (x_1y_1,x_2y_2,x_3y_3)</math>
*<math>\Xi = (\xi_1,\xi_2,\xi_3)</math> denotes an arbitrary curvilinear coordinate <math>(\xi_1,\xi_2,\xi_3)</math>.
*<math>\Delta\Xi = (\Delta\xi_1,\Delta\xi_2,\Delta\xi_3)</math> denotes an arbitrary triplet of differences in the coordinates <math>\xi_1,\xi_2,\xi_3</math>.

Let <math>\omega_0(\Xi)</math> denote the position given by coordinate <math>\Xi = (\xi_1,\xi_2,\xi_3)</math>. 
For each <math>c = 1,2,3</math>, the rate of change in the position <math>\omega_0(\Xi)</math> at coordinate <math>\Xi</math> with respect to <math>\xi_c</math> is given by: <math>\mathbf{u}_{\xi_c}(\Xi) = \frac{\partial\omega_0}{\partial\xi_c}\bigg|_{\Xi}</math>. 
It will be assumed that <math>\mathbf{u}_{\xi_1}</math>, <math>\mathbf{u}_{\xi_2}</math>, and <math>\mathbf{u}_{\xi_3}</math> are all mutually perpendicular, which is the case for Cartesian, cylindrical, and spherical coordinates. 
In addition, for each <math>c = 1,2,3</math>, vector <math>\mathbf{v}_{\xi_c}(\Xi) = \frac{\mathbf{u}_{\xi_c}(\Xi)}{|\mathbf{u}_{\xi_c}(\Xi)|}</math> will also denote a unit length normalization of <math>\mathbf{u}_{\xi_c}(\Xi)</math>.

For the purposes of simplicity, it will be assumed that the coordinates <math>\Xi = (\xi_1, \xi_2, \xi_3)</math> can be any triple of real numbers. Singularities such as the z-axis in cylindrical and spherical coordinates will not be considered by this model. Since singularities are not considered, this model will be considerably simpler than the models given for cylindrical and spherical coordinates.

To model the curvilinear coordinate system, a lattice of nodes <math>\mathcal{N}</math> connected by directed edges <math>\mathcal{E}</math> will be established:

Start by choosing "resolutions" for each coordinate <math>\xi_1, \xi_2, \xi_3</math> denoted by <math>\Delta\xi_1, \Delta\xi_2, \Delta\xi_3</math>. The quantities <math>\Delta\xi_1, \Delta\xi_2, \Delta\xi_3</math> should all be strictly positive, and ideally should be as small as possible.

*For each triplet of integer indices <math>I \in \Z^3</math>, a node <math>n_I \in \mathcal{N}</math> will be created at position <math>\omega_0(n_I) =  \omega_0(I\Delta\Xi)</math>.
*For each <math>c = 1,2,3</math> and for each triplet of integer indices <math>I \in \Z^3</math>, there will exist edge <math>e_{\xi_c,I} \in \mathcal{E}</math> such that <math>t_\bullet(e_{\xi_c,I}) = n_I</math> and <math>t^\bullet(e_{\xi_c,I}) = n_{I+[1]_c}</math>.
*For each <math>I \in \Z^3</math>, the volume of node <math>n_I</math> is <math>\tau_3(n_I) = \prod_{c=1}^3 (\Delta\xi_c|\mathbf{u}_{\xi_c}(I\Delta\Xi)|)</math> <math> = \left(\prod_{c=1}^3|\mathbf{u}_{\xi_c}(I\Delta\Xi)|\right)\left(\prod_{c=1}^3 \Delta\xi_c\right)</math>
*For each <math>c = 1,2,3</math> and <math>I \in \Z^3</math>, the length and area of edge <math>e_{\xi_c,I}</math> are <math>\tau_1(e_{\xi_c,I}) = \Delta\xi_c|\mathbf{u}_{\xi_c}((I+[0.5]_c)\Delta\Xi)|</math> and <math>\tau_2(e_{\xi_c,I}) = \prod_{c'\neq c}\Delta\xi_{c'}|\mathbf{u}_{\xi_{c'}}((I+[0.5]_c)\Delta\Xi)|</math> respectively.

Any node <math>n \in \mathcal{N}</math> for which <math>\tau_3(n) = 0</math> should be deleted along with all connecting edges.

For any edge <math>e \in \mathcal{E}</math> for which <math>\tau_1(e) = 0</math> and <math>\tau_2(e) > 0</math>, then nodes <math>t_\bullet(e)</math> and <math>t^\bullet(e)</math> should be fused into a single node, and edge <math>e</math> should be discarded. For any edge <math>e \in \mathcal{E}</math> for which <math>\tau_2(e) = 0</math> and <math>\tau_1(e) > 0</math>, then edge <math>e</math> should be discarded. Either action may be pursued when <math>\tau_1(e)=\tau_2(e)=0</math>.

Given a scalar field <math>f(\Xi)</math> over the curvilinear coordinate system, <math>f</math> will be approximated by the node based function <math>f_\text{approx} : \mathcal{N} \to \R</math>. Function <math>f_\text{approx}</math> is defined by <math>f_\text{approx}(n_I) = f(I\Delta\Xi)</math> for each <math>I \in \Z^3</math>.

Given a vector field <math>\mathbf{F}(\Xi) = \sum_{c=1}^3 F_{\xi_c}(\Xi)\mathbf{v}_{\xi_c}(\Xi)</math> over the curvilinear coordinates, <math>\mathbf{F}</math> will be approximated by the edge based function <math>F_\text{approx} : \mathcal{E} \to \R</math>. Function <math>F_\text{approx}</math> is defined by <math>F_\text{approx}(e_{\xi_c,I}) = F_{\xi_c}((I+[0.5]_c)\Delta\Xi)</math> for each <math>c = 1,2,3</math> and for each <math>I \in \Z^3</math>.

== Discrete analog to Integrals ==

=== Volume Integrals ===

Let <math>U \subseteq \mathcal{N}</math> denote an arbitrary volume. Given a node based function <math>f</math>, the "volume integral" of <math>f</math> over <math>U</math> is <math>\int_{n \in U} f(n) = \sum_{n \in U} f(n)\tau_3(n)</math> (recall that <math>\tau_3(n)</math> is the volume of node <math>n</math>, and <math>\omega_3(n)</math> is the volume of space represented by <math>n</math>).

If node based function <math>f : \mathcal{N} \to \R</math> is an approximation of a scalar field <math>f_\text{continuous} : \bigg(\bigcup_{n \in \mathcal{N}} \omega_3(n)\bigg) \to \R</math>, where <math>\forall n \in \mathcal{N} : f(n) = f_\text{continuous}(\omega_0(n))</math>, then <math>\int_{n \in U} f(n)</math> is an approximation of the volume integral <math>\iiint_{\mathbf{q} \in \bigcup_{n \in U} \omega_3(n)} f(\mathbf{q})dV</math>.

Recall that <math>U</math> can be described via the node based function <math>\delta_3(n;U) = \left\{\begin{array}{cc} 1 & (n \in U) \\ 0 & (n \notin U) \\ \end{array}\right.</math>. This allows the expression <math>\int_{n \in U} f(n) = \sum_{n \in \mathcal{N}} f(n)\delta_3(n;U)\tau_3(n)</math>. If node based function <math>g</math> denotes a multi-volume, then <math>\int_{n \in g} f(n) = \sum_{n \in \mathcal{N}} f(n)g(n)\tau_3(n) = \Sigma_\mathcal{N}(f,g)</math>.

<math>f</math> most often denotes a density function, and <math>\int_{n \in g} f(n)</math> denotes the total "mass" contained by multi-volume <math>g</math>. The equality <math>\int_{n \in g} f(n) = \Sigma_\mathcal{N}(f,g)</math> means that if <math>f</math> is interpreted as a multi-point, then the total "mass" contained by <math>g</math> is the total intersection weight of <math>f</math> with <math>g</math>.

=== Path Integrals ===

Let <math>P = \langle (e_1,s_1), (e_2,s_2), ..., (e_k,s_k) \rangle</math> denote an arbitrary path. Given an edge based function <math>f</math>, the "path integral" of <math>f</math> along <math>P</math> is <math>\int_{e \in P} f(e) = \sum_{i=1}^k s_if(e_i)\tau_1(e_i)</math> (recall that <math>\tau_1(e)</math> is the length of edge <math>e</math>, and <math>\omega_1(e)</math> is the curve represented by <math>e</math>).

If edge based function <math>f : \mathcal{E} \to \R</math> is an approximation of a vector field <math>\mathbf{F}_\text{continuous} : \bigg(\bigcup_{n \in \mathcal{N}} \omega_3(n)\bigg) \to \R^3</math>, where <math>\forall e \in \mathcal{E} : f(e) = \mathbf{F}_\text{continuous}(\omega_0(t_\bullet(e))) \cdot \frac{\omega_0(t^\bullet(e))-\omega_0(t_\bullet(e))}{|\omega_0(t^\bullet(e))-\omega_0(t_\bullet(e))|}</math>, then <math>\int_{e \in P} f(e)</math> is an approximation of the path integral <math>\int_{\mathbf{q} \in \bigcup_{i=1}^k s_i\omega_1(e_i)} \mathbf{F}(\mathbf{q}) \cdot d\mathbf{q}</math>.

Recall that <math>P</math> can be described via the edge based function <math>\delta_1(e;P) = \sum_{i=1}^k s_i\left\{\begin{array}{cc} 1/\tau_2(e) & (e = e_i) \\ 0 & (e \neq e_i) \\ \end{array}\right.</math>. This allows the expression <math>\int_{e \in P} f(e) = \sum_{e \in \mathcal{E}} f(e)\delta_1(e;P)\tau_2(e)\tau_1(e)</math>. If edge based function <math>g</math> denotes a multi-path, then <math>\int_{e \in g} f(e) = \sum_{e \in \mathcal{E}} f(e)g(e)\tau_2(e)\tau_1(e) = \Sigma_\mathcal{E}(g,f)</math>.

<math>f</math> most often denotes a rate of gain function, and <math>\int_{e \in g} f(e)</math> denotes the total "gain" caused by <math>f</math> along multi-path <math>g</math>. The equality <math>\int_{e \in g} f(e) = \Sigma_\mathcal{E}(g,f)</math> means that if <math>f</math> is interpreted as a multi-surface, then the total "gain" generated along <math>g</math> is the total intersection weight of <math>g</math> with <math>f</math>.

=== Surface Integrals ===

Let <math>S = \{(e_1,s_1), (e_2,s_2), ..., (e_k,s_k)\}</math> denote an arbitrary surface. Given an edge based function <math>f</math>, the "surface integral" of <math>f</math> over <math>S</math> is <math>\int_{e \in S} f(e) = \sum_{(e,s) \in S} sf(e)\tau_2(e)</math> (recall that <math>\tau_2(e)</math> is the cross-sectional area of edge <math>e</math>, and <math>\omega_2(e)</math> is the surface which separates <math>\omega_3(t_\bullet(e))</math> and <math>\omega_3(t^\bullet(e))</math>).

If edge based function <math>f : \mathcal{E} \to \R</math> is an approximation of a vector field <math>\mathbf{F}_\text{continuous} : \bigg(\bigcup_{n \in \mathcal{N}} \omega_3(n)\bigg) \to \R^3</math>, where <math>\forall e \in \mathcal{E} : f(e) = \mathbf{F}_\text{continuous}(\omega_0(t_\bullet(e))) \cdot \frac{\omega_0(t^\bullet(e))-\omega_0(t_\bullet(e))}{|\omega_0(t^\bullet(e))-\omega_0(t_\bullet(e))|}</math>, then <math>\int_{e \in S} f(e)</math> is an approximation of the surface integral <math>\iint_{\mathbf{q} \in \bigcup_{(e,s) \in S} s\omega_2(e)} \mathbf{F}(\mathbf{q}) \cdot \mathbf{dS}</math>.

Recall that <math>S</math> can be described via the edge based function <math>\delta_2(e;S) = \left\{\begin{array}{cc} +1/\tau_1(e) & ((e,+1) \in S) \\ -1/\tau_1(e) & ((e,-1) \in S) \\ 0 & (\text{else}) \\ \end{array}\right.</math>. This allows the expression <math>\int_{e \in S} f(e) = \sum_{e \in \mathcal{E}} f(e)\delta_2(e;S)\tau_2(e)\tau_1(e)</math>. If edge based function <math>g</math> denotes a multi-surface, then <math>\int_{e \in g} f(e) = \sum_{e \in \mathcal{E}} f(e)g(e)\tau_2(e)\tau_1(e) = \Sigma_\mathcal{E}(f,g)</math>.

<math>f</math> most often denotes a flow density function, and <math>\int_{e \in g} f(e)</math> denotes the total "flow" through multi-surface <math>g</math>. The equality <math>\int_{e \in g} f(e) = \Sigma_\mathcal{E}(f,g)</math> means that if <math>f</math> is interpreted as a multi-path, then the total "flow" through <math>g</math> is the total intersection weight of <math>f</math> with <math>g</math>.

== Discrete analog to the Gradient ==

=== Computing the surfaces of volumes ===

[[File:Discrete gradient multi surface inwards orientation.svg|thumb|Given an arbitrary multi-volume denoted by node based function <math>f</math>, the multi-surface which is the inwards oriented surface of the multi-volume is denoted by the edge based function <math>g(e) = \nabla f|_e</math>. When traversing an arbitrary edge, the net number (positive - negative) of transitions through surfaces is the change in <math>f</math>.]]

Consider an arbitrary volume <math>U \subseteq \mathcal{N}</math>. The inwards oriented surface of <math>U</math> is a surface <math>S</math> defined by <math>(e,+1) \in S \iff (t_\bullet(e) \notin U \;\text{and}\; t^\bullet(e) \in U)</math> and <math>(e,-1) \in S \iff (t_\bullet(e) \in U \;\text{and}\; t^\bullet(e) \notin U)</math>. Volume <math>U</math> is denoted by the node based function <math>\delta_3(n;U) = \left\{\begin{array}{cc} 1 & (n \in U) \\ 0 & (n \notin U) \end{array}\right.</math>. The inwards oriented surface <math>S</math> is denoted by the edge based function: 
<math>\delta_2(e;S) = \left\{\begin{array}{cc} +1/\tau_1(e) & ((e,+1) \in S) \\ -1/\tau_1(e) & ((e,-1) \in S) \\ 0 & (\text{else}) \end{array}\right. </math> 
<math> = \left\{\begin{array}{cc} +1/\tau_1(e) & (t_\bullet(e) \notin U \;\text{and}\; t^\bullet(e) \in U) \\ -1/\tau_1(e) & (t_\bullet(e) \in U \;\text{and}\; t^\bullet(e) \notin U) \\ 0 & (\text{else}) \end{array}\right. </math>
<math> = \frac{\delta_3(t^\bullet(e);U) - \delta_3(t_\bullet(e);U)}{\tau_1(e)} </math>

More generally, if node based function <math>f</math> denotes a  multi-volume, then the edge based function <math>g(e) = \frac{f(t^\bullet(e))-f(t_\bullet(e))}{\tau_1(e)}</math> is the inwards oriented multi-surface boundary of <math>f</math>. One required property of the inwards oriented multi-surface boundary is that it must obey linearity: meaning that if multi-volumes <math>f_1</math> and <math>f_2</math> have inwards oriented multi-surface boundaries <math>g_1</math> and <math>g_2</math> respectively, then the inwards oriented multi-surface boundary of <math>f_1 + f_2</math> is <math>g_1 + g_2</math>. In addition, given an arbitrary coefficients <math>c_1,c_2</math>, the inwards oriented multi-surface boundary of <math>c_1f_1+c_2f_2</math> is <math>c_1g_1+c_2g_2</math>. The definition <math>g(e) = \frac{f(t^\bullet(e))-f(t_\bullet(e))}{\tau_1(e)}</math> satisfies this linearity condition.

Given a node based function <math>f</math>, the "gradient" of <math>f</math> is the edge based function <math>\nabla f</math>. For all <math>e \in \mathcal{E}</math>, <math>\nabla f|_e = \frac{f(t^\bullet(e))-f(t_\bullet(e))}{\tau_1(e)}</math>. The gradient denotes both the inwards oriented surface when <math>f</math> is treated as a multi-volume, and the rate of change in <math>f</math> along each edge <math>e</math> in the preferred direction when <math>f</math> is treated as a potential. Note that the gradient satisfies linearity as described above.

=== Computing potential differences from the gradient ===

[[File:Discrete gradient path integral.svg|thumb|The path integral along path <math>P</math> of the gradient of <math>f</math> is the difference in <math>f</math> between the endpoints. The value of <math>f</math> at each node is shown, and along the path the value of <math>s_i\tau_1(e_i)\nabla f|_{e_i}</math> is shown at each edge. (the edge orientations are not shown for simplicity)]]

Let <math>n_0,n_1 \in \mathcal{N}</math> be two nodes that are linked by edge <math>e</math>. Let <math>n_1</math> be reached from <math>n_0</math> by traversing edge <math>e</math> in the direction specified by orientation <math>s \in \{+1,-1\}</math>.

The difference in <math>f</math> between <math>n_0</math> and <math>n_1</math> is:

<math>f(n_1) - f(n_0) = \left\{\begin{array}{cc} f(t^\bullet(e)) - f(t_\bullet(e)) & (s = +1) \\ f(t_\bullet(e)) - f(t^\bullet(e)) & (s = -1) \end{array}\right. = s(f(t^\bullet(e)) - f(t_\bullet(e))) = s\tau_1(e)\frac{f(t^\bullet(e)) - f(t_\bullet(e))}{\tau_1(e)} = s\tau_1(e)\nabla f|_e</math>

When given an arbitrary path <math>P = \langle (e_1,s_1), (e_2,s_2), ..., (e_k,s_k) \rangle</math> that begins at node <math>n_0</math> and ends at node <math>n_1</math>, the difference in <math>f</math> between <math>n_0</math> and <math>n_1</math> can be computed by adding the differences along each edge:

<math>f(n_1) - f(n_0) = \sum_{i=1}^k s_i\tau_1(e_i)\nabla f|_{e_i} = \int_{e \in P} \nabla f|_e </math>

The equation <math> \int_{e \in P} \nabla f|_e = f(n_1) - f(n_0) </math> is referred to as the "gradient theorem".

The following sections will demonstrate that the gradient defined for node based functions is consistent with the gradient for scalar fields in <math>\R^3</math>. This consistency will be shown for Cartesian, cylindrical, and spherical coordinates.

=== The Gradient in Various Coordinate Systems ===

All of Cartesian, cylindrical, and spherical coordinate systems will be modeled using the [[Calculus/Discrete vector calculus#Model of Curvilinear coordinates|discrete model of curvilinear coordinates]]. The same shorthand notation will be used.

Starting with an arbitrary scalar field <math>f(\Xi)</math>, the node based function that approximates <math>f</math> is <math>f_\text{approx}(n_I) = f(I\Delta\Xi)</math> for all <math>I \in \Z^3</math>.

For each <math>c = 1,2,3</math> and <math>I \in \Z^3</math>, the gradient along edge <math>e_{\xi_c,I}</math> is:

<math>\nabla f_\text{approx}|_{e_{\xi_c,I}} = \frac{f_\text{approx}(t^\bullet(e_{\xi_c,I})) - f_\text{approx}(t_\bullet(e_{\xi_c,I}))}{\tau_1(e_{\xi_c,I})}</math>
<math> = \frac{f_\text{approx}(n_{I+[1]_c}) - f_\text{approx}(n_I)}{\Delta\xi_c|\mathbf{u}_{\xi_c}(I\Delta\Xi)|}</math>
<math> = \frac{f(I\Delta\Xi + [\Delta\xi_c]_c) - f(I\Delta\Xi)}{\Delta\xi_c|\mathbf{u}_{\xi_c}((I+[0.5]_c)\Delta\Xi)|}</math>

Given an arbitrary coordinate <math>\Xi</math>, let <math>I(\Xi) \in \Z^3</math> index the node <math>n_{I(\Xi)}</math> whose position <math>I(\Xi)\Delta\Xi</math> is "closest" to <math>\Xi</math>.

The gradient of scalar field <math>f</math> is a vector quantity with components that are multiples of the basis vectors <math>\mathbf{v}_{\xi_1}, \mathbf{v}_{\xi_2}, \mathbf{v}_{\xi_3}</math>. The coefficient of <math>\mathbf{v}_{\xi_1}</math> is the gradient of <math>f_\text{approx}</math> along edges that are parallel to <math>\mathbf{v}_{\xi_1}</math>, and similarly for <math>\mathbf{v}_{\xi_2}</math> and <math>\mathbf{v}_{\xi_3}</math>. The fact that <math>\mathbf{v}_{\xi_1}</math>, <math>\mathbf{v}_{\xi_2}</math>, and <math>\mathbf{v}_{\xi_3}</math> are of unit length and mutually orthogonal (perpendicular) is important.

The gradient of scalar field <math>f</math> at <math>\Xi</math> can now be approximated by the following:

<math>\nabla f|_{\Xi} \approx \sum_{c=1}^3 \nabla f_\text{approx}|_{e_{\xi_c,I(\Xi)}}\mathbf{v}_{\xi_c}(I(\Xi)\Delta\Xi)</math>
<math> = \sum_{c=1}^3 \frac{f(I(\Xi)\Delta\Xi + [\Delta\xi_c]_c) - f(I(\Xi)\Delta\Xi)}{\Delta\xi_c|\mathbf{u}_{\xi_c}((I(\Xi)+[0.5]_c)\Delta\Xi)|}\mathbf{v}_{\xi_c}(I(\Xi)\Delta\Xi)</math>
<math> \approx \sum_{c=1}^3 \frac{f(\Xi + [\Delta\xi_c]_c) - f(\Xi)}{\Delta\xi_c|\mathbf{u}_{\xi_c}(\Xi)|}\mathbf{v}_{\xi_c}(\Xi)</math>

As <math>\Delta\xi_1, \Delta\xi_2, \Delta\xi_3 \to 0^+</math>, it becomes apparent that <math>\nabla f|_{\Xi} = \sum_{c=1}^3 \frac{1}{|\mathbf{u}_{\xi_c}(\Xi)|}\frac{\partial f}{\partial \xi_c}\bigg|_{\Xi}\mathbf{v}_{\xi_c}(\Xi)</math>.

==== The Gradient in Cartesian Coordinates ====

Applying the curvilinear coordinate model to Cartesian coordinates means that <math>\xi_1 = x</math>, <math>\xi_2 = y</math>, and <math>\xi_3 = z</math>. It is also the case that <math>\mathbf{u}_x(x,y,z) = \mathbf{i}</math>, <math>\mathbf{u}_y(x,y,z) = \mathbf{j}</math>, and <math>\mathbf{u}_z(x,y,z) = \mathbf{k}</math>. Applying the formula for the gradient that was derived from the discrete model of curvilinear coordinates gives:

<math>\nabla f|_{(x,y,z)} = \frac{1}{1}\frac{\partial f}{\partial x}\bigg|_{(x,y,z)}\mathbf{i} + \frac{1}{1}\frac{\partial f}{\partial y}\bigg|_{(x,y,z)}\mathbf{j} + \frac{1}{1}\frac{\partial f}{\partial z}\bigg|_{(x,y,z)}\mathbf{k}</math> 
<math> = \frac{\partial f}{\partial x}\mathbf{i} + \frac{\partial f}{\partial y}\mathbf{j} + \frac{\partial f}{\partial z}\mathbf{k}</math>

which is consistent with the formula for the gradient derived for continuous functions over Cartesian coordinates.

==== The Gradient in Cylindrical Coordinates ====

Applying the curvilinear coordinate model to cylindrical coordinates means that <math>\xi_1 = \rho</math>, <math>\xi_2 = \phi</math>, and <math>\xi_3 = z</math>. It is also the case that <math>\mathbf{u}_\rho(\rho,\phi,z) = \hat{\mathbf{\rho}}</math>, <math>\mathbf{u}_\phi(\rho,\phi,z) = \rho\hat{\mathbf{\phi}}</math>, and <math>\mathbf{u}_z(\rho,\phi,z) = \hat{\mathbf{z}}</math>. Applying the formula for the gradient that was derived from the discrete model of curvilinear coordinates gives:

<math>\nabla f|_{(\rho,\phi,z)} = \frac{1}{1}\frac{\partial f}{\partial\rho}\bigg|_{(\rho,\phi,z)}\hat{\mathbf{\rho}} + \frac{1}{\rho}\frac{\partial f}{\partial\phi}\bigg|_{(\rho,\phi,z)}\hat{\mathbf{\phi}} + \frac{1}{1}\frac{\partial f}{\partial z}\bigg|_{(\rho,\phi,z)}\hat{\mathbf{z}}</math> 
<math> = \frac{\partial f}{\partial\rho}\hat{\mathbf{\rho}} + \frac{1}{\rho}\frac{\partial f}{\partial\phi}\hat{\mathbf{\phi}} + \frac{\partial f}{\partial z}\hat{\mathbf{z}}</math>

which is consistent with the formula for the gradient derived for continuous functions over cylindrical coordinates.

==== The Gradient in Spherical Coordinates ====

Applying the curvilinear coordinate model to spherical coordinates means that <math>\xi_1 = r</math>, <math>\xi_2 = \theta</math>, and <math>\xi_3 = \phi</math>. It is also the case that <math>\mathbf{u}_r(r,\theta,\phi) = \hat{\mathbf{r}}</math>, <math>\mathbf{u}_\theta(r,\theta,\phi) = r\hat{\mathbf{\theta}}</math>, and <math>\mathbf{u}_\phi(r,\theta,\phi) = r\sin\theta\hat{\mathbf{\phi}}</math>. Applying the formula for the gradient that was derived from the discrete model of curvilinear coordinates gives:

<math>\nabla f|_{(r,\theta,\phi)} = \frac{1}{1}\frac{\partial f}{\partial r}\bigg|_{(r,\theta,\phi)}\hat{\mathbf{r}} + \frac{1}{r}\frac{\partial f}{\partial\theta}\bigg|_{(r,\theta,\phi)}\hat{\mathbf{\theta}} + \frac{1}{r\sin\theta}\frac{\partial f}{\partial\phi}\bigg|_{(r,\theta,\phi)}\hat{\mathbf{\phi}}</math> 
<math> = \frac{\partial f}{\partial r}\hat{\mathbf{r}} + \frac{1}{r}\frac{\partial f}{\partial\theta}\hat{\mathbf{\theta}} + \frac{1}{r\sin\theta}\frac{\partial f}{\partial\phi}\hat{\mathbf{\phi}}</math>

which is consistent with the formula for the gradient derived for continuous functions over spherical coordinates.

== Discrete analog to the Divergence ==

=== Flow generation density ===

[[File:Discrete divergence multi path.svg|thumb|When an edge based function is depicted as a superposition of paths, the divergence denotes the starting point density minus the end point density. In this image, the flow is depicted using the brown arrows along each edge. More arrows means more flow. The numbers indicate the total flow generation at each node, which is also the net starting point weight.]]

If the edge based function <math>f: \mathcal{E} \to \R</math> denotes the "flow density" across each edge, then for an arbitrary edge <math>e \in \mathcal{E}</math>, the total flow across <math>e</math> in the preferred direction is <math>\tau_2(e)f(e)</math>. For an arbitrary node <math>n \in \mathcal{N}</math>, the total flow being generated at node <math>n</math> is <math>\sum_{e \in \mathcal{E}} \left\{\begin{array}{cc} +\tau_2(e)f(e) & (t_\bullet(e) = n \;\text{and}\; t^\bullet(e) \neq n) \\ -\tau_2(e)f(e) & (t^\bullet(e) = n \;\text{and}\; t_\bullet(e) \neq n) \\ 0 & (\text{else}) \end{array}\right.</math>. The "flow generation density" or "divergence" at node <math>n</math> is:

<math>\nabla \cdot f|_n = \frac{1}{\tau_3(n)}\sum_{e \in \mathcal{E}} \left\{\begin{array}{cc} +\tau_2(e)f(e) & (t_\bullet(e) = n \;\text{and}\; t^\bullet(e) \neq n) \\ -\tau_2(e)f(e) & (t^\bullet(e) = n \;\text{and}\; t_\bullet(e) \neq n) \\ 0 & (\text{else}) \end{array}\right.</math>

Using a similar argument, given a multi-path denoted by edge based function <math>f: \mathcal{E} \to \R</math>, then for an arbitrary edge <math>e \in \mathcal{E}</math>, the total "path weight" that traverses <math>e</math> in the preferred direction is <math>\tau_2(e)f(e)</math> (recall that <math>f(e)</math> is the path weight per unit thickness of <math>e</math>). The total "path weight" that is "left behind" at an arbitrary node <math>n \in \mathcal{N}</math> is <math>\sum_{e \in \mathcal{E}} \left\{\begin{array}{cc} +\tau_2(e)f(e) & (t_\bullet(e) = n \;\text{and}\; t^\bullet(e) \neq n) \\ -\tau_2(e)f(e) & (t^\bullet(e) = n \;\text{and}\; t_\bullet(e) \neq n) \\ 0 & (\text{else}) \end{array}\right.</math>. The path weight that is left behind is the starting point weight minus the ending point weight, no matter how multi-path <math>f</math> is decomposed into a superposition of paths. The starting point weight minus the ending point weight is the "'''net starting point weight'''". By diffusing the net starting point weight over the volume of node <math>n</math>, the node based function that denotes the "'''net starting point'''" is derived. This is the divergence <math>\nabla \cdot f</math>.

It is a relatively simple matter to show that for an a simple path <math>P</math> that starts at node <math>n_0</math> and ends at node <math>n_1</math>, that <math>\nabla \cdot \delta_1(P) = \delta_0(n_0) - \delta_0(n_1)</math> due to the fact that only the end points of <math>P</math> may have a non-zero starting point weight. <math>\delta_0(n_0) - \delta_0(n_1)</math> is a node based function that denotes the superposition of the starting node of <math>P</math> with the final node of <math>P</math> with respective weights of <math>+1</math> and <math>-1</math>.

=== The divergence theorem ===

[[File:Discrete divergence volume integral.svg|thumb|A depiction of the discrete analog of Gauss's divergence theorem. The rate (<math>\tau_2(e)f(e)</math>) and direction of the flow along each edge is shown by the number and direction of the arrows, and each node is labeled with the rate of flow generation (<math>\tau_3(n)(\nabla \cdot f|_n)</math>). A net flow of -4 is flowing out from the highlighted boundary, and the total flow generation contributed by the nodes inside the boundary is also -4.]]

Let <math>U \subseteq \mathcal{N}</math> denote an arbitrary volume, and let <math>S = \{(e,+1) | e \in \mathcal{E} \;\text{and}\; t_\bullet(e) \in U \;\text{and}\; t^\bullet(e) \notin U\} \cup \{(e,-1) | e \in \mathcal{E} \;\text{and}\; t_\bullet(e) \notin U \;\text{and}\; t^\bullet(e) \in U\}</math> denote the outwards oriented surface of <math>U</math>.

The total flow generation inside <math>U</math> is:

<math>\int_{n \in U} \nabla \cdot f|_n = \sum_{n \in U} (\nabla \cdot f|_n)\tau_3(n)</math>
<math> = \sum_{n \in U} \sum_{e \in \mathcal{E}} \left\{\begin{array}{cc} +\tau_2(e)f(e) & (t_\bullet(e) = n \;\text{and}\; t^\bullet(e) \neq n) \\ -\tau_2(e)f(e) & (t^\bullet(e) = n \;\text{and}\; t_\bullet(e) \neq n) \\ 0 & (\text{else}) \end{array}\right.</math>
<math> = \sum_{e \in \mathcal{E}} \sum_{n \in U} \left\{\begin{array}{cc} +\tau_2(e)f(e) & (t_\bullet(e) = n \;\text{and}\; t^\bullet(e) \neq n) \\ -\tau_2(e)f(e) & (t^\bullet(e) = n \;\text{and}\; t_\bullet(e) \neq n) \\ 0 & (\text{else}) \end{array}\right.</math>
<math> = \sum_{e \in \mathcal{E}} \left\{\begin{array}{cc} +\tau_2(e)f(e) & (t_\bullet(e) \in U \;\text{and}\; t^\bullet(e) \notin U) \\ -\tau_2(e)f(e) & (t^\bullet(e) \in U \;\text{and}\; t_\bullet(e) \notin U) \\ 0 & (\text{else}) \end{array}\right.</math>
<math> = \sum_{e \in \mathcal{E}} \tau_2(e)f(e)\left\{\begin{array}{cc} +1 & ((e,+1) \in S) \\ -1 & ((e,-1) \in S) \\ 0 & (\text{else}) \end{array}\right.</math>
<math> = \sum_{(e,s) \in S} sf(e)\tau_2(e)</math>
<math> = \int_{e \in S} f(e)</math>

The total flow generated inside <math>U</math> is the total outwards flow through <math>S</math>: <math>\int_{n \in U} \nabla \cdot f|_n = \int_{e \in S} f(e)</math>. This is the discrete analog to Gauss's Divergence Theorem.

[[File:Point volume intersection vs path surface intersection.svg|thumb|Given an arbitrary multi-volume and a multi-path, the net intersection weight of the path's net starting point with the volume is equal to the net intersection weight of the path with the volume's outwards oriented surface. In this case the net intersection is 3.]]

Interpreting <math>f</math> as a multi-path, and <math>\nabla \cdot f</math> as the net starting point of <math>f</math>, then Gauss's Divergence Theorem can be expressed as: <math>\Sigma_\mathcal{N}(\nabla \cdot f, \delta_3(U)) = \Sigma_\mathcal{E}(f, \delta_2(S))</math>. Due to the fact that <math>S</math> is the outwards oriented surface of <math>U</math>, it must be the case that <math>\delta_2(S) = -\nabla \delta_3(U)</math>, since <math>\nabla \delta_3(U)</math> is the inwards oriented surface of <math>U</math>. This gives: 
 
<math>\Sigma_\mathcal{N}(\nabla \cdot f, \delta_3(U)) = \Sigma_\mathcal{E}(f, -\nabla\delta_3(U))</math>

Replacing <math>U</math> with a multi-volume <math>g</math> gives the generalization:

<math>\Sigma_\mathcal{N}(\nabla \cdot f, g) = \Sigma_\mathcal{E}(f, -\nabla g)</math>

This equation can be interpreted as follows: '''The total intersection of the net starting point of multi-path <math>f</math> with multi-volume <math>g</math> is the net intersection of <math>f</math> with the outwards oriented surface of <math>g</math>. '''

It should also be noted that the equation can be reformatted to read as:

<math>\int_{e \in f} \nabla g|_e = \Sigma_\mathcal{N}(-\nabla \cdot f, g)</math>

When <math>f</math> is a simple path <math>P</math> that starts at node <math>n_0</math> and ends at node <math>n_1</math>, the above equation reduces to the gradient theorem:

<math>\int_{e \in P} \nabla g|_e = g(n_1) - g(n_0)</math>

The gradient theorem and Gauss's Divergence Theorem are therefore equivalent with the general form: <math>\Sigma_\mathcal{N}(\nabla \cdot f, g) = \Sigma_\mathcal{E}(f, -\nabla g)</math>.

The following sections will demonstrate that the divergence defined for edge based functions is consistent with the divergence for vector fields in <math>\mathbb{R}^3</math>. This consistency will be shown for Cartesian, cylindrical, and spherical coordinates.

=== The Divergence in Various Coordinate Systems ===

All of Cartesian, cylindrical, and spherical coordinate systems will be modeled using the [[Calculus/Discrete vector calculus#Model of Curvilinear coordinates|discrete model of curvilinear coordinates]]. The same shorthand notation will be used.

Starting with an arbitrary vector field <math>\mathbf{F}(\Xi) = \sum_{c=1}^3 F_{\xi_c}(\Xi)\mathbf{v}_{\xi_c}(\Xi)</math>, the edge based function that approximates <math>\mathbf{F}</math> is <math>F_\text{approx}(e_{\xi_c,I}) = F_{\xi_c}((I+[0.5]_c)\Delta\Xi)</math> for all <math>c = 1,2,3</math> and <math>I \in \Z^3</math>.

For each <math>I \in \Z^3</math>, the divergence at node <math>n_I</math> is:

<math>\nabla \cdot F_\text{approx}|_{n_I} = \frac{1}{\tau_3(n_I)}\sum_{e \in \mathcal{E}} \left\{\begin{array}{cc} +\tau_2(e)F_\text{approx}(e) & (t_\bullet(e) = n_I \;\text{and}\; t^\bullet(e) \neq n_I) \\ -\tau_2(e)F_\text{approx}(e) & (t^\bullet(e) = n_I \;\text{and}\; t_\bullet(e) \neq n_I) \\ 0 & (\text{else}) \end{array}\right.</math>
<math> = \frac{1}{\prod_{c=1}^3 \Delta\xi_c|\mathbf{u}_{\xi_c}(I\Delta\Xi)|}\sum_{c=1}^3(\tau_2(e_{\xi_c,I})F_\text{approx}(e_{\xi_c,I}) - \tau_2(e_{\xi_c,I-[1]_c})F_\text{approx}(e_{\xi_c,I-[1]_c}))</math>
<math> = \frac{1}{\prod_{c=1}^3 \Delta\xi_c|\mathbf{u}_{\xi_c}(I\Delta\Xi)|}\sum_{c=1}^3\left((\prod_{c'\neq c}\Delta\xi_{c'}|\mathbf{u}_{\xi_{c'}}(I\Delta\Xi + [0.5\Delta\xi_c]_c)|)F_{\xi_c}(I\Delta\Xi + [0.5\Delta\xi_c]_c) - (\prod_{c'\neq c}\Delta\xi_{c'}|\mathbf{u}_{\xi_{c'}}(I\Delta\Xi - [0.5\Delta\xi_c]_c)|)F_{\xi_c}(I\Delta\Xi - [0.5\Delta\xi_c]_c)\right)</math>

Given an arbitrary coordinate <math>\Xi</math>, let <math>I(\Xi) \in \Z^3</math> index the node <math>n_{I(\Xi)}</math> whose position <math>I(\Xi)\Delta\Xi</math> is "closest" to <math>\Xi</math>.

The divergence of vector field <math>\mathbf{F}</math> at <math>\Xi</math> can now be approximated by the following:

<math>\nabla \cdot \mathbf{F}|_\Xi \approx \nabla \cdot F_\text{approx}|_{n_{I(\Xi)}} </math>

<math> = \frac{1}{\prod_{c=1}^3 \Delta\xi_c|\mathbf{u}_{\xi_c}(I(\Xi)\Delta\Xi)|}\sum_{c=1}^3\bigg((\prod_{c'\neq c}\Delta\xi_{c'}|\mathbf{u}_{\xi_{c'}}(I(\Xi)\Delta\Xi + [0.5\Delta\xi_c]_c)|)F_{\xi_c}(I(\Xi)\Delta\Xi + [0.5\Delta\xi_c]_c) </math>
<math> - (\prod_{c'\neq c}\Delta\xi_{c'}|\mathbf{u}_{\xi_{c'}}(I(\Xi)\Delta\Xi - [0.5\Delta\xi_c]_c)|)F_{\xi_c}(I(\Xi)\Delta\Xi - [0.5\Delta\xi_c]_c)\bigg)</math>

<math> \approx \frac{1}{\prod_{c=1}^3 \Delta\xi_c|\mathbf{u}_{\xi_c}(\Xi)|}\sum_{c=1}^3\left((\prod_{c'\neq c}\Delta\xi_{c'}|\mathbf{u}_{\xi_{c'}}(\Xi + [0.5\Delta\xi_c]_c)|)F_{\xi_c}(\Xi + [0.5\Delta\xi_c]_c) - (\prod_{c'\neq c}\Delta\xi_{c'}|\mathbf{u}_{\xi_{c'}}(\Xi - [0.5\Delta\xi_c]_c)|)F_{\xi_c}(\Xi - [0.5\Delta\xi_c]_c)\right)</math>
<math> = \frac{1}{\prod_{c=1}^3 |\mathbf{u}_{\xi_c}(\Xi)|}\sum_{c=1}^3\frac{1}{\Delta\xi_c}\left((\prod_{c'\neq c}|\mathbf{u}_{\xi_{c'}}(\Xi + [0.5\Delta\xi_c]_c)|)F_{\xi_c}(\Xi + [0.5\Delta\xi_c]_c) - (\prod_{c'\neq c}|\mathbf{u}_{\xi_{c'}}(\Xi - [0.5\Delta\xi_c]_c)|)F_{\xi_c}(\Xi - [0.5\Delta\xi_c]_c)\right)</math>

As <math>\Delta\xi_1, \Delta\xi_2, \Delta\xi_3 \to 0^+</math>, it becomes apparent that <math>\nabla \cdot \mathbf{F}|_\Xi = \frac{1}{\prod_{c=1}^3 |\mathbf{u}_{\xi_c}(\Xi)|}\sum_{c=1}^3\frac{\partial}{\partial\xi_c}((\prod_{c'\neq c}|\mathbf{u}_{\xi_{c'}}(\Xi)|)F_{\xi_c}(\Xi))\bigg|_\Xi</math>

==== The Divergence in Cartesian Coordinates ====

Applying the curvilinear coordinate model to Cartesian coordinates means that <math>\xi_1 = x</math>, <math>\xi_2 = y</math>, and <math>\xi_3 = z</math>. It is also the case that <math>\mathbf{u}_x(x,y,z) = \mathbf{i}</math>, <math>\mathbf{u}_y(x,y,z) = \mathbf{j}</math>, and <math>\mathbf{u}_z(x,y,z) = \mathbf{k}</math>. Applying the formula for the divergence that was derived from the discrete model of curvilinear coordinates to vector field <math>\mathbf{F} = F_x\mathbf{i} + F_y\mathbf{j} + F_z\mathbf{k}</math> gives:

<math>\nabla \cdot \mathbf{F}|_{(x,y,z)} = \frac{1}{1 \cdot 1 \cdot 1}(\frac{\partial}{\partial x}(1 \cdot 1 \cdot F_x)\bigg|_{(x,y,z)} + \frac{\partial}{\partial y}(1 \cdot 1 \cdot F_y)\bigg|_{(x,y,z)} + \frac{\partial}{\partial z}(1 \cdot 1 \cdot F_z)\bigg|_{(x,y,z)})</math>
<math> = \frac{\partial F_x}{\partial x} + \frac{\partial F_y}{\partial y} + \frac{\partial F_z}{\partial z}</math>

which is consistent with the formula for the divergence derived for continuous vector fields over Cartesian coordinates.

==== The Divergence in Cylindrical Coordinates ====

Applying the curvilinear coordinate model to cylindrical coordinates means that <math>\xi_1 = \rho</math>, <math>\xi_2 = \phi</math>, and <math>\xi_3 = z</math>. It is also the case that <math>\mathbf{u}_\rho(\rho,\phi,z) = \hat{\mathbf{\rho}}</math>, <math>\mathbf{u}_\phi(\rho,\phi,z) = \rho\hat{\mathbf{\phi}}</math>, and <math>\mathbf{u}_z(\rho,\phi,z) = \hat{\mathbf{z}}</math>. Applying the formula for the divergence that was derived from the discrete model of curvilinear coordinates to vector field <math>\mathbf{F} = F_\rho\hat{\mathbf{\rho}} + F_\phi\hat{\mathbf{\phi}} + F_z\hat{\mathbf{z}}</math> gives:

<math>\nabla \cdot \mathbf{F}|_{(\rho,\phi,z)} = \frac{1}{1 \cdot \rho \cdot 1}(\frac{\partial}{\partial\rho}(\rho \cdot 1 \cdot F_\rho)\bigg|_{(\rho,\phi,z)} + \frac{\partial}{\partial\phi}(1 \cdot 1 \cdot F_\phi)\bigg|_{(\rho,\phi,z)} + \frac{\partial}{\partial z}(1 \cdot \rho \cdot F_z)\bigg|_{(\rho,\phi,z)})</math>
<math> = \frac{1}{\rho}\frac{\partial}{\partial\rho}(\rho F_\rho) + \frac{1}{\rho}\frac{\partial F_\phi}{\partial\phi} + \frac{\partial F_z}{\partial z}</math>

which is consistent with the formula for the divergence derived for continuous vector fields over cylindrical coordinates.

==== The Divergence in Spherical Coordinates ====

Applying the curvilinear coordinate model to spherical coordinates means that <math>\xi_1 = r</math>, <math>\xi_2 = \theta</math>, and <math>\xi_3 = \phi</math>. It is also the case that <math>\mathbf{u}_r(r,\theta,\phi) = \hat{\mathbf{r}}</math>, <math>\mathbf{u}_\theta(r,\theta,\phi) = r\hat{\mathbf{\theta}}</math>, and <math>\mathbf{u}_\phi(r,\theta,\phi) = r\sin\theta\hat{\mathbf{\phi}}</math>. Applying the formula for the divergence that was derived from the discrete model of curvilinear coordinates to vector field <math>\mathbf{F} = F_r\hat{\mathbf{r}} + F_\theta\hat{\mathbf{\theta}} + F_\phi\hat{\mathbf{\phi}}</math> gives:

<math>\nabla \cdot \mathbf{F}|_{(r,\theta,\phi)} = \frac{1}{1 \cdot r \cdot r\sin\theta}(\frac{\partial}{\partial r}(r \cdot r\sin\theta \cdot F_r)\bigg|_{(r,\theta,\phi)} + \frac{\partial}{\partial\theta}(1 \cdot r\sin\theta \cdot F_\theta)\bigg|_{(r,\theta,\phi)} + \frac{\partial}{\partial\phi}(1 \cdot r \cdot F_\phi)\bigg|_{(r,\theta,\phi)})</math>
<math> = \frac{1}{r^2}\frac{\partial}{\partial r}(r^2 F_r) + \frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}(\sin\theta \cdot F_\theta) + \frac{1}{r\sin\theta}\frac{\partial F_\phi}{\partial\phi}</math>

which is consistent with the formula for the divergence derived for continuous vector fields over spherical coordinates.

== Discrete analog to the Curl ==

=== Loops ===

Consider an arbitrary path <math>P = \langle (e_1,s_1), (e_2,s_2), \dots, (e_k,s_k) \rangle</math> that begins and finishes on the same node. <math>P</math> is a closed loop, and the edge based function <math>\delta_1(e;P)</math> that denotes <math>P</math> satisfies <math>\nabla \cdot \delta_1(e;P)|_n = 0</math> for all nodes <math>n \in \mathcal{N}</math>. The divergence-free property of <math>\delta_1(e;P)</math> arises from the fact that closed loops have no starting or ending points and leave every node that they enter.

[[File:Multiloops.svg|thumb|A divergence-free edge based function is a superposition of closed loops.]]

[[File:Basis loops.svg|thumb|A directed graph with a set of basis loops highlighted in different colors.]]

A multi-loop is a superposition of loops, and any edge-based function that denotes a multi-loop is divergence-free. Moreover, any edge-based function that is divergence free can be expressed as a superposition of closed loops and hence denotes a multi-loop.

Given an arbitrary graph <math>\mathcal{G}</math>, a set of "basis" loops can be formed such that every multi-loop can be expressed as a '''unique''' linear superposition of the basis loops. A set of basis loops provides a means of expressing an arbitrary multi-loop as a list of coefficients. Each coefficient is a weight that is assigned to each basis loop. When the basis loops are multiplied by their respective weights and added together, the multi-loop that is represented by the list of coefficients is retrieved.

The set of all edge based functions is a [[Linear Algebra|vector space]] with dimension <math>|\mathcal{E}|</math>, and each edge based function is a vector whose entries are the list of values assigned to each edge. The set of all multi-loops is a subspace of the set of all edge based functions, and each multi-loop can be interpreted as a vector that is a linear combination of the loops from a set of [[Linear Algebra/Basis|basis]] loops.

Three approaches for generating a set <math>\mathcal{B}</math> of basis loops for an arbitrary graph <math>\mathcal{G}</math> will now be given:

To start, the following definitions are required:
*A '''path component''' of graph <math>\mathcal{G}</math> is an "island" of <math>\mathcal{G}</math> where any two nodes on the island can be linked by a continuous path, and any edge that is connected to a node in the path component is part of the same component and vice versa. 
*A '''bridge''' is an edge that when removed breaks its host path component into two path components. 
*An '''arc''', also referred to as a '''simple path''', is a path where with the exception of the start and end nodes does not visit or traverse the same node or edge more than once.
*A '''simple loop''' is a closed loop that does visit a node or edge twice on the same cycle. 
*A '''tree''' is a graph where any two nodes can be linked by exactly one simple path.
*A '''spanning tree''' is a sub-graph of a path component <math>\mathcal{G}</math> that is both a tree and which includes every node from <math>\mathcal{G}</math>.

'''Generating Basis Loops Approach 1''':

[[File:Generating basis loops 1.svg|thumb|A depiction of how approach 1 generates basis loops by repeatedly including arcs. The original loop is shown in red, and successive colors add each arc.]]

To start, all edges that are "bridges" should be removed. After all bridges have been removed, any nodes that are not connected to any edges and are themselves path components are themselves removed. Let <math>\mathcal{G}_1, \mathcal{G}_2, ..., \mathcal{G}_k</math> denote the resultant path components of <math>\mathcal{G}</math>. Let the set of basis loops <math>\mathcal{B}</math> start out empty. Let <math>N(\mathcal{B})</math> denote the set of nodes visited by the loops from <math>\mathcal{B}</math>. Let <math>E(\mathcal{B})</math> denote the set of edges that are traversed by loops from <math>\mathcal{B}</math>.

The following steps are now repeated for each path component <math>\mathcal{G}_i</math> for each <math>i = 1, 2, ..., k</math>:

Start by choosing an arbitrary simple loop <math>P</math> from <math>\mathcal{G}_i</math>. Add loop <math>P</math> to <math>\mathcal{B}</math>. For as long as there are nodes or edges in <math>\mathcal{G}_i</math> that are not contained by <math>N(\mathcal{B})</math> or <math>E(\mathcal{B})</math>, find an arc <math>P'</math> that starts and ends at nodes from <math>N(\mathcal{B})</math>. With the exception of the starting and ending nodes, arc <math>P'</math> must not visit a node or edge from <math>N(\mathcal{B})</math> or <math>E(\mathcal{B})</math>. The end nodes of <math>P'</math> can be linked by a path through the nodes and edges from <math>N(\mathcal{B})</math> and <math>E(\mathcal{B})</math> to form a simple loop <math>P''</math> that is then added to <math>\mathcal{B}</math>.

The above process forms a set <math>\mathcal{B}</math> of basis loops for <math>\mathcal{G}</math>.

Let <math>H</math> denote the number of path components of <math>\mathcal{G}</math>. As bridges and single nodes are removed from <math>\mathcal{G}</math>, the quantity <math>|\mathcal{E}|-|\mathcal{N}|+H</math> remains constant. After all bridges and single nodes have been removed, path component <math>\mathcal{G}_i</math> yields <math>E_i-N_i+1</math> basis loops where <math>N_i</math> and <math>E_i</math> are the number of nodes and edges from <math>\mathcal{G}_i</math> respectively. The total number of basis loops is therefore <math>|\mathcal{E}|-|\mathcal{N}|+H</math>. The vector space of multi-loops has <math>|\mathcal{E}|-|\mathcal{N}|+H</math> dimensions.

'''Generating Basis Loops Approach 2''':

[[File:Generating basis loops 2.svg|thumb|A depiction of how approach 2 generates basis loops by using a spanning tree. The solid edges are included in the spanning tree. Each dashed edge corresponds to a basis loop formed from the dashed edge and the unique path through the spanning tree.]]

Let <math>\mathcal{G}_1, \mathcal{G}_2, ..., \mathcal{G}_H</math> denote the path components of graph <math>\mathcal{G}</math>. <math>H</math> is the number of path components of <math>\mathcal{G}</math>. For each <math>i = 1, 2, ..., H</math>, let <math>\mathcal{G}_{T,i}</math> denote a "spanning tree" for path component <math>\mathcal{G}_i</math>. Each edge <math>e</math> from <math>\mathcal{G}_i</math> that is not part of the spanning tree <math>\mathcal{G}_{T,i}</math> corresponds to a basis loop. This basis loop is formed from edge <math>e</math> and the unique arc that connects the end points of <math>e</math> through the spanning tree <math>\mathcal{G}_{T,i}</math>.

For path component <math>\mathcal{G}_i</math>, let <math>N_i</math> and <math>E_i</math> respectively denote the number of nodes and edges from <math>\mathcal{G}_i</math>. The spanning tree <math>\mathcal{G}_{T,i}</math> contains <math>N_i-1</math> edges, so <math>E_i-N_i+1</math> edges from <math>\mathcal{G}_i</math> are not part of the spanning tree, and hence correspond to basis loops. The total number of basis loops from <math>\mathcal{G}</math> is <math>\sum_{i=1}^H (E_i-N_i+1) = |\mathcal{E}|-|\mathcal{N}|+H</math>.

'''Generating Basis Loops Approach 3''':

A third approach to generating a set of basis loops for graph <math>\mathcal{G}</math> proceeds by identifying the loops that pass through each node. To start, the set of nodes <math>\mathcal{N}</math> is ordered to form <math>\mathcal{N} = \{n_1,n_2,\dots,n_{|\mathcal{N}|}\}</math>.

For each node <math>n_i</math>, a set of basis loops is assigned to <math>n_i</math>. These basis loops are confined to the nodes <math>\{n_i,n_{i+1},\dots,n_{|\mathcal{N}|}\}</math> and are required to pass through <math>n_i</math>. These loops will be chosen such that no nonzero linear superposition of these loops will result in an edge based function where all edges that connect to <math>n_i</math> are assigned 0.

Each edge <math>e \in \mathcal{E}</math> has 2 endings: <math>t_\bullet(e)</math> and <math>t^\bullet(e)</math>. With respect to node <math>n_i</math>, let <math>T_i</math> denote the set of all edge endings that connect to <math>n_i</math>, excluding any edges that connect to nodes from <math>\{n_1,n_2,\dots,n_{i-1}\}</math>. The set <math>T_i</math> can now be partitioned into subsets according to the following equivalence relation: two edge endings <math>t_1, t_2 \in T_i</math> belong to the same partition if and only if there exists a simple loop that is confined to the nodes <math>\{n_i,n_{i+1},\dots,n_{|\mathcal{N}|}\}</math> and that visits node <math>n_i</math> exactly once, arriving through <math>t_1</math> or <math>t_2</math> and leaving through the other ending. Let <math>T_i = U_{i,1} \cup U_{i,2} \cup \dots U_{i,k_i}</math> denote the partitioning of <math>T_i</math>.

For each partition <math>U_{i,j}</math>, choose an ending <math>t_{i,j} \in U_{i,j}</math> that will serve as the "prime ending". Each other ending <math>t</math> from <math>U_{i,j}</math> will correspond to a basis loop that leaves <math>n_i</math> through <math>t</math> and returns to <math>n_i</math> through <math>t_{i,j}</math>, never visiting any nodes from <math>\{n_1,n_2,\dots,n_{i-1}\}</math>.

For each <math>i = 1, 2, ..., |\mathcal{N}|</math> and <math>j = 1, 2, ..., k_i</math>, each ending from <math>U_{i,j} \setminus \{t_{i,j}\}</math> corresponds to exactly 1 basis loop. Each edge <math>e \in \mathcal{E}</math> will contribute at most 1 ending to the total set <math>S = \bigcup_{i=1}^{|\mathcal{N}|}\bigcup_{j=1}^{k_i} (U_{i,j} \setminus \{t_{i,j}\})</math>, including edges that start and end on the same node. To count the number of basis loops, it is necessary to count the number <math>F</math> of edges that do not contribute any endings to <math>S</math>. These edges are the edges that contain a "prime ending" <math>t_{i,j}</math> and do not begin and end on the same node.

To determine <math>F</math>, let <math>\mathcal{G}_i</math> denote graph <math>\mathcal{G}</math> where the nodes <math>\{n_1,n_2,\dots,n_{i-1}\}</math> and any connected edges have been removed. <math>H_i</math> will denote the number of path components of <math>\mathcal{G}_i</math>. <math>H_1 = H</math> where <math>H</math> is the number of path components of <math>\mathcal{G}</math>, and <math>H_{|\mathcal{N}|} = 1</math>. Start with <math>\mathcal{G}_{i+1}</math> and reintroduce node <math>n_i</math> and any connected edges to restore <math>\mathcal{G}_i</math>. Introducing <math>n_i</math> creates a new path component. Connecting <math>n_i</math> to a disconnected path component via an edge will reduce the number of path components by 1, but will also create an edge that contains a prime ending. Connecting <math>n_i</math> to an already connected path component via an edge will not change the number of path components, nor create an edge with a prime ending. The number of prime endings that are connected to node <math>n_i</math> not counting edges that start and end on <math>n_i</math> is <math>H_{i+1}-H_i+1</math> (the +1 compensates for the path component that <math>n_i</math> introduces by default).

<math> F = \sum_{i=1}^{|\mathcal{N}|-1} (H_{i+1}-H_i+1) = (H_{|\mathcal{N}|}-H_1) + (|\mathcal{N}|-1) = (1-H)+(|\mathcal{N}|-1) = |\mathcal{N}|-H</math>

so the number of basis loops is <math>|\mathcal{E}|-F = |\mathcal{E}|-|\mathcal{N}|+H</math> which is in agreement with the other approaches.

=== Surfaces ===

[[File:Closed surface and closed path.svg|thumb|Given an arbitrary oriented closed surface and an arbitrary oriented closed path, this image demonstrates how the number of times the path leaves the surface is equal to the number of times the path enters the surface.]]

Consider an arbitrary surface <math>S = \{(e_1,s_1),(e_2,s_2),...,(e_k,s_k)\}</math>. Surface <math>S</math> is "closed" if it has the following property: any path <math>P</math> that is also closed must cross surface <math>S</math> in the preferred direction and against the preferred direction an equal number of times. In mathematical terms, this means that given any closed path <math>P</math> that <math>\Sigma_{\mathcal{E}}(P,S) = \int_{e \in P} \delta_2(e;S) = \int_{e \in S} \delta_1(e;P) = \sum_{e \in \mathcal{E}} \delta_1(e;P)\delta_2(e;S)\tau_2(e)\tau_1(e) = 0</math>.

[[File:Multi closed surfaces.svg|thumb|Any conservative edge based function (denoted by arrows in this image) is the superposition of closed surfaces (denoted by the dashed lines).]]

An edge based function <math>f</math> is "conservative" if any path integral of <math>f</math> along an arbitrary '''closed''' path <math>P</math> is 0: <math>\int_{e \in P} f(e) = 0</math>. In a similar manner to how any divergence free edge based function denotes a multi-loop, any conservative edge based function denotes a "multi-closed surface", which is the superposition of any number of closed surfaces.

[[File:Graph Surface Model.svg|thumb|An example of multiple surfaces sharing the same boundary.]]

A surface that is not closed is said to have a "boundary". Non-closed surfaces that share a common boundary can be grouped into "equivalence classes" and treated as though they are identical. It is relatively simple to observe that if surfaces <math>S_1</math> and <math>S_2</math> share a common boundary then the surface denoted by the edge based function <math>\delta_2(e;S_1)-\delta_2(e;S_2)</math> is a closed surface since surfaces <math>S_1</math> and <math>-S_2</math> form the two halves of a "clam shell" that when added form a closed surface (<math>-S_2</math> is <math>S_2</math> with the orientation reversed). The edge based function <math>\delta_2(e;S_1)-\delta_2(e;S_2)</math> is therefore conservative. The difference between edge based functions that denote multi-surfaces with a common boundary denotes a multi-closed surface, and is hence conservative. If edge based function <math>f</math> denotes an arbitrary multi-surface, then when given any conservative edge based function <math>g</math>, the edge based function <math>f+g</math> denotes a multi-surface with the same boundary as <math>f</math>.

Given an arbitrary edge based function <math>f</math> that denotes a multi-surface, <math>\partial f</math> will denote the boundary of <math>f</math>, and <math>\mathfrak{B}[f]</math> will denote the set of all edge based functions that have the same boundary as <math>f</math>. For any edge based function <math>g</math>, <math>\partial g = \partial f</math> if and only if <math>g \in \mathfrak{B}[f]</math>. Moreover, <math>\partial g = \partial f</math> if and only if <math>g - f</math> is conservative.

Multi-surfaces, and their corresponding edge based functions, that all share a common boundary will be grouped into a single "equivalence class" and treated as if they are identical. An edge based function <math>f</math> "represents" its equivalence class <math>\mathfrak{B}[f]</math>. An equivalence class <math>\mathfrak{B}</math> effectively denotes the common boundary shared by all surfaces <math>f \in \mathfrak{B}</math>, and will be referred to as a "boundary".

The property of an edge based function being conservative is linear, which means that given conservative edge based functions <math>f</math> and <math>g</math> and coefficients <math>c</math> and <math>d</math>, it will be the case that the linear superposition <math>cf + dg</math> is also conservative. The linearity of being conservative implies that: given edge based functions <math>f_1,f_2</math> and <math>g_1, g_2</math> where <math>\partial f_1 = \partial f_2</math> and <math>\partial g_1 = \partial g_2</math>, then <math>\partial(cf_1 + dg_1) = \partial(cf_2 + dg_2)</math> for any coefficients <math>c, d</math>. The previous property implies that the linear superposition of boundaries <math>\mathfrak{B}_1</math> and <math>\mathfrak{B}_2</math> can be performed unambiguously by choosing "representative" surfaces <math>f \in \mathfrak{B}_1</math> and <math>g \in \mathfrak{B}_2</math> and generating a surface <math>cf + dg</math> that represents the linear superposition <math>c\mathfrak{B}_1 + d\mathfrak{B}_2</math>.

The set of all conservative edge based functions (the set of all closed surfaces) is a vector space, and the set of all equivalence classes/boundaries is a vector space.

A set of basis closed surfaces can be generated by the following simple process: Let <math>\mathcal{G}_1, \mathcal{G}_2, ..., \mathcal{G}_H</math> denote the path components of <math>\mathcal{G}</math>. For each <math>i = 1, 2, ..., H</math>, a set of basis closed surfaces for path component <math>\mathcal{G}_i</math> can be produced by wrapping each node except for 1 in a basis surface as shown below. Any "bubble" can be generated by adding the bubbles that individually wrap each node inside the bubble. A surface that only wraps the remaining node can be generated by adding all of the bubbles that wrap the other nodes in the path component, and then reversing the orientation. The number of basis closed surfaces is <math>N_i - 1</math> where <math>N_i</math> is the number of nodes from path component <math>\mathcal{G}_i</math>. The total number of basis closed surfaces for graph <math>\mathcal{G}</math> is <math>\sum_{i=1}^H (N_i - 1) = |\mathcal{N}| - H</math>.

A set of basis equivalence classes/boundaries can be generated by the following simple process: Let <math>\mathcal{G}_1, \mathcal{G}_2, ..., \mathcal{G}_H</math> denote the path components of <math>\mathcal{G}</math>. For each <math>i = 1,2,...,H</math>, let <math>\mathcal{G}_{T,i}</math> denote a "spanning tree" for path component <math>\mathcal{G}_i</math>. Each edge <math>e</math> from <math>\mathcal{G}_i</math> that is not part of the spanning tree <math>\mathcal{G}_{T,i}</math> corresponds to a basis equivalence class which consists of all edge based functions that can be formed by adding a conservative edge based function to the edge based function that is <math>1</math> only at edge <math>e</math> and <math>0</math> for all other edges. The number of basis equivalence classes for path component <math>\mathcal{G}_i</math> is <math>E_i-(N_i-1)</math> where <math>N_i</math> and <math>E_i</math> are the number of nodes and edges from path component <math>\mathcal{G}_i</math> respectively. <math>N_i-1</math> is the number of edges from spanning tree <math>\mathcal{G}_{T,i}</math>. The total number of basis equivalence classes is <math>\sum_{i=1}^H (E_i-N_i+1) = |\mathcal{E}|-|\mathcal{N}|+H</math>. This is the same as the number of basis closed loops.

The fact that the previously mentioned set of basis of equivalence classes is in fact a basis, which means that every equivalence class is a unique linear superposition of the basis equivalence classes, can be confirmed by the following observations: Let <math>f</math> be an arbitrary edge based function. A conservative edge based function <math>g</math> can be created by letting <math>g = f</math> for all edges that are in the spanning trees, and choosing values for <math>g</math> at the edges that are not in the spanning trees such that <math>g</math> is conservative. The difference <math>f-g</math> is 0 for all edges in the spanning trees. The coefficients for each basis equivalence class are the values that <math>f-g</math> attains at each non spanning tree edge. Moreover, any non-zero choice of coefficients for the basis equivalence classes will always create a non-conservative edge based function.

In the images below, the left image shows a set of basis closed surfaces, the middle image shows a set of basis non-closed surfaces, and the right image shows how any edge based function/surface (shown in blue) is the linear combination of basis non-closed surfaces plus a closed surface (shown in black).

{{gallery
|width = 300
|File:Generating basis closed surfaces.svg|A set of basis closed surfaces for the shown graph. Any closed surface is a linear combination of the basis surfaces. A surface that wraps the bottom right node can be formed by adding all of the shown basis surfaces, and then reversing the orientation.
|File:Generating basis non closed surfaces.svg|A set of basis open surfaces for the shown graph. The edges in the spanning tree do not correspond to any basis surface. Any surface is a unique linear combination of the shown open surfaces plus a closed surface.
|File:Linear decomposition of a non closed surface.svg|A non-closed surface, shown in blue, is the linear superposition of basis non-closed surfaces and a closed surface. Observe how the blue arrows on the edges are the sum of the black arrows (opposites cancel out).
}}

=== Loop-Surface Duality ===

The number of dimensions of the set of multi-loops is <math>|\mathcal{E}|-|\mathcal{N}|+H</math>. The number of dimensions of the set of surface equivalence classes is also <math>|\mathcal{E}|-|\mathcal{N}|+H</math>. This equivalence implies that it is possible to form a '''[[w:Bijection|bijective]]''' (one to one and onto) '''[[w:Linear map|linear mapping]]''' between the set of multi-loops and surface equivalence classes. In simpler terms, any multi-loop corresponds to a unique equivalence class of surfaces, and any equivalence class of surfaces corresponds to a unique multi-loop. The correspondence is linear, which means that given multi-loops <math>f_1</math> and <math>f_2</math> which respectively correspond to equivalence classes <math>\mathfrak{B}_1</math> and <math>\mathfrak{B}_2</math>, and coefficients <math>c</math> and <math>d</math>, then multi-loop <math>cf_1 + df_2</math> corresponds to equivalence class <math>c\mathfrak{B}_1 + d\mathfrak{B}_2</math>.

Previously, when given a multi-surface denoted by edge based function <math>g</math>, the boundary of <math>g</math> was a hypothetical boundary assumed to be shared by all edge based functions whose difference with <math>g</math> was a conservative edge based function. The boundary <math>\partial g</math> is now an edge based function that denotes the "counter-clockwise oriented boundary of <math>g</math>". Since the boundary of any surface is a loop, it will be required that <math>\partial g</math> is a multi-loop: <math>\nabla \cdot (\partial g) = 0</math>. Mapping each surface to its boundary should also be linear: given surfaces <math>g_1, g_2</math> and coefficients <math>c, d</math>, then <math>\partial(cg_1+dg_2) = c\partial g_1 + d\partial g_2</math>.

The bijective linear mapping between the set of multi-loops and surface equivalence classes maps each multi-loop <math>f</math> to the equivalence class of surfaces which have <math>f</math> as their common boundary.

[[File:Path surface intersections.svg|thumb|A path is shown intersecting a surface multiple times. The boundary of the surface has a counter-clockwise orientation. The net number of times that the surface is crossed in the preferred direction is N = +1-1+1+1+1 = 3]]

Given path <math>P = \langle (e_{P,1},s_{P,1}), (e_{P,2},s_{P,2}), ..., (e_{P,k_P},s_{P,k_P}) \rangle</math> and surface <math>S = \{ (e_{S,1},s_{S,1}), (e_{S,2},s_{S,2}), ..., (e_{S,k_S},s_{S,k_S}) \}</math>, the net number of times that path <math>P</math> intersects surface <math>S</math> in the preferred direction is <math>\Sigma_\mathcal{E}(P,S) = \sum_{i=1}^{k_P}\sum_{j=1}^{k_S} \left\{\begin{array}{cc} s_{P,i}s_{S,j} & (e_{P,i} = e_{S,j}) \\ 0 & (e_{P,i} \neq e_{S,j}) \\ \end{array}\right. = \sum_{e \in \mathcal{E}} \delta_1(e;P)\delta_2(e;S)\tau_2(e)\tau_1(e) = \int_{e \in P} \delta_2(e;S) = \int_{e \in S} \delta_1(e;P)</math>. This formula generalizes to edge based functions that denote multi-paths and multi-surfaces. Given edge based function <math>f</math> that denotes a multi-path and edge based function <math>g</math> that denotes a multi-surface, then the net number of times <math>f</math> intersects <math>g</math> in the preferred direction is <math>\Sigma_\mathcal{E}(f,g) = \sum_{e \in \mathcal{E}} f(e)g(e)\tau_2(e)\tau_1(e) = \int_{e \in f} g(e) = \int_{e \in g} f(e)</math>.

It is also important to note that <math>\Sigma_\mathcal{E}(f,g)</math> is linear with respect to <math>f</math> and <math>g</math>: <math>\Sigma_\mathcal{E}(cf_1+df_2,g) = c\Sigma_\mathcal{E}(f_1,g) + d\Sigma_\mathcal{E}(f_2,g)</math> and <math>\Sigma_\mathcal{E}(f,cg_1+dg_2) = c\Sigma_\mathcal{E}(f,g_1)+d\Sigma_\mathcal{E}(f,g_2)</math>. It is also important to note that if <math>f</math> is a multi-loop, and <math>g</math> is a multi-closed surface (conservative), then <math>\Sigma_\mathcal{E}(f,g) = 0</math>. This latter property implies that given surfaces <math>g_1, g_2</math> with a common boundary (<math>\partial g_1 = \partial g_2</math>), then <math>\Sigma_\mathcal{E}(f,g_1) = \Sigma_\mathcal{E}(f,g_2)</math> since <math>g_2-g_1</math> is conservative.

[[File:Loop surface duality.svg|thumb|Two surfaces, each with a counter-clockwise oriented boundary, are shown. The net number of times each boundary intersects the other surface is the same. The red boundary passes through the green surface in the preferred direction 2 times, and the green boundary passes through the red surface in the preferred direction 2 times.]]

One requirement of the mapping of a boundary to each surface equivalence class is that given multi-surfaces <math>g_1</math> and <math>g_2</math>, the net number of times the counter clockwise oriented boundary <math>\partial g_1</math> passes through <math>g_2</math> in the preferred direction is identical to the net number of times <math>\partial g_2</math> passes through <math>g_1</math> in the preferred direction: <math>\sum_{e \in \mathcal{E}} (\partial g_1)|_e g_2(e) \tau_2(e)\tau_1(e) = \Sigma_\mathcal{E}(\partial g_1,g_2) = \Sigma_\mathcal{E}(\partial g_2,g_1) = \sum_{e \in \mathcal{E}} (\partial g_2)|_e g_1(e) \tau_2(e)\tau_1(e)</math>. See the 2nd image on the right for an example of this property.

Since <math>\Sigma_\mathcal{E}(f,g)</math> is linear with respect to <math>f</math> and <math>g</math>, and <math>\partial g</math> is linear with respect to <math>g</math>, the condition <math>\Sigma_\mathcal{E}(\partial g_1,g_2) = \Sigma_\mathcal{E}(\partial g_2,g_1)</math> is linear with respect to <math>g_1</math> and <math>g_2</math>:
: <math>\Sigma_\mathcal{E}(\partial g_{1,1},g_2) = \Sigma_\mathcal{E}(\partial g_2,g_{1,1}) \;\text{and}\; \Sigma_\mathcal{E}(\partial g_{1,2},g_2) = \Sigma_\mathcal{E}(\partial g_2,g_{1,2}) \implies \Sigma_\mathcal{E}(\partial(cg_{1,1}+dg_{1,2}),g_2) = \Sigma_\mathcal{E}(\partial g_2,cg_{1,1}+dg_{1,2})</math>
: <math>\Sigma_\mathcal{E}(\partial g_1,g_{2,1}) = \Sigma_\mathcal{E}(\partial g_{2,1},g_1) \;\text{and}\; \Sigma_\mathcal{E}(\partial g_1,g_{2,2}) = \Sigma_\mathcal{E}(\partial g_{2,2},g_1) \implies \Sigma_\mathcal{E}(\partial g_1,cg_{2,1}+dg_{2,2}) = \Sigma_\mathcal{E}(\partial(cg_{2,1}+dg_{2,2}),g_1)</math>

Since the condition <math>\Sigma_\mathcal{E}(\partial g_1,g_2) = \Sigma_\mathcal{E}(\partial g_2,g_1)</math> is linear, mapping each multi-loop to the surface equivalence class that have the multi-loop as their common boundary is done by determining a set of basis multi-loops <math>\{f_1,f_2,...,f_{|\mathcal{E}|-|\mathcal{N}|+H}\}</math> and a set of basis surface equivalence classes <math>\{\mathfrak{B}[g_1],\mathfrak{B}[g_2],...,\mathfrak{B}[g_{|\mathcal{E}|-|\mathcal{N}|+H}]\}</math> such that for all <math>i,j \in \{1, 2, ..., |\mathcal{E}|-|\mathcal{N}|+H\}</math>, it is the case that <math>\Sigma_\mathcal{E}(f_i,g_j) = \Sigma_\mathcal{E}(f_j,g_i)</math>. For each <math>i = 1, 2, ..., |\mathcal{E}|-|\mathcal{N}|+H</math>, <math>\mathfrak{B}[g_i]</math> corresponds to <math>f_i</math> and the rest of the mapping is determined via linearity. Given an arbitrary multi-surface <math>g</math>, expressing the equivalence class <math>\mathfrak{B}[g]</math> as a linear combination of basis surface equivalence classes: <math>\mathfrak{B}[g] = \sum_{i=1}^{|\mathcal{E}|-|\mathcal{V}|+H} c_i\mathfrak{B}[g_i]</math> means that <math>\partial g = \sum_{i=1}^{|\mathcal{E}|-|\mathcal{V}|+H} c_if_i</math>.

One such choice of basis multi-loops and basis surface equivalence classes can be determined as follows: Let <math>\mathcal{G}_1,\mathcal{G}_2,...,\mathcal{G}_H</math> denote the path components of <math>\mathcal{G}</math>. For each <math>i = 1,2,...,H</math>, let <math>\mathcal{G}_{T,i}</math> denote a "spanning tree" for path component <math>\mathcal{G}_i</math>. Each edge <math>e</math> that is not part of any spanning tree corresponds to both a basis loop and a basis surface equivalence class. The basis loop <math>f_e</math> is the loop that passes through edge <math>e</math> and the unique path through a spanning tree that links the endpoints of <math>e</math>. A surface <math>g_e</math> that is part of the basis surface equivalence class is the surface that contains edge <math>e</math> alone. It is easy to check that <math>\Sigma_\mathcal{E}(f_e,g_e) = 1</math>, and when given different edges <math>e</math> and <math>e'</math>, that <math>\Sigma_\mathcal{E}(f_e,g_{e'}) = \Sigma_\mathcal{E}(f_{e'},g_e) = 0</math>. This is only one example of a pair of a set of basis multi-loops and a set of basis surface equivalence classes (this example also does not make sense, as the prescribed boundary of each basis surface passes through the surface itself).

==== Example ====

This example will demonstrate the duality of loops and surfaces using an example from electrical engineering: two simple circuits that are magnetically linked to themselves and to each other.

Three images describing the scenario are shown below:
:The left image depicts a circuit diagram that denotes the circuits that are the subject of this example. Each loop consists of a voltage source in series with a resistor, and a self inductance followed by a mutual inductance. At the bottom of the image, the "magnetic circuit" is shown. The magnetic circuit carries magnetic flux instead of electric current, and helps to quantify the self and mutual inductance between the two loops of current. 
:The center image depicts the orientation of the magnetic circuit relative to the electric circuit. As can be seen, the magnetic flux through the first loop is measured relative to "down", while the magnetic flux through the second loop is measured relative to "up".   
:The right image depicts the directed graph that model the electric circuit. The set of nodes is <math>\mathcal{N} = \{n_1,n_2,n_3,n_4,n_5,n_6\}</math> and the set of edges is <math>\mathcal{E} = \{e_1,e_2,e_3,e_4,e_5,e_6,e_7\}</math>. 
:: The origin and terminal nodes of each edge are: <math>(t_\bullet(e_1),t^\bullet(e_1)) = (n_1,n_2)</math>; <math>(t_\bullet(e_2),t^\bullet(e_2)) = (n_2,n_1)</math>; <math>(t_\bullet(e_3),t^\bullet(e_3)) = (n_3,n_4)</math>; <math>(t_\bullet(e_4),t^\bullet(e_4)) = (n_4,n_3)</math>; <math>(t_\bullet(e_5),t^\bullet(e_5)) = (n_5,n_6)</math>; <math>(t_\bullet(e_6),t^\bullet(e_6)) = (n_6,n_5)</math>; and <math>(t_\bullet(e_7),t^\bullet(e_7)) = (n_6,n_5)</math>.
:: There are 4 basis loops: <math>P_1 = \langle (e_1,+1), (e_2,+1) \rangle</math>; <math>P_2 = \langle (e_3,+1), (e_4,+1) \rangle</math>; <math>P_3 = \langle (e_5,+1), (e_7,+1) \rangle</math>; and <math>P_4 = \langle (e_7,-1), (e_6,+1) \rangle</math>. These loops are the respective counterclockwise boundaries of the following 4 surfaces: <math>S_1 = \{(e_5,+1)\}</math>; <math>S_2 = \{(e_6,+1)\}</math>; <math>S_3 = \{(e_2,+1)\}</math>; and <math>S_4 = \{(e_4,+1)\}</math>. It can be determined that the paths intersect the surfaces the following numbers of times:

{| class = "wikitable"
|- 
| <math>\Sigma_\mathcal{E}(P_1,S_1) = 0</math> || <math>\Sigma_\mathcal{E}(P_1,S_2) = 0</math> || <math>\Sigma_\mathcal{E}(P_1,S_3) = 1</math> || <math>\Sigma_\mathcal{E}(P_1,S_4) = 0</math> 
|- 
| <math>\Sigma_\mathcal{E}(P_2,S_1) = 0</math> || <math>\Sigma_\mathcal{E}(P_2,S_2) = 0</math> || <math>\Sigma_\mathcal{E}(P_2,S_3) = 0</math> || <math>\Sigma_\mathcal{E}(P_2,S_4) = 1</math> 
|- 
| <math>\Sigma_\mathcal{E}(P_3,S_1) = 1</math> || <math>\Sigma_\mathcal{E}(P_3,S_2) = 0</math> || <math>\Sigma_\mathcal{E}(P_3,S_3) = 0</math> || <math>\Sigma_\mathcal{E}(P_3,S_4) = 0</math> 
|- 
| <math>\Sigma_\mathcal{E}(P_4,S_1) = 0</math> || <math>\Sigma_\mathcal{E}(P_4,S_2) = 1</math> || <math>\Sigma_\mathcal{E}(P_4,S_3) = 0</math> || <math>\Sigma_\mathcal{E}(P_4,S_4) = 0</math> 
|}

Note that the above table is symmetric along its diagonal.

{|
| [[File:Electric and magnetic circuit.svg|thumb|300px|An electric and magnetic circuit that will be subject of this example.]]
| [[File:Electric and magnetic circuit realtive positions.svg|thumb|300px|The magnetic circuit is upright relative to the electric circuit. "Up" in the magnetic circuit protrudes from the page in the electric circuit.]]
| [[File:Electric and magnetic circuit modelling graph.svg|thumb|300px|A directed graph that will be used to model the electric and magnetic circuits, along with chosen basis loops and their corresponding surfaces.]]
|}

The goal will be to derive expressions for the self inductances <math>L_1, L_2</math>, and the mutual inductance <math>M_{1,2}</math>.

An edge based function <math>B</math> will denote the magnetic field density at each edge in the edge's preferred orientation.

A total current of <math>I_1</math> flows around loop <math>P_1</math>. This current flows outwards through surface <math>S_3</math>, so due to Ampere's law, the total path integral of <math>B</math> around <math>P_3</math> is <math>\int_{e \in P_3} B(e) = \mu_0I_1</math> (<math>\mu_0 = 4\pi \times 10^{-7}\text{Tm/A}</math> is the magnetic permeability of free space). Therefore: <math>\tau_1(e_5)B(e_5) + \tau_1(e_7)B(e_7) = \mu_0I_1</math>

A total current of <math>I_2</math> flows around loop <math>P_2</math>. This current flows outwards through surface <math>S_4</math>, so due to Ampere's law, the total path integral of <math>B</math> around <math>P_4</math> is <math>\int_{e \in P_4} B(e) = \mu_0I_2</math>. Therefore: <math>\tau_1(e_6)B(e_6) - \tau_1(e_7)B(e_7) = \mu_0I_2</math>

The divergence of the magnetic field at any node is 0: <math>\nabla \cdot B = 0</math>. Therefore: <math>\tau_2(e_5)B(e_5)-\tau_2(e_6)B(e_6)-\tau_2(e_7)B(e_7) = 0</math>

The following equations describe the magnetic field which is present in the magnetic circuit:

<math>\tau_1(e_5)B(e_5) + \tau_1(e_7)B(e_7) = \mu_0I_1</math>

<math>\tau_1(e_6)B(e_6) - \tau_1(e_7)B(e_7) = \mu_0I_2</math>

<math>\tau_2(e_5)B(e_5)-\tau_2(e_6)B(e_6)-\tau_2(e_7)B(e_7) = 0</math>

Solving the above equations yields:

<math>B(e_5) = \mu_0\frac{(\tau_2(e_6)\tau_1(e_7)+\tau_1(e_6)\tau_2(e_7))I_1 + \tau_2(e_6)\tau_1(e_7)I_2}{\tau_2(e_5)\tau_1(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_2(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_1(e_6)\tau_2(e_7)}</math>

<math>B(e_6) = \mu_0\frac{\tau_2(e_5)\tau_1(e_7)I_1 + (\tau_2(e_5)\tau_1(e_7)+\tau_1(e_5)\tau_2(e_7))I_2}{\tau_2(e_5)\tau_1(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_2(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_1(e_6)\tau_2(e_7)}</math>

<math>B(e_7) = \mu_0\frac{\tau_2(e_5)\tau_1(e_6)I_1 - \tau_1(e_5)\tau_2(e_6)I_2}{\tau_2(e_5)\tau_1(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_2(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_1(e_6)\tau_2(e_7)}</math>

The magnetic flux through the first loop (surface <math>S_1</math>) is:

<math>\Phi_1 = \tau_2(e_5)B(e_5) = \mu_0\frac{(\tau_2(e_5)\tau_2(e_6)\tau_1(e_7)+\tau_2(e_5)\tau_1(e_6)\tau_2(e_7))I_1 + \tau_2(e_5)\tau_2(e_6)\tau_1(e_7)I_2}{\tau_2(e_5)\tau_1(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_2(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_1(e_6)\tau_2(e_7)}</math>

The magnetic flux through the second loop (surface <math>S_2</math>) is:

<math>\Phi_2 = \tau_2(e_6)B(e_6) = \mu_0\frac{\tau_2(e_5)\tau_2(e_6)\tau_1(e_7)I_1 + (\tau_2(e_5)\tau_2(e_6)\tau_1(e_7)+\tau_1(e_5)\tau_2(e_6)\tau_2(e_7))I_2}{\tau_2(e_5)\tau_1(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_2(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_1(e_6)\tau_2(e_7)}</math>

This means that the self inductances of the first and second loop are respectively:

<math>L_1 = \mu_0\frac{\tau_2(e_5)\tau_2(e_6)\tau_1(e_7)+\tau_2(e_5)\tau_1(e_6)\tau_2(e_7)}{\tau_2(e_5)\tau_1(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_2(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_1(e_6)\tau_2(e_7)}</math> (the coefficient of <math>I_1</math> in the expression for <math>\Phi_1</math>)

and

<math>L_2 = \mu_0\frac{\tau_2(e_5)\tau_2(e_6)\tau_1(e_7)+\tau_1(e_5)\tau_2(e_6)\tau_2(e_7)}{\tau_2(e_5)\tau_1(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_2(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_1(e_6)\tau_2(e_7)}</math> (the coefficient of <math>I_2</math> in the expression for <math>\Phi_2</math>)

The mutual inductance is:

<math>M_{1,2} = \mu_0\frac{\tau_2(e_5)\tau_2(e_6)\tau_1(e_7)}{\tau_2(e_5)\tau_1(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_2(e_6)\tau_1(e_7) + \tau_1(e_5)\tau_1(e_6)\tau_2(e_7)}</math> (the coefficient of <math>I_2</math> in the expression for <math>\Phi_1</math>, or the coefficient of <math>I_1</math> in the expression for <math>\Phi_2</math>)

=== The Curl ===

==== Circulation Density ====

Given an arbitrary edge based '''rate-of-gain''' function <math>f</math>, and an arbitrary multi-loop <math>P</math> (<math>\nabla \cdot P|_n = 0</math>), the "'''circulation'''" of <math>f</math> around <math>P</math> is <math>\int_{e \in P} f(e) = \sum_{e \in \mathcal{E}} f(e)P(e)\tau_2(e)\tau_1(e)</math>. With <math>f</math> fixed, the circulation is a linear function of <math>P</math>: <math>\int_{e \in cP_1 + dP_2} f(e) = c\int_{e \in P_1} f(e) + d\int_{e \in P_2} f(e)</math>.

[[File:2D loop decomposition.svg|thumb|A large loop can be decomposed into an ensemble of tiny loops.]]

[[File:3D loop decomposition.svg|thumb|A large loop can be decomposed into an ensemble of tiny loops.]]

As seen in the images to the right, a large loop can be decomposed into a sum of smaller loops, and since the circulation is linear with respect to the loop that is being followed, the total circulation around the large loop is simply the sum of the circulations around the little loops.

Each edge <math>e \in \mathcal{E}</math> is its own little surface: <math>\{(e,+1)\}</math>, and any surface is a linear combination of these basis surfaces. Any surface <math>S</math> can be expressed as the following linear combination: <math>\forall e \in \mathcal{E} : S(e) = \sum_{e' \in \mathcal{E}} S(e')\tau_1(e')\delta_2(e;\{(e',+1)\})</math>. The length term <math>\tau_1(e')</math> cancels out the factor of <math>\frac{1}{\tau_1(e')}</math> in <math>\delta_2(e;\{(e',+1)\})</math>. If the circulation of <math>f</math> around the counter clockwise boundary of each basis surface <math>\{(e',+1)\}</math> is computed, then the total circulation around the counter clockwise boundary of <math>S</math> can be computed via the same linear combination:

<math>\int_{e \in \partial S} f(e) = \sum_{e' \in \mathcal{E}} S(e')\tau_1(e')\int_{e \in \partial(\{(e',+1)\})} f(e)</math>

The circulation per unit area of surface <math>\{(e',+1)\}</math> is <math>\frac{1}{\tau_2(e')}\int_{e \in \partial(\{(e',+1)\})} f(e)</math>. This "circulation density" at edge <math>e'</math> will be referred to as the "curl" and denoted by the edge based function <math>\nabla \times f|_{e'} = \frac{1}{\tau_2(e')}\int_{e \in \partial(\{(e',+1)\})} f(e)</math>.

The total circulation around the counter clockwise boundary of <math>S</math> is <math>\int_{e \in \partial S} f(e) = \sum_{e' \in \mathcal{E}} S(e')\tau_1(e')\int_{e \in \partial(\{(e',+1)\})} f(e) = \sum_{e' \in \mathcal{E}} S(e')\tau_1(e')\tau_2(e')(\nabla \times f|_{e'}) = \int_{e \in S} \nabla \times f|_e</math>. This is the discrete analog to Stokes' Theorem.

==== Computing the counter-clockwise boundary of a surface ====

[[File:Curl as the boundary.svg|thumb|The counter-clockwise boundary of the red surface is generated by the curl around each of the blue loops that wrap around the boundary.]]

[[File:Curl as the boundary 2.svg|thumb|Given an oriented surface, a cross-section of which is shown here, the circulation around a tiny loop is only non-zero at the boundaries. The curl is oriented in the direction of the counter-clockwise boundary.]]

Let <math>S</math> be an arbitrary multi-surface. Since <math>S</math> can be expressed as the linear combination <math>\forall e \in \mathcal{E} : S(e) = \sum_{e' \in \mathcal{E}} S(e')\tau_1(e')\delta_2(e;\{(e',+1)\})</math>, the counterclockwise boundary <math>\partial S</math> is the linear combination of boundaries: <math>\forall e \in \mathcal{E} : \partial S|_e = \sum_{e' \in \mathcal{E}} S(e')\tau_1(e')(\partial(\{(e',+1)\})|_e)</math>.

Given arbitrary edges <math>e_1, e_2 \in \mathcal{E}</math>, it is the case that <math>\Sigma_{\mathcal{E}}(\partial(\{(e_1,+1)\}),\{(e_2,+1)\}) = \Sigma_{\mathcal{E}}(\partial(\{(e_2,+1)\}),\{(e_1,+1)\})</math> due to the requirement placed on the boundaries of surfaces. This restriction can also be expressed as: <math>\partial(\{(e_1,+1)\})|_{e_2}\tau_2(e_2) = \partial(\{(e_2,+1)\})|_{e_1}\tau_2(e_1)</math>.

The linear combination that describes <math>\partial S</math> becomes:

<math>\forall e \in \mathcal{E} : \partial S|_e = \sum_{e' \in \mathcal{E}} S(e')\tau_1(e')(\partial(\{(e',+1)\})|_e)</math>
<math> = \sum_{e' \in \mathcal{E}} S(e')\frac{\tau_1(e')}{\tau_2(e)}(\partial(\{(e',+1)\})|_e\tau_2(e))</math>
<math> = \sum_{e' \in \mathcal{E}} S(e')\frac{\tau_1(e')}{\tau_2(e)}(\partial(\{(e,+1)\})|_{e'}\tau_2(e'))</math>
<math> = \frac{1}{\tau_2(e)}\sum_{e' \in \mathcal{E}} S(e')\partial(\{(e,+1)\})|_{e'}\tau_1(e')\tau_2(e')</math> 
<math> = \frac{1}{\tau_2(e)}\Sigma_{\mathcal{E}}(\partial(\{(e,+1)\}),S)</math>
<math> = \frac{1}{\tau_2(e)}\int_{e' \in \partial(\{(e,+1)\})} S(e')</math>
<math> = \nabla \times S|_e</math>

The counterclockwise boundary of multi-surface <math>S</math> is therefore the curl: <math>\nabla \times S</math>.

Stoke's theorem becomes:

<math>\int_{e \in \nabla \times S} f(e) = \int_{e \in S} (\nabla \times f)|_e</math> which is equivalent to:
<math>\Sigma_{\mathcal{E}}(\nabla \times S, f) = \Sigma_{\mathcal{E}}(\nabla \times f, S)</math> which is an exact rephrasing of the condition <math>\Sigma_{\mathcal{E}}(\partial S, f) = \Sigma_{\mathcal{E}}(\partial f, S)</math>.

=== The Curl in Various Coordinate Systems ===

All of Cartesian, cylindrical, and spherical coordinate systems will be modeled using the [[Calculus/Discrete vector calculus#Model of Curvilinear coordinates|discrete model of curvilinear coordinates]]. The same shorthand notation will be used. However, the graph will be modified as follows:

Similar to the notation from the [[Calculus/Discrete vector calculus#Model of Curvilinear coordinates|discrete model of curvilinear coordinates]], the following notation will be used:

* Given any <math>c \in \{1,2,3\}</math>, <math>c + 1 = \left\{\begin{array}{cc} c+1 & (c \in \{1,2\}) \\ c-2 & (c = 3) \end{array}\right.</math> and <math>c + 2 = \left\{\begin{array}{cc} c+2 & (c = 1) \\ c-1 & (c \in \{2,3\}) \end{array}\right.</math>
* Given any number <math>k \in \R</math>, then <math>\Z + k</math> denotes the set of integers where <math>k</math> has been added to each integer.
* Given any number <math>k \in \R</math>, and subset <math>C \subseteq \{1,2,3\}</math>, then <math>[k]_C</math> denotes the triple <math>[k]_C = \left(\left\{\begin{array}{cc} k & (1 \in C) \\ 0 & (1 \notin C) \end{array}\right.,\left\{\begin{array}{cc} k & (2 \in C) \\ 0 & (2 \notin C) \end{array}\right.,\left\{\begin{array}{cc} k & (3 \in C) \\ 0 & (3 \notin C) \end{array}\right.\right)</math>
* Given any two triples of numbers <math>(x_1,x_2,x_3), (y_1,y_2,y_3) \in \R^3</math>, then <math>(x_1,x_2,x_3) + (y_1,y_2,y_3) = (x_1+y_1,x_2+y_2,x_3+y_3)</math>
* Given any two triples of numbers <math>(x_1,x_2,x_3), (y_1,y_2,y_3) \in \R^3</math>, then <math>(x_1,x_2,x_3)(y_1,y_2,y_3) = (x_1y_1,x_2y_2,x_3y_3)</math>

After choosing a resolution for each coordinate: <math>\Delta\Xi = (\Delta\xi_1, \Delta\xi_2, \Delta\xi_3)</math>, the graph <math>\mathcal{G} = \langle \mathcal{N}, \mathcal{E} \rangle</math> has the following nodes and edges:

* For each <math>I \in \Z^3 \cup (\Z + 0.5)^3</math>, there exists a node <math>n_I \in \mathcal{N}</math> at position <math>\omega_0(n_I) = \omega_0(I\Delta\Xi)</math>. The volume of <math>n_I</math> is <math>\tau_3(n_I) = \prod_{c=1}^3 (\Delta\xi_c|\mathbf{u}_{\xi_c}(I\Delta\Xi)|)</math>.
* For each <math>c = 1,2,3</math> and <math>I \in \Z^3 \cup (\Z + 0.5)^3</math>, there exists an edge <math>e_{\xi_c,I} \in \mathcal{E}</math> such that <math>t_\bullet(e_{\xi_c,I}) = n_I</math> and <math>t^\bullet(e_{\xi_c,I}) = n_{I+[1]_c}</math>. The length of <math>e_{\xi_c,I}</math> is <math>\tau_1(e_{\xi_c,I}) = \Delta\xi_c|\mathbf{u}_{\xi_c}((I+[0.5]_c)\Delta\Xi)|</math>. The thickness of <math>e_{\xi_c,I}</math> is <math>\tau_2(e_{\xi_c,I}) = \prod_{c' \neq c} (\Delta\xi_{c'}|\mathbf{u}_{\xi_{c'}}((I+[0.5]_c)\Delta\Xi)|)</math>.

The graph consists of two interleaved 3D-grids, with nodes <math>n_I</math> where <math>I \in \Z^3</math> forming one grid and nodes <math>n_I</math> where <math>I \in (\Z+0.5)^3</math> forming the other grid.

Lastly, for any edge <math>e_{\xi_c,I} \in \mathcal{E}</math>, the counter-clockwise boundary of the surface <math>\{(e_{\xi_c,I},+1)\}</math> is the path: <math>\langle (e_{\xi_{c+1},I + [0.5]_c - [0.5]_{c+1} - [0.5]_{c+2}},+1) , (e_{\xi_{c+2},I + [0.5]_c + [0.5]_{c+1} - [0.5]_{c+2}},+1) , (e_{\xi_{c+1},I + [0.5]_c - [0.5]_{c+1} + [0.5]_{c+2}},-1) , (e_{\xi_{c+2},I + [0.5]_c - [0.5]_{c+1} - [0.5]_{c+2}},-1) \rangle</math>. This path is a square that passes through nodes <math>n_{I + [0.5]_c - [0.5]_{c+1} - [0.5]_{c+2}}; n_{I + [0.5]_c + [0.5]_{c+1} - [0.5]_{c+2}}; n_{I + [0.5]_c + [0.5]_{c+1} + [0.5]_{c+2}}; n_{I + [0.5]_c - [0.5]_{c+1} + [0.5]_{c+2}}</math>

Given an arbitrary vector field <math>\mathbf{F}(\Xi) = \sum_{c=1}^3 F_{\xi_c}(\Xi)\mathbf{v}_{\xi_c}(\Xi)</math>, the edge based function that approximates <math>\mathbf{F}</math> is <math>F_\text{approx}(e_{\xi_c,I}) = F_{\xi_c}((I+[0.5]_c)\Delta\Xi)</math> for all <math>c \in \{1,2,3\}</math> and <math>I \in \Z^3 \cup (\Z + 0.5)^3</math>.

For each <math>e_{\xi_c,I} \in \mathcal{E}</math>, the curl (circulation density/counter-clockwise boundary) of <math>F_\text{approx}</math> is:

<math>\nabla \times F_\text{approx}|_{e_{\xi_c,I}} = \frac{1}{\tau_2(e_{\xi_c,I})}\int_{e \in \partial(\{(e_{\xi_c,I},+1)\})} F_\text{approx}(e)</math>
 
<math> = \frac{1}{\tau_2(e_{\xi_c,I})}\bigg(</math>
<math> F_\text{approx}(e_{\xi_{c+1},I + [0.5]_c - [0.5]_{c+1} - [0.5]_{c+2}})\tau_1(e_{\xi_{c+1},I + [0.5]_c - [0.5]_{c+1} - [0.5]_{c+2}}) </math>
<math> + F_\text{approx}(e_{\xi_{c+2},I + [0.5]_c + [0.5]_{c+1} - [0.5]_{c+2}})\tau_1(e_{\xi_{c+2},I + [0.5]_c + [0.5]_{c+1} - [0.5]_{c+2}}) </math>
<math> - F_\text{approx}(e_{\xi_{c+1},I + [0.5]_c - [0.5]_{c+1} + [0.5]_{c+2}})\tau_1(e_{\xi_{c+1},I + [0.5]_c - [0.5]_{c+1} + [0.5]_{c+2}}) </math>
<math> - F_\text{approx}(e_{\xi_{c+2},I + [0.5]_c - [0.5]_{c+1} - [0.5]_{c+2}})\tau_1(e_{\xi_{c+2},I + [0.5]_c - [0.5]_{c+1} - [0.5]_{c+2}})\bigg) </math>

<math> = \frac{1}{|\mathbf{u}_{\xi_{c+1}}((I+[0.5]_c)\Delta\Xi)||\mathbf{u}_{\xi_{c+2}}((I+[0.5]_c)\Delta\Xi)|\Delta\xi_{c+1}\Delta\xi_{c+2}}\bigg(</math>
<math> F_{\xi_{c+1}}((I + [0.5]_c)\Delta\Xi - [0.5\Delta\xi_{c+2}]_{c+2})|\mathbf{u}_{\xi_{c+1}}((I + [0.5]_c)\Delta\Xi - [0.5\Delta\xi_{c+2}]_{c+2})|\Delta\xi_{c+1} </math>
<math> + F_{\xi_{c+2}}((I + [0.5]_c)\Delta\Xi + [0.5\Delta\xi_{c+1}]_{c+1})|\mathbf{u}_{\xi_{c+2}}((I + [0.5]_c)\Delta\Xi + [0.5\Delta\xi_{c+1}]_{c+1})|\Delta\xi_{c+2} </math>
<math> - F_{\xi_{c+1}}((I + [0.5]_c)\Delta\Xi + [0.5\Delta\xi_{c+2}]_{c+2})|\mathbf{u}_{\xi_{c+1}}((I + [0.5]_c)\Delta\Xi + [0.5\Delta\xi_{c+2}]_{c+2})|\Delta\xi_{c+1} </math>
<math> - F_{\xi_{c+2}}((I + [0.5]_c)\Delta\Xi - [0.5\Delta\xi_{c+1}]_{c+1})|\mathbf{u}_{\xi_{c+2}}((I + [0.5]_c)\Delta\Xi - [0.5\Delta\xi_{c+1}]_{c+1})|\Delta\xi_{c+2} </math>
<math>\bigg)</math>

Using <math>I\Delta\Xi</math> and <math>(I + [0.5]_c)\Delta\Xi</math> as approximations for <math>\Xi</math>, gives:

<math>\nabla \times \mathbf{F}|_\Xi \approx \sum_{c=1}^3 \frac{1}{|\mathbf{u}_{\xi_{c+1}}(\Xi)||\mathbf{u}_{\xi_{c+2}}(\Xi)|\Delta\xi_{c+1}\Delta\xi_{c+2}}\bigg(</math>
<math> F_{\xi_{c+1}}(\Xi - [0.5\Delta\xi_{c+2}]_{c+2})|\mathbf{u}_{\xi_{c+1}}(\Xi - [0.5\Delta\xi_{c+2}]_{c+2})|\Delta\xi_{c+1} </math>
<math> + F_{\xi_{c+2}}(\Xi + [0.5\Delta\xi_{c+1}]_{c+1})|\mathbf{u}_{\xi_{c+2}}(\Xi + [0.5\Delta\xi_{c+1}]_{c+1})|\Delta\xi_{c+2} </math>
<math> - F_{\xi_{c+1}}(\Xi + [0.5\Delta\xi_{c+2}]_{c+2})|\mathbf{u}_{\xi_{c+1}}(\Xi + [0.5\Delta\xi_{c+2}]_{c+2})|\Delta\xi_{c+1} </math>
<math> - F_{\xi_{c+2}}(\Xi - [0.5\Delta\xi_{c+1}]_{c+1})|\mathbf{u}_{\xi_{c+2}}(\Xi - [0.5\Delta\xi_{c+1}]_{c+1})|\Delta\xi_{c+2} </math>
<math>\bigg)\mathbf{v}_{\xi_c}(\Xi)</math>

<math> = \sum_{c=1}^3 \frac{1}{|\mathbf{u}_{\xi_{c+1}}(\Xi)||\mathbf{u}_{\xi_{c+2}}(\Xi)|}\bigg(</math>
<math> \frac{1}{\Delta\xi_{c+1}}(F_{\xi_{c+2}}(\Xi + [0.5\Delta\xi_{c+1}]_{c+1})|\mathbf{u}_{\xi_{c+2}}(\Xi + [0.5\Delta\xi_{c+1}]_{c+1})| - F_{\xi_{c+2}}(\Xi - [0.5\Delta\xi_{c+1}]_{c+1})|\mathbf{u}_{\xi_{c+2}}(\Xi - [0.5\Delta\xi_{c+1}]_{c+1})|)</math>
<math> - \frac{1}{\Delta\xi_{c+2}}(F_{\xi_{c+1}}(\Xi + [0.5\Delta\xi_{c+2}]_{c+2})|\mathbf{u}_{\xi_{c+1}}(\Xi + [0.5\Delta\xi_{c+2}]_{c+2})| - F_{\xi_{c+1}}(\Xi - [0.5\Delta\xi_{c+2}]_{c+2})|\mathbf{u}_{\xi_{c+1}}(\Xi - [0.5\Delta\xi_{c+2}]_{c+2})|)</math>
<math>\bigg)\mathbf{v}_{\xi_c}(\Xi)</math>

As <math>\Delta\xi_1, \Delta\xi_2, \Delta\xi_3 \to 0^+</math>, it is the case that: 
<math>\nabla \times \mathbf{F}|_{\Xi} = \sum_{c=1}^3 \frac{1}{|\mathbf{u}_{\xi_{c+1}}(\Xi)||\mathbf{u}_{\xi_{c+2}}(\Xi)|}\left(\frac{\partial}{\partial\xi_{c+1}}(F_{\xi_{c+2}}(\Xi)|\mathbf{u}_{\xi_{c+2}}(\Xi)|)\bigg|_{\Xi} - \frac{\partial}{\partial\xi_{c+2}}(F_{\xi_{c+1}}(\Xi)|\mathbf{u}_{\xi_{c+1}}(\Xi)|)\bigg|_{\Xi}\right)\mathbf{v}_{\xi_c}(\Xi)</math>

==== The Curl in Cartesian Coordinates ====

Applying the curvilinear coordinate model to Cartesian coordinates means that <math>\xi_1 = x</math>, <math>\xi_2 = y</math>, and <math>\xi_3 = z</math>. It is also the case that <math>\mathbf{u}_x(x,y,z) = \mathbf{i}</math>, <math>\mathbf{u}_y(x,y,z) = \mathbf{j}</math>, and <math>\mathbf{u}_z(x,y,z)=\mathbf{k}</math>. Applying the formula for the curl that was derived from the discrete model of curvilinear coordinates to vector field <math>\mathbf{F} = F_x\mathbf{i} + F_y\mathbf{j} + F_z\mathbf{k}</math> gives:

<math>\nabla \times \mathbf{F}|_{(x,y,z)} = \frac{1}{1 \cdot 1}\left(\frac{\partial}{\partial y}(1 \cdot F_z)\bigg|_{(x,y,z)} - \frac{\partial}{\partial z}(1 \cdot F_y)\bigg|_{(x,y,z)}\right)\mathbf{i} </math> <math> + \frac{1}{1 \cdot 1}\left(\frac{\partial}{\partial z}(1 \cdot F_x)\bigg|_{(x,y,z)} - \frac{\partial}{\partial x}(1 \cdot F_z)\bigg|_{(x,y,z)}\right)\mathbf{j} </math> <math> + \frac{1}{1 \cdot 1}\left(\frac{\partial}{\partial x}(1 \cdot F_y)\bigg|_{(x,y,z)} - \frac{\partial}{\partial y}(1 \cdot F_x)\bigg|_{(x,y,z)}\right)\mathbf{k}</math>

<math> = \left(\frac{\partial F_z}{\partial y} - \frac{\partial F_y}{\partial z}\right)\mathbf{i} + \left(\frac{\partial F_x}{\partial z} - \frac{\partial F_z}{\partial x}\right)\mathbf{j} + \left(\frac{\partial F_y}{\partial x} - \frac{\partial F_x}{\partial y}\right)\mathbf{k}</math>

==== The Curl in Cylindrical Coordinates ====

Applying the curvilinear coordinate model to cylindrical coordinates means that <math>\xi_1 = \rho</math>, <math>\xi_2 = \phi</math>, and <math>\xi_3 = z</math>. It is also the case that <math>\mathbf{u}_\rho(\rho,\phi,z) = \hat{\mathbf{\rho}}</math>, <math>\mathbf{u}_\phi(\rho,\phi,z) = \rho\hat{\mathbf{\phi}}</math>, and <math>\mathbf{u}_z(\rho,\phi,z) = \hat{\mathbf{z}}</math>. Applying the formula for the curl that was derived from the discrete model of curvilinear coordinates to vector field <math>\mathbf{F} = F_\rho\hat{\mathbf{\rho }} + F_\phi\hat{\mathbf{\phi}} + F_z\hat{\mathbf{z}}</math> gives:

<math>\nabla \times \mathbf{F}|_{(\rho,\phi,z)} = \frac{1}{\rho \cdot 1}\left(\frac{\partial}{\partial\phi}(1 \cdot F_z) - \frac{\partial}{\partial z}(\rho \cdot F_\phi)\right)\hat{\mathbf{\rho}} </math><math> + \frac{1}{1 \cdot 1}\left(\frac{\partial}{\partial z}(1 \cdot F_\rho) - \frac{\partial}{\partial\rho}(1 \cdot F_z)\right)\hat{\mathbf{\phi}} </math><math> + \frac{1}{1 \cdot \rho}\left(\frac{\partial}{\partial\rho}(\rho \cdot F_\phi) - \frac{\partial}{\partial\phi}(1 \cdot F_\rho)\right)\hat{\mathbf{z}}</math>

<math> = \left(\frac{1}{\rho}\frac{\partial F_z}{\partial\phi} - \frac{\partial F_\phi}{\partial z}\right)\hat{\mathbf{\rho}} + \left(\frac{\partial F_\rho}{\partial z} - \frac{\partial F_z}{\partial\rho}\right)\hat{\mathbf{\phi}} + \frac{1}{\rho}\left(\frac{\partial}{\partial\rho}(\rho \cdot F_\phi) - \frac{\partial F_\rho}{\partial\phi}\right)\hat{\mathbf{z}}</math>

==== The Curl in Spherical Coordinates ====

Applying the curvilinear coordinate model to spherical coordinates means that <math>\xi_1 = r</math>, <math>\xi_2 = \theta</math>, and <math>\xi_3 = \phi</math>. It is also the case that <math>\mathbf{u}_r(r,\theta,\phi) = \hat{\mathbf{r}}</math>, <math>\mathbf{u}_\theta(r,\theta,\phi) = r\hat{\mathbf{\theta}}</math>, and <math>\mathbf{u}_\phi(r,\theta,\phi) = r\sin\theta\hat{\mathbf{\phi}}</math>. Applying the formula for the curl that was derived from the discrete model of curvilinear coordinates to vector field <math>\mathbf{F} = F_r\hat{\mathbf{r}} + F_\theta\hat{\mathbf{\theta}} + F_\phi\hat{\mathbf{\phi}}</math> gives:

<math>\nabla \times \mathbf{F}|_{(r,\theta,\phi)} = </math><math> \frac{1}{r \cdot r\sin\theta}\left(\frac{\partial}{\partial\theta}(r\sin\theta \cdot F_\phi) - \frac{\partial}{\partial\phi}(r \cdot F_\theta)\right)\hat{\mathbf{r}} </math><math> + \frac{1}{r\sin\theta \cdot 1}\left(\frac{\partial}{\partial\phi}(1 \cdot F_r) - \frac{\partial}{\partial r}(r\sin\theta \cdot F_\phi)\right)\hat{\mathbf{\theta}} </math><math> + \frac{1}{1 \cdot r}\left(\frac{\partial}{\partial r}(r \cdot F_\theta) - \frac{\partial}{\partial\theta}(1 \cdot F_r)\right)\hat{\mathbf{\phi}}</math>

<math> = \frac{1}{r\sin\theta}\left(\frac{\partial}{\partial\theta}(\sin\theta \cdot F_\phi) - \frac{\partial F_\theta}{\partial\phi}\right)\hat{\mathbf{r}} + \frac{1}{r}\left(\frac{1}{\sin\theta}\frac{\partial F_r}{\partial\phi} - \frac{\partial}{\partial r}(r \cdot F_\phi)\right)\hat{\mathbf{\theta}} + \frac{1}{r}\left(\frac{\partial}{\partial r}(r \cdot F_\theta) - \frac{\partial F_r}{\partial\theta}\right)\hat{\mathbf{\phi}}</math>

== Intersections ==

=== Identities ===

== Converting Paths to Surfaces and Surfaces to Paths ==

{{BookCat}}