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

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

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

It is however, accurate to approximate ${\displaystyle \Delta y(i)}$ with ${\displaystyle \Delta z(i)=2x_{i}\Delta x}$ since the error relative to ${\displaystyle \Delta x}$ does vanish as ${\displaystyle \Delta x\to 0^{+}}$: it is the case that ${\displaystyle {\frac {\Delta y(i)-\Delta z(i)}{\Delta x}}=\Delta x\to 0}$ as ${\displaystyle \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 ${\displaystyle \Delta x}$ approach 0 as ${\displaystyle \Delta x\to 0}$ as opposed to the absolute error, the following example will be used:

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

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

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

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

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

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

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

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

An arbitrary function ${\displaystyle 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 ${\displaystyle \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: ${\displaystyle \{n_{1}\mapsto f(n_{1}),n_{2}\mapsto f(n_{2}),\dots \}}$

An arbitrary function ${\displaystyle 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 ${\displaystyle \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: ${\displaystyle \{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 ${\displaystyle {\mathcal {N}}=\{n_{1},n_{2},n_{3},n_{4},n_{5}\}}$ and the edge set is ${\displaystyle {\mathcal {E}}=\{e_{1},e_{2},e_{3},e_{4},e_{5},e_{6}\}}$. The origin and terminal nodes of each edge are: ${\displaystyle (t_{\bullet }(e_{1}),t^{\bullet }(e_{1}))=(n_{1},n_{2})}$; ${\displaystyle (t_{\bullet }(e_{2}),t^{\bullet }(e_{2}))=(n_{1},n_{3})}$; ${\displaystyle (t_{\bullet }(e_{3}),t^{\bullet }(e_{3}))=(n_{2},n_{4})}$; ${\displaystyle (t_{\bullet }(e_{4}),t^{\bullet }(e_{4}))=(n_{4},n_{3})}$; ${\displaystyle (t_{\bullet }(e_{5}),t^{\bullet }(e_{5}))=(n_{4},n_{5})}$; and ${\displaystyle (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 ${\displaystyle n\in {\mathcal {N}}}$ will correspond to both a point ${\displaystyle \mathbf {q} (n)}$ in space, and a volume ${\displaystyle \Omega (n)}$ that contains ${\displaystyle \mathbf {q} (n)}$. ${\displaystyle V(n)}$ will denote the volume of ${\displaystyle \Omega (n)}$, also referred to as the volume of node ${\displaystyle n}$. The volume of each node must be positive: ${\displaystyle V(n)>0}$ for all ${\displaystyle n\in {\mathcal {N}}}$. In the example directed graph, an example choice of volumes for each node is ${\displaystyle V=\{n_{1}\mapsto 3,n_{2}\mapsto 2,n_{3}\mapsto 2,n_{4}\mapsto 1,n_{5}\mapsto 4\}}$.

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

${\displaystyle l(e)}$ will denote the distance ${\displaystyle |\mathbf {q} (t^{\bullet }(e))-\mathbf {q} (t_{\bullet }(e))|}$ from ${\displaystyle \mathbf {q} (t_{\bullet }(e))}$ to ${\displaystyle \mathbf {q} (t^{\bullet }(e))}$, which will be referred to as the length of edge ${\displaystyle e}$. The length of each edge must be positive: ${\displaystyle l(e)>0}$ for all ${\displaystyle e\in {\mathcal {E}}}$. ${\displaystyle A(e)}$ will denote the area of ${\displaystyle \sigma (e)}$ projected onto a plane that is perpendicular to the displacement ${\displaystyle \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 ${\displaystyle e}$. The cross-sectional area of each edge must be positive: ${\displaystyle A(e)>0}$ for all ${\displaystyle e\in {\mathcal {E}}}$. In the example directed graph, an example choice of lengths for each edge is ${\displaystyle 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 ${\displaystyle 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 ${\displaystyle {\mathcal {G}}}$ will be characterized by:

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

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

Point ${\displaystyle n_{0}}$ can be denoted by the node based function ${\displaystyle \delta _{0}(n;n_{0}):{\mathcal {N}}\to \mathbb {R} }$ where for each node ${\displaystyle n\in {\mathcal {N}}}$ it is the case that ${\displaystyle \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 ${\displaystyle {\frac {1}{V(n)}}}$ arises from the fact that the point is "spread out" over the volume of ${\displaystyle n}$.

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

### Volumes and multi-volumes

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

Volume ${\displaystyle N}$ can be denoted by a node based membership function: ${\displaystyle \delta _{3}(n;N):{\mathcal {N}}\to \mathbb {R} }$ where for each node ${\displaystyle n\in {\mathcal {N}}}$, it is the case that ${\displaystyle \delta _{3}(n;N)=\left\{{\begin{array}{cc}1&(n\in N)\\0&(n\notin N)\\\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 ${\displaystyle 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:

${\displaystyle 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 ${\displaystyle \{n_{1}\},\{n_{3}\},\{n_{4},n_{5}\}}$ with the respective weights ${\displaystyle -2,4,0.5}$.

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

### Paths and multi-paths

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

A path ${\displaystyle P}$ can be denoted by the edge based function ${\displaystyle \delta _{1}(e;P):{\mathcal {E}}\to \mathbb {R} }$ where for each edge ${\displaystyle e\in {\mathcal {E}}}$, it is the case that ${\displaystyle \delta _{1}(e;P)=\sum _{i=1}^{k}s_{i}\delta _{1}(e;e_{i})}$ where ${\displaystyle \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 ${\displaystyle {\frac {1}{A(e)}}}$ arises from the fact that the path is "spread out" over the thickness of ${\displaystyle 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 ${\displaystyle A(e)=1}$ for all ${\displaystyle e\in {\mathcal {E}}}$, the edge based function ${\displaystyle 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:

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

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

### 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 ${\displaystyle S}$ is a collection of edges and orientations: ${\displaystyle S=\{(e_{1},s_{1}),(e_{2},s_{2}),\dots ,(e_{k},s_{k})\}}$ where all of the edges ${\displaystyle e_{1},e_{2},\dots ,e_{k}\in {\mathcal {E}}}$ are distinct and ${\displaystyle s_{i}\in \{+1,-1\}}$ is an indicator that describes the surface's orientation. Each edge ${\displaystyle e_{i}}$ is an edge that passes through the surface, and ${\displaystyle s_{i}}$ is ${\displaystyle +1}$ if ${\displaystyle e_{i}}$ passes through the surface in the preferred direction, and ${\displaystyle s_{i}}$ is ${\displaystyle -1}$ if ${\displaystyle e_{i}}$ passes through the surface in the opposite direction. The combined surface denoted by ${\displaystyle S}$ is ${\displaystyle \sigma (S)=\bigcup _{i=1}^{k}s_{i}\sigma (e_{i})}$ (surface ${\displaystyle -\sigma (e_{i})}$ is surface ${\displaystyle \sigma (e_{i})}$ with the orientation reversed) and the total area is ${\displaystyle A(S)=\sum _{i=1}^{k}A(e_{i})}$.

A surface ${\displaystyle S}$ can be denoted by the edge based function ${\displaystyle \delta _{2}(e;S):{\mathcal {E}}\to \mathbb {R} }$ where for each edge ${\displaystyle e\in {\mathcal {E}}}$, it is the case that ${\displaystyle \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 ${\displaystyle {\frac {1}{l(e)}}}$ arises from the fact that the surface is "spread out" over the length of ${\displaystyle 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 ${\displaystyle l(e)=1}$ for all ${\displaystyle e\in {\mathcal {E}}}$, the edge based function ${\displaystyle 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:

${\displaystyle 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 ${\displaystyle \{(e_{1},+1)\},\{(e_{2},+1)\},}$ and ${\displaystyle \{(e_{6},+1)\}}$ with the respective weights ${\displaystyle 2,2,-1}$.

${\displaystyle 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 ${\displaystyle \{(e_{1},+1),(e_{2},+1),(e_{6},-1)\}}$ and ${\displaystyle \{(e_{1},+1),(e_{2},+1)\}}$ with the respective weights ${\displaystyle 1}$ and ${\displaystyle 1}$.

### Model of Cartesian coordinates

The discrete model used for Cartesian coordinates will consist of an infinite 3 dimensional grid of nodes that spans ${\displaystyle \mathbb {R} ^{3}}$. To establish the model, the following is needed:

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

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

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

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

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

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

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

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

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

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

### Model of Cylindrical coordinates

The discrete model used for Cylindrical coordinates will consist of an infinite 3 dimensional cylindrical grid of nodes that spans ${\displaystyle \mathbb {R} ^{3}}$. To establish the model, the following is needed:

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

Each node ${\displaystyle n\in {\mathcal {N}}}$ will correspond to a unique triplet ${\displaystyle (i,j,k)}$ of integers, and will be denoted by ${\displaystyle n_{(i,j,k)}}$. The cylindrical coordinate denoted by ${\displaystyle n_{(i,j,k)}}$ will be ${\displaystyle \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 ${\displaystyle \phi }$ is cyclical), not every triplet ${\displaystyle (i,j,k)}$ will be allowed. The following three domains will be defined: ${\displaystyle I_{0}}$, ${\displaystyle I_{1}}$, and ${\displaystyle I_{2}}$.

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

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

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

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

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

There exists an edge ${\displaystyle e_{\phi ,(i,j,k)}}$ for each ${\displaystyle (i,j,k)\in I_{0}}$ such that ${\displaystyle t_{\bullet }(e_{\phi ,(i,j,k)})=n_{(i,j,k)}}$ and ${\displaystyle t^{\bullet }(e_{\phi ,(i,j,k)})=\left\{{\begin{array}{cc}n_{(i,j+1,k)}&(j

There exists an edge ${\displaystyle e_{z,(i,j,k)}}$ for each ${\displaystyle (i,j,k)\in I_{1}}$ such that ${\displaystyle t_{\bullet }(e_{z,(i,j,k)})=n_{(i,j,k)}}$ and ${\displaystyle 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 ${\displaystyle l_{\phi }(i_{\text{real}})=(i_{\text{real}}\Delta \rho )\Delta \phi }$ for any nonnegative real number ${\displaystyle i_{\text{real}}}$.

For each ${\displaystyle (i,j,k)\in I_{1}}$ the cylindrical coordinate that corresponds to ${\displaystyle n_{(i,j,k)}}$ is ${\displaystyle \mathbf {q} (n_{(i,j,k)})=(i\Delta \rho ,j\Delta \phi ,k\Delta z)}$, and the space (described using cylindrical coordinates) associated with ${\displaystyle n_{(i,j,k)}}$ is ${\displaystyle \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 ${\displaystyle (i,j,k)\in I_{1}}$ the volume of node ${\displaystyle n_{(i,j,k)}}$ is ${\displaystyle 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 ${\displaystyle (i,j,k)\in I_{2}}$ the length and area of edge ${\displaystyle e_{\rho ,(i,j,k)}}$ is ${\displaystyle l(e_{\rho ,(i,j,k)})=\Delta \rho }$ and ${\displaystyle A(e_{\rho ,(i,j,k)})\approx l_{\phi }(i+{\frac {1}{2}})\cdot \Delta z}$ respectively.

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

For each ${\displaystyle (i,j,k)\in I_{1}}$ the length and area of edge ${\displaystyle e_{z,(i,j,k)}}$ is ${\displaystyle l(e_{z,(i,j,k)})=\Delta z}$ and ${\displaystyle 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 ${\displaystyle f(\rho ,\phi ,z)}$ over cylindrical coordinates, ${\displaystyle f}$ will be approximated by the node based function ${\displaystyle f_{\text{approx}}:{\mathcal {N}}\to \mathbb {R} }$. Function ${\displaystyle f_{\text{approx}}}$ is defined by ${\displaystyle f_{\text{approx}}(n_{(i,j,k)})=f(i\Delta \rho ,j\Delta \phi ,k\Delta z)}$ for each ${\displaystyle (i,j,k)\in I_{1}}$.

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

${\displaystyle \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.}$

${\displaystyle \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)}$

${\displaystyle \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

The discrete model used for Spherical coordinates will consist of an infinite 3 dimensional spherical grid of nodes that spans ${\displaystyle \mathbb {R} ^{3}}$. To establish the model, the following is needed:

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

Each node ${\displaystyle n\in {\mathcal {N}}}$ will correspond to a unique triplet ${\displaystyle (i,j,k)}$ of integers, and will be denoted by ${\displaystyle n_{(i,j,k)}}$. The spherical coordinate denoted by ${\displaystyle n_{(i,j,k)}}$ will be ${\displaystyle \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 ${\displaystyle \phi }$ is cyclical), not every triplet ${\displaystyle (i,j,k)}$ will be allowed. The following three domains will be defined: ${\displaystyle I_{n}}$, ${\displaystyle I_{r}}$, ${\displaystyle I_{\theta }}$, and ${\displaystyle I_{\phi }}$.

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

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

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

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

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

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

There exists an edge ${\displaystyle e_{\theta ,(i,j,k)}}$ for each ${\displaystyle (i,j,k)\in I_{\theta }}$ such that ${\displaystyle 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 ${\displaystyle t^{\bullet }(e_{\theta ,(i,j,k)})=\left\{{\begin{array}{cc}n_{(i,j+1,k)}&(j

There exists an edge ${\displaystyle e_{\phi ,(i,j,k)}}$ for each ${\displaystyle (i,j,k)\in I_{\phi }}$ such that ${\displaystyle t_{\bullet }(e_{\phi ,(i,j,k)})=n_{(i,j,k)}}$ and ${\displaystyle t^{\bullet }(e_{\phi ,(i,j,k)})=\left\{{\begin{array}{cc}n_{(i,j,k+1)}&(k

Next, the volume of each node and the lengths and area of each edge will be calculated. To simplify notation, let ${\displaystyle l_{\theta }(i_{\text{real}})=(i_{\text{real}}\Delta r)\Delta \theta }$ for any nonnegative real number ${\displaystyle i_{\text{real}}}$, and let ${\displaystyle 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 ${\displaystyle i_{\text{real}}}$ and any real number ${\displaystyle j_{\text{real}}\in [0,N_{\theta }]}$.

For each ${\displaystyle (i,j,k)\in I_{n}}$ the spherical coordinate that corresponds to ${\displaystyle n_{(i,j,k)}}$ is ${\displaystyle \mathbf {q} (n_{(i,j,k)})=(i\Delta r,j\Delta \theta ,k\Delta \phi )}$, and the space (described using spherical coordinates) associated with ${\displaystyle n_{(i,j,k)}}$ is ${\displaystyle \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}}\;00\;{\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 ${\displaystyle (i,j,k)\in I_{n}}$ the volume of node ${\displaystyle n_{(i,j,k)}}$ is ${\displaystyle 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}}\;00\;{\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 ${\displaystyle (i,j,k)\in I_{r}}$ the length and area of edge ${\displaystyle e_{r,(i,j,k)}}$ is ${\displaystyle l(e_{r,(i,j,k)})=\Delta r}$ and ${\displaystyle 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 respectively.

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

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

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

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

${\displaystyle \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.}$

${\displaystyle \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.}$

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

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

• ${\displaystyle I=(i,j,k)}$ denotes an arbitrary triplet of integer indices from ${\displaystyle \mathbb {Z} ^{3}}$.
• For each ${\displaystyle c=1,2,3}$: ${\displaystyle 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.}$
• ${\displaystyle \Xi =(\xi _{1},\xi _{2},\xi _{3})}$ denotes an arbitrary curvilinear coordinate ${\displaystyle (\xi _{1},\xi _{2},\xi _{3})}$.
• For each ${\displaystyle c=1,2,3}$: ${\displaystyle \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.}$
• ${\displaystyle \Delta \Xi =(\Delta \xi _{1},\Delta \xi _{2},\Delta \xi _{3})}$ denotes an arbitrary triplet of differences in the coordinates ${\displaystyle \xi _{1},\xi _{2},\xi _{3}}$.
• ${\displaystyle I\Delta \Xi =(i\Delta \xi _{1},j\Delta \xi _{2},k\Delta \xi _{3})}$

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

For the purposes of simplicity, it will be assumed that the coordinates ${\displaystyle \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 ${\displaystyle {\mathcal {N}}}$ connected by directed edges ${\displaystyle {\mathcal {E}}}$ will be established:

Start by choosing "resolutions" for each coordinate ${\displaystyle \xi _{1},\xi _{2},\xi _{3}}$ denoted by ${\displaystyle \Delta \xi _{1},\Delta \xi _{2},\Delta \xi _{3}}$. The quantities ${\displaystyle \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 ${\displaystyle I\in \mathbb {Z} ^{3}}$, a node ${\displaystyle n_{I}\in {\mathcal {N}}}$ will be created at position ${\displaystyle \mathbf {q} (n_{I})=\mathbf {q} (I\Delta \Xi )}$.

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

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

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

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

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

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

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

## Discrete analog to Integrals

### Volume Integrals

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

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

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

### Path Integrals

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

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

Recall that ${\displaystyle P}$ can be described via the edge based function ${\displaystyle \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 ${\displaystyle \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 ${\displaystyle g}$ denotes a multi-path, then ${\displaystyle \int _{e\in g}f(e)=\sum _{e\in {\mathcal {E}}}f(e)g(e)A(e)l(e)}$.

### Surface Integrals

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

If edge based function ${\displaystyle f:{\mathcal {E}}\to \mathbb {R} }$ is an approximation of a vector field ${\displaystyle \mathbf {F} _{\text{continuous}}:{\bigg (}\bigcup _{n\in {\mathcal {N}}}\Omega (n){\bigg )}\to \mathbb {R} ^{3}}$, where ${\displaystyle \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 ${\displaystyle \int _{e\in S}f(e)}$ is an approximation of the surface integral ${\displaystyle \iint _{\mathbf {q} \in \bigcup _{(e,s)\in S}s\sigma (e)}\mathbf {F} (\mathbf {q} )\cdot \mathbf {dS} }$.

Recall that ${\displaystyle S}$ can be described via the edge based function ${\displaystyle \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 ${\displaystyle \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 ${\displaystyle g}$ denotes a multi-surface, then ${\displaystyle \int _{e\in g}f(e)=\sum _{e\in {\mathcal {E}}}f(e)g(e)A(e)l(e)}$.

## Discrete analog to the Gradient

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

The path integral along path ${\displaystyle P}$ of the gradient of ${\displaystyle f}$ is the difference in ${\displaystyle f}$ between the endpoints. The value of ${\displaystyle f}$ at each node is shown, and along the path the value of ${\displaystyle 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 ${\displaystyle P=\langle (e_{1},s_{1}),(e_{2},s_{2}),...,(e_{k},s_{k})\rangle }$, the path integral of the gradient along ${\displaystyle P}$ is:

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

The condition that ${\displaystyle P}$ is unbroken means that:

${\displaystyle \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.}$

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

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

${\displaystyle 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:

${\displaystyle \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.}$ ${\displaystyle =\sum _{i=1}^{k}\left\{{\begin{array}{cc}f(n_{i+1})-f(n_{i})&(s_{i}=+1)\\f(n_{i+1})-f(n_{i})&(s_{i}=-1)\end{array}}\right.}$ ${\displaystyle =\sum _{i=1}^{k}(f(n_{i+1})-f(n_{i}))}$ ${\displaystyle =f(n_{k+1})-f(n_{1})}$

The path integral of the gradient of ${\displaystyle f}$ along ${\displaystyle P}$ is simply the difference of ${\displaystyle f}$ between the final node and the starting node: ${\displaystyle \int _{e\in P}\nabla f|_{e}=f(n_{k+1})-f(n_{1})}$.

Given an arbitrary multi-volume denoted by node based function ${\displaystyle f}$, the multi-surface which is the outwards oriented surface of the multi-volume is denoted by the edge based function ${\displaystyle g(e)=-\nabla f|_{e}}$. When traversing an arbitrary edge, the net number of outwards transitions through surfaces is the negative change in ${\displaystyle f}$.

The gradient can also be used to compute the surfaces of volumes:

Consider an arbitrary volume ${\displaystyle U\subseteq {\mathcal {N}}}$. The outwards oriented surface ${\displaystyle S}$ of ${\displaystyle U}$ is defined by: for every directed edge ${\displaystyle e\in {\mathcal {E}}}$, it is the case that ${\displaystyle (e,+1)\in S\iff (t_{\bullet }(e)\in U\;{\text{and}}\;t^{\bullet }(e)\notin U)}$ and ${\displaystyle (e,-1)\in S\iff (t_{\bullet }(e)\notin U\;{\text{and}}\;t^{\bullet }(e)\in U)}$.

Volume ${\displaystyle U}$ can be denoted by the node based function ${\displaystyle \delta _{3}(n;U)=\left\{{\begin{array}{cc}1&(n\in U)\\0&(n\notin U)\end{array}}\right.}$, and the outwards oriented surface ${\displaystyle S}$ can be denoted by the edge based function ${\displaystyle \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 gradient of ${\displaystyle \delta _{3}(n;U)}$ is:

${\displaystyle \nabla \delta _{3}(n;U)|_{e}={\frac {\delta _{3}(t^{\bullet }(e);U)-\delta _{3}(t_{\bullet }(e);U)}{l(e)}}}$ ${\displaystyle ={\frac {1}{l(e)}}\left\{{\begin{array}{cc}1&(t^{\bullet }(e)\in U\;{\text{and}}\;t_{\bullet }(e)\notin U)\\-1&(t^{\bullet }(e)\notin U\;{\text{and}}\;t_{\bullet }(e)\in U)\\0&({\text{else}})\end{array}}\right.}$ ${\displaystyle ={\frac {1}{l(e)}}\left\{{\begin{array}{cc}1&((e,-1)\in S)\\-1&((e,+1)\in S)\\0&({\text{else}})\end{array}}\right.}$ ${\displaystyle =-\left\{{\begin{array}{cc}+1/l(e)&((e,+1)\in S)\\-1/l(e)&((e,-1)\in S)\\0&({\text{else}})\end{array}}\right.}$ ${\displaystyle =-\delta _{2}(e;S)}$

Therefore the outwards oriented surface or multi-surface of a volume or multi-volume denoted by an arbitrary node based function ${\displaystyle f}$ is denoted by the edge based function ${\displaystyle -\nabla f}$.

The following sections will demonstrate that the gradient defined for node based functions is consistent with the gradient for scalar fields in ${\displaystyle \mathbb {R} ^{3}}$. 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 discrete model of curvilinear coordinates. The same shorthand notation will be used.

Starting with an arbitrary scalar field ${\displaystyle f(\Xi )}$, the node based function that approximates ${\displaystyle f}$ is ${\displaystyle f_{\text{approx}}(n_{I})=f(I\Delta \Xi )}$ for all ${\displaystyle I\in \mathbb {Z} ^{3}}$.

For each ${\displaystyle c=1,2,3}$ and ${\displaystyle I\in \mathbb {Z} ^{3}}$, the gradient along edge ${\displaystyle e_{\xi _{c},I}}$ is:

${\displaystyle \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}))}{l(e_{\xi _{c},I})}}}$ ${\displaystyle ={\frac {f_{\text{approx}}(n_{I+1_{c}})-f_{\text{approx}}(n_{I})}{\Delta \xi _{c}|\mathbf {u} _{\xi _{c}}(I\Delta \Xi )|}}}$ ${\displaystyle ={\frac {f(I\Delta \Xi +\Delta \xi _{c})-f(I\Delta \Xi )}{\Delta \xi _{c}|\mathbf {u} _{\xi _{c}}(I\Delta \Xi )|}}}$

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

The gradient of scalar field ${\displaystyle f}$ is a vector quantity with components that are multiples of the basis vectors ${\displaystyle \mathbf {v} _{\xi _{1}},\mathbf {v} _{\xi _{2}},\mathbf {v} _{\xi _{3}}}$. The coefficient of ${\displaystyle \mathbf {v} _{\xi _{1}}}$ is the gradient of ${\displaystyle f_{\text{approx}}}$ along edges that are parallel to ${\displaystyle \mathbf {v} _{\xi _{1}}}$, and similarly for ${\displaystyle \mathbf {v} _{\xi _{2}}}$ and ${\displaystyle \mathbf {v} _{\xi _{3}}}$. The fact that ${\displaystyle \mathbf {v} _{\xi _{1}}}$, ${\displaystyle \mathbf {v} _{\xi _{2}}}$, and ${\displaystyle \mathbf {v} _{\xi _{3}}}$ are of unit length and mutually orthogonal (perpendicular) is important.

The gradient of scalar field ${\displaystyle f}$ at ${\displaystyle \Xi }$ can now be approximated by the following:

${\displaystyle \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 )}$ ${\displaystyle =\sum _{c=1}^{3}{\frac {f(I(\Xi )\Delta \Xi +\Delta \xi _{c})-f(I(\Xi )\Delta \Xi )}{\Delta \xi _{c}|\mathbf {u} _{\xi _{c}}(I(\Xi )\Delta \Xi )|}}\mathbf {v} _{\xi _{c}}(I(\Xi )\Delta \Xi )}$ ${\displaystyle \approx \sum _{c=1}^{3}{\frac {f(\Xi +\Delta \xi _{c})-f(\Xi )}{\Delta \xi _{c}|\mathbf {u} _{\xi _{c}}(\Xi )|}}\mathbf {v} _{\xi _{c}}(\Xi )}$

As ${\displaystyle \Delta \xi _{1},\Delta \xi _{2},\Delta \xi _{3}\to 0^{+}}$, it becomes apparent that ${\displaystyle \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 )}$.

#### The Gradient in Cartesian Coordinates

Applying the curvilinear coordinate model to Cartesian coordinates means that ${\displaystyle \xi _{1}=x}$, ${\displaystyle \xi _{2}=y}$, and ${\displaystyle \xi _{3}=z}$. It is also the case that ${\displaystyle \mathbf {u} _{x}(x,y,z)=\mathbf {i} }$, ${\displaystyle \mathbf {u} _{y}(x,y,z)=\mathbf {j} }$, and ${\displaystyle \mathbf {u} _{z}(x,y,z)=\mathbf {k} }$. Applying the formula for the gradient that was derived from the discrete model of curvilinear coordinates gives:

${\displaystyle \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} }$ ${\displaystyle ={\frac {\partial f}{\partial x}}\mathbf {i} +{\frac {\partial f}{\partial y}}\mathbf {j} +{\frac {\partial f}{\partial z}}\mathbf {k} }$

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 ${\displaystyle \xi _{1}=\rho }$, ${\displaystyle \xi _{2}=\phi }$, and ${\displaystyle \xi _{3}=z}$. It is also the case that ${\displaystyle \mathbf {u} _{\rho }(\rho ,\phi ,z)={\hat {\mathbf {\rho } }}}$, ${\displaystyle \mathbf {u} _{\phi }(\rho ,\phi ,z)=\rho {\hat {\mathbf {\phi } }}}$, and