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 $x$ as an example, takes on discrete quantities separated by small intervals of width $\Delta x$. $\Delta x$ is a measure of the model's "coarseness" (when $\Delta x$ is increasing) or "fineness" (when $\Delta x$ is decreasing). When handling quantities at the scale of $\Delta x$, an approximation for which the error relative to $\Delta x$ vanishes as $\Delta x$ becomes infinitely small, becomes exact when $x$ is treated as a continuous variable.

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

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

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

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 $\Delta x$ approach 0 as $\Delta x\to 0$ as opposed to the absolute error, the following example will be used:

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

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

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

$|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|$

$\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}}}$

$\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)$ $=\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)$ $=(b-a)\max _{j}\left({\frac {|f(x_{j}^{*},\Delta x_{j})-g(x_{j}^{*},\Delta x_{j})|}{\Delta x_{j}}}\right)$

For each $j=1,2,\dots ,N$, it is assumed that ${\frac {|f(x_{j}^{*},\Delta x_{j})-g(x_{j}^{*},\Delta x_{j})|}{\Delta x_{j}}}\to 0$ as $\Delta x_{j}\to 0^{+}$. Therefore $\max _{j}\left({\frac {|f(x_{j}^{*},\Delta x_{j})-g(x_{j}^{*},\Delta x_{j})|}{\Delta x_{j}}}\right)\to 0$ and hence $|F-G|\to 0$ as $N\to +\infty$ and $\max(\Delta x_{i})\to 0^{+}$. Therefore $F_{\infty }=G_{\infty }$.

The left image shows the sums $F$ and $G$. The white area of each strip depicts $f$, and the red area of each strip depicts $g-f$. The total area highlighted in red depicts the difference between $G$ and $F$. As $N$ 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" ${\mathcal {G}}=\langle {\mathcal {N}},{\mathcal {E}}\rangle$ consists of a set ${\mathcal {N}}$ of nodes/vertices and a set ${\mathcal {E}}$ of directed edges. For each edge $e\in {\mathcal {E}}$, $e$ has an origin node $t_{\bullet }(e)\in {\mathcal {N}}$ that is the beginning of edge $e$, and a terminal node $t^{\bullet }(e)\in {\mathcal {N}}$ that is the end of edge $e$. It will also be allowed that edges can form loops ($t_{\bullet }(e)=t^{\bullet }(e)$) and that multiple edges can have both the same beginning and end ($t_{\bullet }(e_{1})=t_{\bullet }(e_{2})$ and $t^{\bullet }(e_{1})=t^{\bullet }(e_{2})$), though this will be irrelevant for the cases that will be considered.

An arbitrary function $f:{\mathcal {N}}\to \mathbb {R}$ over the set of nodes will be referred to as a "node based function", and is effectively a scalar field. The notation $\langle n\in {\mathcal {N}}:f(n)\rangle$ 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: $\{n_{1}\mapsto f(n_{1}),n_{2}\mapsto f(n_{2}),\dots \}$

An arbitrary function $f:{\mathcal {E}}\to \mathbb {R}$ over the set of edges will be referred to as an "edge based function", and is effectively a scalar or vector field. The notation $\langle e\in {\mathcal {E}}:f(e)\rangle$ 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: $\{e_{1}\mapsto f(e_{1}),e_{2}\mapsto f(e_{2}),\dots \}$. 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.

An example directed graph.

An example directed graph is shown to the right. The node set is ${\mathcal {N}}=\{n_{1},n_{2},n_{3},n_{4},n_{5}\}$ and the edge set is ${\mathcal {E}}=\{e_{1},e_{2},e_{3},e_{4},e_{5},e_{6}\}$. The origin and terminal nodes of each edge are: $(t_{\bullet }(e_{1}),t^{\bullet }(e_{1}))=(n_{1},n_{2})$; $(t_{\bullet }(e_{2}),t^{\bullet }(e_{2}))=(n_{1},n_{3})$; $(t_{\bullet }(e_{3}),t^{\bullet }(e_{3}))=(n_{2},n_{4})$; $(t_{\bullet }(e_{4}),t^{\bullet }(e_{4}))=(n_{4},n_{3})$; $(t_{\bullet }(e_{5}),t^{\bullet }(e_{5}))=(n_{4},n_{5})$; and $(t_{\bullet }(e_{6}),t^{\bullet }(e_{6}))=(n_{5},n_{1})$.

Interpreting nodes and edges

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 $n\in {\mathcal {N}}$ will correspond to both a point $\mathbf {q} (n)$ in space, and a volume $\Omega (n)$ that contains $\mathbf {q} (n)$. $V(n)$ will denote the volume of $\Omega (n)$, also referred to as the volume of node $n$. The volume of each node must be positive: $V(n)>0$ for all $n\in {\mathcal {N}}$. In the example directed graph, an example choice of volumes for each node is $V=\{n_{1}\mapsto 3,n_{2}\mapsto 2,n_{3}\mapsto 2,n_{4}\mapsto 1,n_{5}\mapsto 4\}$.

