# Representing a Stencil of a Finite Difference Operator with Mathematica's Graphics3D

I have the following finite difference operator: $$Lu_{ijk}:= du_{ijk} +c(u_{i-1,j,k} + u_{i+1,j,k} + u_{i,j-1,k} + u_{i,j+1,k} + u_{i,j,k-1} + u_{i,j,k+1})\\ -u_{i-1,j+1,k}-u_{i-1,j-1,k} - u_{i-1,j,k-1} - u_{i-1,j,k+1} \\ -u_{i+1,j+1,k}-u_{i+1,j-1,k} - u_{i+1,j,k-1} - u_{i+1,j,k+1} \\ - u_{i,j+1,k-1} - u_{i,j-1,k-1} - u_{i, j+1,k+1} - u_{i,j-1,k+1},$$ where $d$ and $c$ are known constants. This is way too complex of an expression to expect people to understand, so I want to represent it by a stencil:

For example, I'd put a sphere of one color (representing $du_{ijk}$) at the origin, a sphere of another color at $(-1,0,0)$, (representing $cu_{i-1,j,k}$), and so on, and finally spheres of a final color at $(-1,1,0)$ (representing $-u_{i-1,j+1,k}$). My preliminary Mathematica commands are:

Graphics3D[{Specularity[White, 10], Red, Sphere[{0, 0, 0}],
Blue, Sphere[{0, 0, 1}], Sphere[{1, 0, 0}], Sphere[{-1, 0, 0}],
Sphere[{0, -1, 0}], Sphere[{0, 0, -1}],
Green, Sphere[{1, 1, 0}], Sphere[{1, 0, 1}]}, Boxed -> False]


but these spheres are way too big and cover each other up. How can this be fixed?

-
add comment

## 3 Answers

Take it in steps:

1. Extract the coefficients and locations into an appropriate data structure.

2. Use that data structure to create the graphics.

By examining the FullForm of the original expression, we can cobble a rule to find the key data: the coefficients $c$, $d$, and $-1$ and the offsets to the indexes. First, the expression itself:

s = Subscript;
exp = d s[u, i, j, k] +
c (s[u, i - 1, j, k] + s[u, i + 1, j, k] + s[u, i, j - 1, k] +
s[u, i, j + 1, k] + s[u, i, j, k - 1] + s[u, i, j, k + 1]) -
s[u, i - 1, j + 1, k] - s[u, i - 1, j - 1, k] - s[u, i - 1, j, k - 1] - s[u, i - 1, j, k + 1] -
s[u, i + 1, j + 1, k] - s[u, i + 1, j - 1, k] - s[u, i + 1, j, k - 1] - s[u, i + 1, j, k + 1] -
s[u, i, j + 1, k - 1] - s[u, i, j - 1, k - 1] -
s[u, i, j + 1, k + 1] - s[u, i, j - 1, k + 1]


Expand pairs the coefficients with the subscripts and Cases extracts the essential information:

data = Cases[Expand@exp, Times[c_, s[u, i_, j_, k_]] :> {c, i, j, k}];
data // MatrixForm


$$\left( \begin{array}{cccc} -1 & -1+i & -1+j & k \\ -1 & -1+i & j & -1+k \\ c & -1+i & j & k \\ \cdots \\ -1 & 1+i & 1+j & k \end{array} \right)$$

We are really interested in the offsets to the central index $(i,j,k)$, so one more step to extract them (via Replace this time) will be helpful. After doing it, let's group the offsets by common coefficient using GatherBy:

spec = GatherBy[{First@#,
Replace[Rest@#,  {Plus[x_?NumberQ, i_] ->  x, x_Symbol -> 0}, 1]} & /@ data, First]


To illustrate, here is what the first few elements of the first entry in spec look like:

$$\left( \begin{array}{cc} -1 & \{-1,-1,0\} \\ -1 & \{-1,0,-1\} \\ -1 & \{-1,0,1\} \\ \cdots \\ -1 & \{1,1,0\} \end{array} \right)$$

(You might be happier just entering the data in this format, or something close to it, at the outset: it's easier than entering all those subscripts.)

Choose some colors:

colors = Array[Hue[# / Length@spec, .8, .8] &, Length@spec];


The rest is easy. Let's make sure to include some visual cues such as thin lines connecting the base point to its neighbors, for otherwise this will look only like a random jumble of balls.

Graphics3D[ {
Table[{Specularity[White, 10],
GrayLevel[0.7], Tube[{{0, 0, 0}, Last@#}, 0.025] & /@ spec[[i]],
colors[[i]], Sphere[Last@#, .2] & /@ spec[[i]]}, {i, 1, Length@spec}]},
Boxed -> False, Axes -> True, AxesLabel -> {"i", "j", "k"} ]


(I leave the creation of a color key as an exercise :-).)

-
If I can nitpick... I prefer horizontal and vertical connections so that it is clearer which point is located where w.r.t. the center. –  rubenvb Jul 10 at 8:30
@rubenvb That's a great idea. How would you implement it? In this diagram, renderings of such connections would all run into one another. I can think of some possibilities, such as introducing a light 3D grid, or a set of nested transparent cubes, or other such visual references, but perhaps you have a specific form of visualization in mind? –  whuber Jul 10 at 13:10
no idea really, I am terrible at drawing with Mathematica. I don't see how they would run into each other though. Something like this comes to mind though. –  rubenvb Jul 10 at 13:56
@rubenvb Thank you: I see what you mean. There are arbitrary elements to that diagram--it could be drawn in many equivalent ways. However, something like it could be implemented in Mathematica, provided we settle on rules to determine which sets of line segments to draw. –  whuber Jul 10 at 14:06
add comment

Although you haven't exactly asked this, you might like to generate your graphic automatically by applying pattern matching on your difference operator. The basic idea is as below:

Clear[i, j, k];

op = Plus @@
MapThread[Subscript[u, i - #1, j - #2, k - #3] &,
RotateRight[{0, -1, 1, 0, 0, 0, 0}, #] & /@ {0, 2, 4}]


Giving $op = u_{i-1,j,k}+u_{i,j-1,k}+u_{i,j,k-1}+u_{i,j,k}+u_{i,j,k+1}+u_{i,j+1,k}+u_{i+1,j,k}$

Graphics3D[
Level[op /.
Subscript[_, i_, j_, k_] :> Sphere[{i, j, k}, .1] /. {i -> 0,
j -> 0, k -> 0}, {-3}] /. {s : Sphere[{0, 0, 0}, _] :>
Sequence[Red, s], s_Sphere :> Sequence[Blue, s]}, Axes -> True]


Disclaimer: I'm a fairly new and rather unsophisticated Mathematica user, so my code is probably crap and likely to break if you breathe on it too hard, but anyway, there you have it.

-
add comment

Not sure if this is what you're after. Anyway:

r = .1;
Graphics3D[{Specularity[White, 10],
Red, Sphere[{0, 0, 0}, r],
Blue, Sphere[{0, 0, 1}, r], Sphere[{1, 0, 0}, r], Sphere[{-1, 0, 0}, r],
Sphere[{0, -1, 0}, r], Sphere[{0, 0, -1}, r],
Green, Sphere[{1, 1, 0}, r], Sphere[{1, 0, 1}, r]}, Boxed -> False,
Axes -> True]


-
Oh, that was really easy! Thanks! –  Nick Thompson May 1 '13 at 23:15
add comment