A nice way to visualize a solution of the heat equation on a three dimensional cube

Having nothing to do on the train today, I fired up Mathematica with the intention to solve a differential equation on a three dimensional domain just to see how well it worked. It turned out to be easy as π. Just type in the equations. And better yet... With the slightest amount of Mathematica magic, I was able to visualize the solution for a fixed time.

Since we are not four dimensional (or more) beings, it is hard four us to look inside three dimensional objects. Unless, of course, they have some kind of openings. An easy way to open up a cube is to make it into lots of small cubes, like this.

It is just a lot of Cuboid's. The color function which I called heatColor is essentially a wrapper around RGBColor. It has two parameters min and max, for temperatures T close to min we choose litte red and a lot of blue. Vice versa for T close to max.

In order to solve the partial differential equation I use NDSolve. The boundary conditions are zero on all sides of the cube except the bottom (z=0). I just chose something that would be consistent with the other boundary conditions as well as with the inital condition. The inital condition is also chosen as anything simple that will be zero when x, y or z get close to zero or one. The result from NDSolve comes in the form of an interpolated function, so you can evaluate the numerical solution in any point that you wish!

You can play around with the number of cubes and the opacity to find your personal optimum.

The code:

solution =
NDSolve[{D[u[x, y, z, t], t] ==
D[u[x, y, z, t], x, x] + D[u[x, y, z, t], y, y] +
D[u[x, y, z, t], z, z],
u[0, y, z, t] == 0,
u[1, y, z, t] == 0,
u[x, 0, z, t] == 0,
u[x, 1, z, t] == 0,
u[x, y, 0, t] == t x y (1 - x) (1 - y),
u[x, y, 1, t] == 0,
u[x, y, z, 0] == x y z (1 - x) (1 - y) (1 - z)},
u,
{x, 0, 1}, {y, 0, 1}, {z, 0, 1}, {t, 0, 2}]


(* Let v be the first - and only - solution *)
v = (u /. solution)[[1]]


heatColor[T_, min_, max_] :=
With[{c = (T - min)/(max - min)}, RGBColor[{c, 0, 1 - c}]]


Graphics3D[Table[
{heatColor[v[i/5, j/5, k/5, 1], 0, 0.03],
Opacity[.50],
Cuboid[{i, j, k}/5, {i, j, k}/5 + .10]},
{i, 0, 5}, {j, 0, 5}, {k, 0, 5}], Axes -> True,
AxesLabel -> {x, y, z}]

Programming for multi-core processors

CPUs come with an increasingly number of cores. As a programmer, you just have to relate to this in some way. One way would be to make sure you really understand parallel programming to utilize the higher responsiveness and speedup made possible with multi-core CPUs. But parallel programming is complicated and it's quite common not see the speedup you hoped for. Here's an insight I just had:

If you parallelize your code yourself, you will also have to maintain and improve it yourself. Conversely, if you rely on automatic programming constructs, the provider will continue to improve the parallelization for you!

For instance, by using the Cilk extension of the C programming language to expose parallelism in your program, you can recompile it with an improved version of Cilk to get a better speedup. If you would have used for instance MPI, you would have to reprogram your code by yourself to get a faster program. We both know it; it's not gonna happen, it's way too boring to deal with old code :)

With this argument in mind, I think the only reason to learn about parallel programming is because it's fun. It's not necessary in order to remain competitive in a multi-core future. I guess this calls for a blog post on different ways to exploit multi-cores CPUs without having to do the actual parallelization yourself...

Bézier Curves in Mathematica

Picking up Applied Geometry for Computer Graphics and CAD by Duncan March, I read a chapter on Bézier curves. They are polynomial curves defined by control points that in a very intuitive way let a designer shape the curve. Bézier curves are really a killer app for the Manipulate and Locator commands in Mathematica. A cubic polynomial need four Locators (control points) in order to be defined:

You need a cubic polynomial in order to be able to create a loop or a cusp:

Is there any other interesting phenomena that becomes possible if we allow higer order polynomials? I don't know but here's a heart or the contour of a fox head. The LocatorAutoCreate option let us add Locators just by holding the alt key and click in the graphics area.

To do this in Mathematica (version 6 or more), define the Bernstein polynomials and the Bézier curve as on page 131 in March's book. Then in the Manipulate construction, you can make a list of Locators called b. Here's the code:

bernstein[i_, n_, t_] := n! (1 - t)^(n - i) t^i/((n - i)! i!)

bezier[t_, b_List] := 

 Module[{n = Length[b] - 1}, 

  Sum[b[[i + 1]] bernstein[i, n, t], {i, 0, n}]]

Manipulate[

 Module[

  {bez},

  bez[t_] = bezier[t, b];

  ListLinePlot[Table[bez[t], {t, 0, 1, 0.01}], 

   PlotRange -> {{0, 1}, {0, 5}}]

  ],

 {

  {b, Table[{i, 0.0}, {i, 0, 1, 1/3.0}]},

  Locator,

  LocatorAutoCreate -> True

  }

 ]

It is crucial that the bezier function takes an arbitrary length list b as parameter. Otherwise it wouldn't work to use LocatorAutoCreate. Consider doing this in Java, how many cups of coffee you would need. To finish, in version 7 you can use an in-built function called BezierCurve instead of defining your own. But it's not necessary from a performance perspective, the response is immediate with my version on a single-core 1.7 MHz Pentium M.