Each edge $e\in {\mathcal {E}}$ will correspond to both a curve $C(e)$ originating from $\mathbf {q} (t_{\bullet }(e))$ and terminating on $\mathbf {q} (t^{\bullet }(e))$, and a surface $\sigma (e)$ that separates $\Omega (t_{\bullet }(e))$ from $\Omega (t^{\bullet }(e))$ oriented from $\Omega (t_{\bullet }(e))$ to $\Omega (t^{\bullet }(e))$.

$l(e)$ will denote the distance $|\mathbf {q} (t^{\bullet }(e))-\mathbf {q} (t_{\bullet }(e))|$ from $\mathbf {q} (t_{\bullet }(e))$ to $\mathbf {q} (t^{\bullet }(e))$, which will be referred to as the length of edge $e$. The length of each edge must be positive: $l(e)>0$ for all $e\in {\mathcal {E}}$. $A(e)$ will denote the area of $\sigma (e)$ projected onto a plane that is perpendicular to the displacement $\mathbf {q} (t^{\bullet }(e))-\mathbf {q} (t_{\bullet }(e))$, which will be referred to as the cross-sectional area, thickness, or width of edge $e$. The cross-sectional area of each edge must be positive: $A(e)>0$ for all $e\in {\mathcal {E}}$. In the example directed graph, an example choice of lengths for each edge is $l=\{e_{1}\mapsto 1,e_{2}\mapsto 1,e_{3}\mapsto 2,e_{4}\mapsto 1,e_{5}\mapsto 3,e_{6}\mapsto 4\}$, and an example choice of areas for each edge is $A=\{e_{1}\mapsto 1,e_{2}\mapsto 2,e_{3}\mapsto 1,e_{4}\mapsto 2,e_{5}\mapsto 3,e_{6}\mapsto 3\}$.

In summary an arbitrary directed graph ${\mathcal {G}}$ will be characterized by:

• the set of nodes ${\mathcal {N}}$
• the set of edges ${\mathcal {E}}$
• the start $t_{\bullet }(e)$ and end $t^{\bullet }(e)$ of each edge $e\in {\mathcal {E}}$
• the volume $V(n)$ of each node $n\in {\mathcal {N}}$
• the length $l(e)$ and area $A(e)$ of each edge $e\in {\mathcal {E}}$

The point $\mathbf {q} (n)$ and space $\Omega (n)$ associated with each node $n\in {\mathcal {N}}$, and the curve $C(e)$ and surface $\sigma (e)$ associated with each edge $e\in {\mathcal {E}}$, 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 $n_{0}\in {\mathcal {N}}$, which corresponds to the space point $\mathbf {q} (n_{0})$.

Point $n_{0}$ can be denoted by the node based function $\delta _{0}(n;n_{0}):{\mathcal {N}}\to \mathbb {R}$ where for each node $n\in {\mathcal {N}}$ it is the case that $\delta _{0}(n;n_{0})=\left\{{\begin{array}{cc}1/V(n)&(n=n_{0})\\0&(n\neq n_{0})\\\end{array}}\right.$. The term ${\frac {1}{V(n)}}$ arises from the fact that the point is "spread out" over the volume of $n$.

A superposition of points is a "multi-point". Any node based function that returns real values can denote a multi-point.

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

Volumes and multi-volumes

The example directed graph with the example multi-volume.

The discrete analog to volumes using a directed graph is that a volume $U$ is a collection of nodes: $U\subseteq {\mathcal {N}}$. The combined volume region denoted by $U$ is $\Omega (U)=\bigcup _{n\in U}\Omega (n)$ and the total volume amount is $V(U)=\sum _{n\in U}V(n)$.

Volume $U$ can be denoted by a node based membership function: $\delta _{3}(n;U):{\mathcal {N}}\to \mathbb {R}$ where for each node $n\in {\mathcal {N}}$, it is the case that $\delta _{3}(n;U)=\left\{{\begin{array}{cc}1&(n\in U)\\0&(n\notin U)\\\end{array}}\right.$

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 $f=\{n_{1}\mapsto -2,n_{2}\mapsto 0,n_{3}\mapsto 4,n_{4}\mapsto 0.5,n_{5}\mapsto 0.5\}$ can be expressed as the following superpositions:

$f(n)=-2\delta _{3}(n;\{n_{1}\})+4\delta _{3}(n;\{n_{3}\})+0.5\delta _{3}(n;\{n_{4},n_{5}\})$ is the superposition of volumes $\{n_{1}\},\{n_{3}\},\{n_{4},n_{5}\}$ with the respective weights $-2,4,0.5$.

$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}\})$ is the superposition of volumes $\{n_{1}\},\{n_{3},n_{4},n_{5}\},\{n_{3}\}$ with the respective weights $-2,0.5,3.5$.

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

The example directed graph with the example multi-path.

The discrete analog to paths using a directed graph is that a path $P$ is a series of edges and orientations: $P=\langle (e_{1},s_{1}),(e_{2},s_{2}),\dots ,(e_{k},s_{k})\rangle$ where $e_{i}\in {\mathcal {E}}$ is the $i^{\text{th}}$ edge and $s_{i}\in \{+1,-1\}$ is an indicator that is $+1$ if $e_{i}$ is traversed in the preferred direction, and is $-1$ if $e_{i}$ is traversed in the opposite direction. The path must be unbroken: for each $i\in \{1,2,\dots ,k-1\}$ it must be that $\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.$. The combined path denoted by $P$ is $C(P)=\bigcup _{i=1}^{k}s_{i}C(e_{i})$ (path $-C(e_{i})$ is path $C(e_{i})$ the orientation reversed) and the total length is $l(P)=\sum _{i=1}^{k}l(e_{i})$.

A path $P$ can be denoted by the edge based function $\delta _{1}(e;P):{\mathcal {E}}\to \mathbb {R}$ where for each edge $e\in {\mathcal {E}}$, it is the case that $\delta _{1}(e;P)=\sum _{i=1}^{k}s_{i}\delta _{1}(e;e_{i})$ where $\delta _{1}(e;e_{i})=\left\{{\begin{array}{cc}1/A(e)&(e=e_{i})\\0&(e\neq e_{i})\\\end{array}}\right.$. The term ${\frac {1}{A(e)}}$ arises from the fact that the path is "spread out" over the thickness of $e$.

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 $A(e)=1$ for all $e\in {\mathcal {E}}$, the edge based function $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\}$ can be expressed as the following superpositions:

$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 )$ is the superposition of paths $\langle (e_{1},+1)\rangle ,\langle (e_{2},+1)\rangle ,\langle (e_{3},+1)\rangle$, and $\langle (e_{4},+1)\rangle$ with the respective weights $1,-1,1,2$.

$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 )$ is the superposition of paths $\langle (e_{2},-1),(e_{1},+1),(e_{3},+1),(e_{4},+1)\rangle$, and $\langle (e_{4},+1)\rangle$ with the respective weights $1,1$.

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.

Surfaces and multi-surfaces

The example directed graph with the example multi-surface.

The discrete analog to oriented surfaces using a directed graph is that an oriented surface $S$ is a collection of edges and orientations: $S=\{(e_{1},s_{1}),(e_{2},s_{2}),\dots ,(e_{k},s_{k})\}$ where all of the edges $e_{1},e_{2},\dots ,e_{k}\in {\mathcal {E}}$ are distinct and $s_{i}\in \{+1,-1\}$ is an indicator that describes the surface's orientation. Each edge $e_{i}$ is an edge that passes through the surface, and $s_{i}$ is $+1$ if $e_{i}$ passes through the surface in the preferred direction, and $s_{i}$ is $-1$ if $e_{i}$ passes through the surface in the opposite direction. The combined surface denoted by $S$ is $\sigma (S)=\bigcup _{i=1}^{k}s_{i}\sigma (e_{i})$ (surface $-\sigma (e_{i})$ is surface $\sigma (e_{i})$ with the orientation reversed) and the total area is $A(S)=\sum _{i=1}^{k}A(e_{i})$.

A surface $S$ can be denoted by the edge based function $\delta _{2}(e;S):{\mathcal {E}}\to \mathbb {R}$ where for each edge $e\in {\mathcal {E}}$, it is the case that $\delta _{2}(e;S)=\left\{{\begin{array}{cc}+1/l(e)&((e,+1)\in S)\\-1/l(e)&((e,-1)\in S)\\0&({\text{else}})\\\end{array}}\right.$. The term ${\frac {1}{l(e)}}$ arises from the fact that the surface is "spread out" over the length of $e$.

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 $l(e)=1$ for all $e\in {\mathcal {E}}$, the edge based function $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\}$ can be expressed as the following superpositions:

$f(e)=2\delta _{2}(e;\{(e_{1},+1)\})+2\delta _{2}(e;\{(e_{2},+1)\})-\delta _{2}(e;\{(e_{6},+1)\})$ is the superposition of surfaces $\{(e_{1},+1)\},\{(e_{2},+1)\},$ and $\{(e_{6},+1)\}$ with the respective weights $2,2,-1$.

$f(e)=\delta _{2}(e;\{(e_{1},+1),(e_{2},+1),(e_{6},-1)\})+\delta _{2}(e;\{(e_{1},+1),(e_{2},+1)\})$ is the superposition of surfaces $\{(e_{1},+1),(e_{2},+1),(e_{6},-1)\}$ and $\{(e_{1},+1),(e_{2},+1)\}$ with the respective weights $1$ and $1$.

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.

Intersections

Multi-point multi-volume intersections

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 amount 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 $n_{0}\in {\mathcal {N}}$ and volume $U\subseteq {\mathcal {N}}$, let $\Sigma _{\mathcal {N}}(n_{0},U)=\left\{{\begin{array}{cc}1&(n_{0}\in U)\\0&(n_{0}\notin U)\end{array}}\right.$ denote an indicator function that returns 1 if point $n_{0}$ is contained by volume $U$ and returns 0 if otherwise.

Node based function $\delta _{0}(n;n_{0})=\left\{{\begin{array}{cc}1/V(n)&(n=n_{0})\\0&(n\neq n_{0})\end{array}}\right.$ denotes the point $n_{0}$, and node based function $\delta _{3}(n;U)=\left\{{\begin{array}{cc}1&(n\in U)\\0&(n\notin U)\end{array}}\right.$ denotes the volume $U$.

The indicator function $\Sigma _{\mathcal {N}}(n_{0},U)$ can be computed via the following sum: $\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)V(n)$. The term $V(n)$ cancels out the ${\frac {1}{V(n)}}$ in $\delta _{0}(n;n_{0})$.

$\Sigma _{\mathcal {N}}(n_{0},U)$ counts the number of times point $n_{0}$ occurs in volume $U$ which can only be 0 or 1 for simple points and volumes. This count can be generalized to sets of points. Let $W\subseteq {\mathcal {N}}$ denote a collection of points. The collection of points can be denoted by the node based function $\delta _{0}(n;W)=\left\{{\begin{array}{cc}1/V(n)&(n\in W)\\0&(n\notin W)\end{array}}\right.$. The number of points from $W$ that are contained by volume $U$ can again be computed by the sum: $\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)V(n)$.

In the general case, given a node based function $f$ that denotes a multi-point, and a node based function $g$ that denotes a multi-volume, then the total weight of all the occurrences of $f$ being inside of $g$ is $\Sigma _{\mathcal {N}}(f,g)=\sum _{n\in {\mathcal {N}}}f(n)g(n)V(n)$.

It should also be noted that $\Sigma _{\mathcal {N}}(f,g)$ is commutative: $\Sigma _{\mathcal {N}}(f,g)=\Sigma _{\mathcal {N}}(g,f)$

Multi-path multi-surface intersections

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 $P=\langle (e_{P,1},s_{P,1}),(e_{P,2},s_{P,2}),...,(e_{P,k_{P}},s_{P,k_{P}})\rangle$ and surface $S=\{(e_{S,1},s_{S,1}),(e_{S,2},s_{S,2}),...,(e_{S,k_{S}},s_{S,k_{S}})\}$, let $\Sigma _{\mathcal {E}}(P,S)$ denote the number of times path $P$ passes through surface $S$ in the preferred direction minus the number of times path $P$ passes through surface $S$ in the opposite direction. For an arbitrary $i=1,2,...,k_{P}$, and an arbitrary $j=1,2,...,k_{S}$, if $e_{P,i}=e_{S,j}$ then the $i^{\text{th}}$ link of path $P$ coincides with the $j^{\text{th}}$ segment of surface $S$. If $s_{P,i}s_{S,j}=+1$, then $s_{P,i}=s_{S,j}$ and $P$ passes through $S$ in the preferred direction. If $s_{P,i}s_{S,j}=-1$, then $s_{P,i}\neq s_{S,j}$ and $P$ passes through $S$ in the opposite direction.

$\Sigma _{\mathcal {E}}(P,S)$ can be computed via the following sum: $\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)A(e)l(e)$. The term $A(e)$ cancels out the ${\frac {1}{A(e)}}$ from $\delta _{1}(e;P)$, and the term $l(e)$ cancels out the ${\frac {1}{l(e)}}$ from $\delta _{2}(e;S)$.

In the general case, given an edge based function $f$ that denotes a multi-path, and an edge based function $g$ that denotes a multi-surface, then the total weight of all the intersections of $f$ with $g$ is $\Sigma _{\mathcal {E}}(f,g)=\sum _{e\in {\mathcal {E}}}f(e)g(e)A(e)l(e)$.

It should also be noted that $\Sigma _{\mathcal {E}}(f,g)$ is commutative: $\Sigma _{\mathcal {E}}(f,g)=\Sigma _{\mathcal {E}}(g,f)$

Lumped coordinate system models

Model of Cartesian coordinates

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 $\mathbb {R} ^{3}$. To establish the model, the following is needed:

• A spacing $\Delta x$ for the x-coordinates.
• A spacing $\Delta y$ for the y-coordinates.
• A spacing $\Delta z$ for the z-coordinates.

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

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

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

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

The Cartesian coordinate associated with node $n_{(i,j,k)}$ is $\mathbf {q} (n_{(i,j,k)})=(i\Delta x,j\Delta y,k\Delta z)$, and the space associated with node $n_{(i,j,k)}$ is $\Omega (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]$.

The volume of node $n_{(i,j,k)}$ is $V(n_{(i,j,k)})=\Delta x\Delta y\Delta z$.

The length of edge $e_{x,(i,j,k)}$ is $l(e_{x,(i,j,k)})=\Delta x$. The length of edge $e_{y,(i,j,k)}$ is $l(e_{y,(i,j,k)})=\Delta y$. The length of edge $e_{z,(i,j,k)}$ is $l(e_{z,(i,j,k)})=\Delta z$.

The area of edge $e_{x,(i,j,k)}$ is $A(e_{x,(i,j,k)})=\Delta y\Delta z$. The area of edge $e_{y,(i,j,k)}$ is $A(e_{y,(i,j,k)})=\Delta z\Delta x$. The area of edge $e_{z,(i,j,k)}$ is $A(e_{z,(i,j,k)})=\Delta x\Delta y$.

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

Given a vector field $\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}$ over Cartesian coordinates, $\mathbf {F}$ will be approximated by the edge based function $F_{\text{approx}}:{\mathcal {E}}\to \mathbb {R}$. Function $F_{\text{approx}}$ is defined by $F_{\text{approx}}(e_{x,(i,j,k)})=F_{x}(i\Delta x,j\Delta y,k\Delta z)$; $F_{\text{approx}}(e_{y,(i,j,k)})=F_{y}(i\Delta x,j\Delta y,k\Delta z)$; and $F_{\text{approx}}(e_{z,(i,j,k)})=F_{z}(i\Delta x,j\Delta y,k\Delta z)$.

Model of Cylindrical coordinates

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 $\mathbb {R} ^{3}$. To establish the model, the following is needed:

• A spacing $\Delta \rho$ for the coordinate $\rho$.
• A large positive integer $N_{\phi }$ that will yield a spacing $\Delta \phi ={\frac {2\pi }{N_{\phi }}}$ for the coordinate $\phi$.
• A spacing $\Delta z$ for the coordinate $z$.

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

Clearly, due to the nature of cylindrical coordinates (such as the fact that coordinate $\phi$ is cyclical), not every triplet $(i,j,k)$ will be allowed. The following three domains will be defined: $I_{0}$, $I_{1}$, and $I_{2}$.

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

$(i,j,k)\in I_{1}\iff (i\geq 0\;{\text{and}}\;(i=0\implies j=0)\;{\text{and}}\;0\leq j<N_{\phi })$

$(i,j,k)\in I_{2}\iff (i\geq 0\;{\text{and}}\;0\leq j<N_{\phi })$

There exists a node $n_{(i,j,k)}$ for each $(i,j,k)\in I_{1}$.

There exists an edge $e_{\rho ,(i,j,k)}$ for each $(i,j,k)\in I_{2}$ such that $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.$ and $t^{\bullet }(e_{\rho ,(i,j,k)})=n_{(i+1,j,k)}$

There exists an edge $e_{\phi ,(i,j,k)}$ for each $(i,j,k)\in I_{0}$ such that $t_{\bullet }(e_{\phi ,(i,j,k)})=n_{(i,j,k)}$ and $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.$

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

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

For each $(i,j,k)\in I_{1}$ the cylindrical coordinate that corresponds to $n_{(i,j,k)}$ is $\mathbf {q} (n_{(i,j,k)})=(i\Delta \rho ,j\Delta \phi ,k\Delta z)$, and the space (described using cylindrical coordinates) associated with $n_{(i,j,k)}$ is $\Omega (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.$

For each $(i,j,k)\in I_{1}$ the volume of node $n_{(i,j,k)}$ is $V(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.$

For each $(i,j,k)\in I_{2}$ the length and area of edge $e_{\rho ,(i,j,k)}$ is $l(e_{\rho ,(i,j,k)})=\Delta \rho$ and $A(e_{\rho ,(i,j,k)})\approx l_{\phi }(i+{\frac {1}{2}})\cdot \Delta z$ respectively.

For each $(i,j,k)\in I_{0}$ the length and area of edge $e_{\phi ,(i,j,k)}$ is $l(e_{\phi ,(i,j,k)})\approx l_{\phi }(i)$ and $A(e_{\phi ,(i,j,k)})=\Delta \rho \cdot \Delta z$ respectively.

For each $(i,j,k)\in I_{1}$ the length and area of edge $e_{z,(i,j,k)}$ is $l(e_{z,(i,j,k)})=\Delta z$ and $A(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.$ respectively.

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

Given a vector field $\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} }}$ over cylindrical coordinates, $\mathbf {F}$ will be approximated by the edge based function $F_{\text{approx}}:{\mathcal {E}}\to \mathbb {R}$. Function $F_{\text{approx}}$ is defined by:

$\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.$

$\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)$

$\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)$

Model of Spherical coordinates

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 $\mathbb {R} ^{3}$. To establish the model, the following is needed:

• A spacing $\Delta r$ for the coordinate $r$.
• A large positive integer $N_{\theta }$ that will yield a spacing $\Delta \theta ={\frac {\pi }{N_{\theta }}}$ for the coordinate $\theta$.
• A large positive integer $N_{\phi }$ that will yield a spacing $\Delta \phi ={\frac {2\pi }{N_{\phi }}}$ for the coordinate $\phi$.

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

Clearly, due to the nature of spherical coordinates (such as the fact that coordinate $\phi$ is cyclical), not every triplet $(i,j,k)$ will be allowed. The following three domains will be defined: $I_{n}$, $I_{r}$, $I_{\theta }$, and $I_{\phi }$.

$(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))$

$(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))$

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

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

There exists a node $n_{(i,j,k)}$ for each $(i,j,k)\in I_{n}$.

There exists an edge $e_{r,(i,j,k)}$ for each $(i,j,k)\in I_{r}$ such that $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.$ and $t^{\bullet }(e_{r,(i,j,k)})=n_{(i+1,j,k)}$

There exists an edge $e_{\theta ,(i,j,k)}$ for each $(i,j,k)\in I_{\theta }$ such that $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.$ and $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.$

There exists an edge $e_{\phi ,(i,j,k)}$ for each $(i,j,k)\in I_{\phi }$ such that $t_{\bullet }(e_{\phi ,(i,j,k)})=n_{(i,j,k)}$ and $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.$

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

For each $(i,j,k)\in I_{n}$ the spherical coordinate that corresponds to $n_{(i,j,k)}$ is $\mathbf {q} (n_{(i,j,k)})=(i\Delta r,j\Delta \theta ,k\Delta \phi )$, and the space (described using spherical coordinates) associated with $n_{(i,j,k)}$ is $\Omega (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.$

For each $(i,j,k)\in I_{n}$ the volume of node $n_{(i,j,k)}$ is $V(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.$

For each $(i,j,k)\in I_{r}$ the length and area of edge $e_{r,(i,j,k)}$ is $l(e_{r,(i,j,k)})=\Delta r$ and $A(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.$ respectively.

For each $(i,j,k)\in I_{\theta }$ the length and area of edge $e_{\theta ,(i,j,k)}$ is $l(e_{\theta ,(i,j,k)})\approx l_{\theta }(i)$ and $A(e_{\theta ,(i,j,k)})\approx \Delta r\cdot l_{\phi }(i,j+{\frac {1}{2}})$ respectively.

For each $(i,j,k)\in I_{\phi }$ the length and area of edge $e_{\phi ,(i,j,k)}$ is $l(e_{\phi ,(i,j,k)})\approx l_{\phi }(i,j)$ and $A(e_{\phi ,(i,j,k)})=\Delta r\cdot l_{\theta }(i)$ respectively.

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

Given a vector field $\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 } }}$ over spherical coordinates, $\mathbf {F}$ will be approximated by the edge based function $F_{\text{approx}}:{\mathcal {E}}\to \mathbb {R}$. Function $F_{\text{approx}}$ is defined by:

$\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.$

$\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.$

$\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 )$

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 $\xi _{1},\xi _{2},\xi _{3}$.

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

• $I=(i,j,k)$ denotes an arbitrary triplet of integer indices from $\mathbb {Z} ^{3}$.

• For each $c=1,2,3$: $I+1_{c}=\left\{{\begin{array}{cc}(i+1,j,k)&(c=1)\\(i,j+1,k)&(c=2)\\(i,j,k+1)&(c=3)\end{array}}\right.$

• $\Xi =(\xi _{1},\xi _{2},\xi _{3})$ denotes an arbitrary curvilinear coordinate $(\xi _{1},\xi _{2},\xi _{3})$.

• For each $c=1,2,3$: $\Xi +\Delta \xi _{c}=\left\{{\begin{array}{cc}(\xi _{1}+\Delta \xi _{1},\xi _{2},\xi _{3})&(c=1)\\(\xi _{1},\xi _{2}+\Delta \xi _{2},\xi _{3})&(c=2)\\(\xi _{1},\xi _{2},\xi _{3}+\Delta \xi _{3})&(c=3)\end{array}}\right.$

• $\Delta \Xi =(\Delta \xi _{1},\Delta \xi _{2},\Delta \xi _{3})$ denotes an arbitrary triplet of differences in the coordinates $\xi _{1},\xi _{2},\xi _{3}$.

• $I\Delta \Xi =(i\Delta \xi _{1},j\Delta \xi _{2},k\Delta \xi _{3})$

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

For the purposes of simplicity, it will be assumed that the coordinates $\Xi =(\xi _{1},\xi _{2},\xi _{3})$ 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 of nodes ${\mathcal {N}}$ connected by directed edges ${\mathcal {E}}$ will be established:

Start by choosing "resolutions" for each coordinate $\xi _{1},\xi _{2},\xi _{3}$ denoted by $\Delta \xi _{1},\Delta \xi _{2},\Delta \xi _{3}$. The quantities $\Delta \xi _{1},\Delta \xi _{2},\Delta \xi _{3}$ should all be strictly positive, and ideally should be as small as possible.

For each triplet of integer indices $I\in \mathbb {Z} ^{3}$, a node $n_{I}\in {\mathcal {N}}$ will be created at position $\mathbf {q} (n_{I})=\mathbf {q} (I\Delta \Xi )$.

For each $c=1,2,3$ and for each triplet of integer indices $I\in \mathbb {Z} ^{3}$, there will exist edge $e_{\xi _{c},I}\in {\mathcal {E}}$ such that $t_{\bullet }(e_{\xi _{c},I})=n_{I}$ and $t^{\bullet }(e_{\xi _{c},I})=n_{I+1_{c}}$.

For each $I\in \mathbb {Z} ^{3}$, the volume of node $n_{I}$ is $V(n_{I})=\prod _{c=1}^{3}(\Delta \xi _{c}|\mathbf {u} _{\xi _{c}}(I\Delta \Xi )|)$ $=\left(\prod _{c=1}^{3}|\mathbf {u} _{\xi _{c}}(I\Delta \Xi )|\right)\left(\prod _{c=1}^{3}\Delta \xi _{c}\right)$

For each $c=1,2,3$, the length and area of edge $e_{\xi _{c},I}$ are $l(e_{\xi _{c},I})=\Delta \xi _{c}|\mathbf {u} _{\xi _{c}}(I\Delta \Xi )|$ and $A(e_{\xi _{c},I})=\prod _{c'\neq c}\Delta \xi _{c'}|\mathbf {u} _{\xi _{c'}}(I\Delta \Xi )|$ respectively.

Any node $n\in {\mathcal {N}}$ for which $V(n)=0$ should be deleted along with all connecting edges.

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

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

Given a vector field $\mathbf {F} (\Xi )=\sum _{c=1}^{3}F_{\xi _{c}}(\Xi )\mathbf {v} _{\xi _{c}}(\Xi )$ over the curvilinear coordinates, $\mathbf {F}$ will be approximated by the edge based function $F_{\text{approx}}:{\mathcal {E}}\to \mathbb {R}$. Function $F_{\text{approx}}$ is defined by $F_{\text{approx}}(e_{\xi _{c},I})=F_{\xi _{c}}(I\Delta \Xi )$ for each $c=1,2,3$ and for each $I\in \mathbb {Z} ^{3}$.

Discrete analog to Integrals

Volume Integrals

Let $U\subseteq {\mathcal {N}}$ denote an arbitrary volume. Given a node based function $f$, the "volume integral" of $f$ over $U$ is $\int _{n\in U}f(n)=\sum _{n\in U}f(n)V(n)$ (recall that $V(n)$ is the volume of node $n$, and $\Omega (n)$ is the volume of space represented by $n$).

If node based function $f:{\mathcal {N}}\to \mathbb {R}$ is an approximation of a scalar field $f_{\text{continuous}}:{\bigg (}\bigcup _{n\in {\mathcal {N}}}\Omega (n){\bigg )}\to \mathbb {R}$, where $\forall n\in {\mathcal {N}}:f(n)=f_{\text{continuous}}(\mathbf {q} (n))$, then $\int _{n\in U}f(n)$ is an approximation of the volume integral $\iiint _{\mathbf {q} \in \bigcup _{n\in U}\Omega (n)}f(\mathbf {q} )dV$.

Recall that $U$ can be described via the node based function $\delta _{3}(n;U)=\left\{{\begin{array}{cc}1&(n\in U)\\0&(n\notin U)\\\end{array}}\right.$. This allows the expression $\int _{n\in U}f(n)=\sum _{n\in {\mathcal {N}}}f(n)\delta _{3}(n;U)V(n)$. If node based function $g$ denotes a multi-volume, then $\int _{n\in g}f(n)=\sum _{n\in {\mathcal {N}}}f(n)g(n)V(n)=\Sigma _{\mathcal {N}}(f,g)$.

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

Path Integrals

Let $P=\langle (e_{1},s_{1}),(e_{2},s_{2}),...,(e_{k},s_{k})\rangle$ denote an arbitrary path. Given an edge based function $f$, the "path integral" of $f$ along $P$ is $\int _{e\in P}f(e)=\sum _{i=1}^{k}s_{i}f(e_{i})l(e_{i})$ (recall that $l(e)$ is the length of edge $e$, and $C(e)$ is the curve represented by $e$).

If edge based function $f:{\mathcal {E}}\to \mathbb {R}$ is an approximation of a vector field $\mathbf {F} _{\text{continuous}}:{\bigg (}\bigcup _{n\in {\mathcal {N}}}\Omega (n){\bigg )}\to \mathbb {R} ^{3}$, where $\forall e\in {\mathcal {E}}:f(e)=\mathbf {F} _{\text{continuous}}(\mathbf {q} (t_{\bullet }(e)))\cdot {\frac {\mathbf {q} (t^{\bullet }(e))-\mathbf {q} (t_{\bullet }(e))}{|\mathbf {q} (t^{\bullet }(e))-\mathbf {q} (t_{\bullet }(e))|}}$, then $\int _{e\in P}f(e)$ is an approximation of the path integral $\int _{\mathbf {q} \in \bigcup _{i=1}^{k}s_{i}C(e_{i})}\mathbf {F} (\mathbf {q} )\cdot d\mathbf {q}$.

Recall that $P$ can be described via the edge based function $\delta _{1}(e;P)=\sum _{i=1}^{k}s_{i}\left\{{\begin{array}{cc}1/A(e)&(e=e_{i})\\0&(e\neq e_{i})\\\end{array}}\right.$. This allows the expression $\int _{e\in P}f(e)=\sum _{e\in {\mathcal {E}}}f(e)\delta _{1}(e;P)A(e)l(e)$. If edge based function $g$ denotes a multi-path, then $\int _{e\in g}f(e)=\sum _{e\in {\mathcal {E}}}f(e)g(e)A(e)l(e)=\Sigma _{\mathcal {E}}(g,f)$.

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

Surface Integrals

Let $S=\{(e_{1},s_{1}),(e_{2},s_{2}),...,(e_{k},s_{k})\}$ denote an arbitrary surface. Given an edge based function $f$, the "surface integral" of $f$ over $S$ is $\int _{e\in S}f(e)=\sum _{(e,s)\in S}sf(e)A(e)$ (recall that $A(e)$ is the cross-sectional area of edge $e$, and $\sigma (e)$ is the surface which separates $\Omega (t_{\bullet }(e))$ and $\Omega (t^{\bullet }(e))$).

If edge based function $f:{\mathcal {E}}\to \mathbb {R}$ is an approximation of a vector field $\mathbf {F} _{\text{continuous}}:{\bigg (}\bigcup _{n\in {\mathcal {N}}}\Omega (n){\bigg )}\to \mathbb {R} ^{3}$, where $\forall e\in {\mathcal {E}}:f(e)=\mathbf {F} _{\text{continuous}}(\mathbf {q} (t_{\bullet }(e)))\cdot {\frac {\mathbf {q} (t^{\bullet }(e))-\mathbf {q} (t_{\bullet }(e))}{|\mathbf {q} (t^{\bullet }(e))-\mathbf {q} (t_{\bullet }(e))|}}$, then $\int _{e\in S}f(e)$ is an approximation of the surface integral $\iint _{\mathbf {q} \in \bigcup _{(e,s)\in S}s\sigma (e)}\mathbf {F} (\mathbf {q} )\cdot \mathbf {dS}$.

Recall that $S$ can be described via the edge based function $\delta _{2}(e;S)=\left\{{\begin{array}{cc}+1/l(e)&((e,+1)\in S)\\-1/l(e)&((e,-1)\in S)\\0&({\text{else}})\\\end{array}}\right.$. This allows the expression $\int _{e\in S}f(e)=\sum _{e\in {\mathcal {E}}}f(e)\delta _{2}(e;S)A(e)l(e)$. If edge based function $g$ denotes a multi-surface, then $\int _{e\in g}f(e)=\sum _{e\in {\mathcal {E}}}f(e)g(e)A(e)l(e)=\Sigma _{\mathcal {E}}(f,g)$.

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

Discrete analog to the Gradient

Given a node based function $f$, the "gradient" of $f$ is the edge based function $\nabla f$. For all $e\in {\mathcal {E}}$, $\nabla f|_{e}={\frac {f(t^{\bullet }(e))-f(t_{\bullet }(e))}{l(e)}}$. The gradient denotes the rate of change in $f$ along each edge $e$ in the preferred direction.

The path integral along path $P$ of the gradient of $f$ is the difference in $f$ between the endpoints. The value of $f$ at each node is shown, and along the path the value of $s_{i}l(e_{i})\nabla f|_{e_{i}}$ is shown at each edge. (the edge orientations are not shown for simplicity)

When given an arbitrary path $P=\langle (e_{1},s_{1}),(e_{2},s_{2}),...,(e_{k},s_{k})\rangle$, the path integral of the gradient along $P$ is:

$\int _{e\in P}\nabla f|_{e}=\sum _{i=1}^{k}s_{i}l(e_{i})\nabla f|_{e_{i}}$ $=\sum _{i=1}^{k}s_{i}l(e_{i}){\frac {f(t^{\bullet }(e_{i}))-f(t_{\bullet }(e_{i}))}{l(e_{i})}}$ $=\sum _{i=1}^{k}s_{i}(f(t^{\bullet }(e_{i}))-f(t_{\bullet }(e_{i})))$

The condition that $P$ is unbroken means that:

$\forall i\in \{1,2,\dots ,k-1\}:\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.$

$n_{1},n_{2},\dots ,n_{k+1}$ will denote the sequence of nodes visited by $P$:

$n_{1}=\left\{{\begin{array}{cc}t_{\bullet }(e_{1})&(s_{1}=+1)\\t^{\bullet }(e_{1})&(s_{1}=-1)\end{array}}\right.$

For each $i\in \{2,3,\dots ,k\}$, it is the case that $n_{i}=\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.$

$n_{k+1}=\left\{{\begin{array}{cc}t^{\bullet }(e_{k})&(s_{k}=+1)\\t_{\bullet }(e_{k})&(s_{k}=-1)\end{array}}\right.$

hence:

$\int _{e\in P}\nabla f|_{e}=\sum _{i=1}^{k}\left\{{\begin{array}{cc}f(t^{\bullet }(e_{i}))-f(t_{\bullet }(e_{i}))&(s_{i}=+1)\\f(t_{\bullet }(e_{i}))-f(t^{\bullet }(e_{i}))&(s_{i}=-1)\end{array}}\right.$