Scientific Computing: Interactive Modules

Michael Heath and company present a bunch of educational Java modules relevant to Scientific Computing. They are fun to play around with in a browser.

As an interesting side note, we used the book on which many of the modules are based, in our Algorithms 1 undergrad class.

Links

1. In defense of inefficient code: Mike Croucher makes a case for working-but-not-particularly-efficient code written by practitioners in high-level languages (Matlab, Mathematica, python). As he puts it:

It comes down to this. CPU time is cheap. Very cheap. Human time, particularly specialised human time, is expensive.

and also:

In my opinion, high level programming languages such as Mathematica, MATLAB and Python have democratised scientific programming. Now, almost anyone who can think logically can turn their scientific ideas into working code. I’ve seen people who have had no formal programming training at all whip up models, get results and move on with their research. Let’s be clear here – It’s results that matter not how you coded them.

I often use these high-level languages to rapidly prototype new ideas. If and when required, it can always be translated into C++ or Fortran. Also from a practical perspective inefficient code that runs 10 times slower than highly optimized code is acceptable, if it still takes only 10 seconds to run.

2. Geeky jokes on Tanya Khovanova's Blog:

I like this one:

I just learned that 4,416,237 people got married in the US in 2010. Not to nitpick, but shouldn’t it be an even number?

Obviously, you can understand the audience she caters to from this comment.

Why numerical differentiation may be trickier than you think?

Here is a link to a fascinating presentation by Harvey Stein ("Risky Measures of Risk: Error Analysis of Numerical Differentiation").

He makes a very "visual" case for why one needs to think carefully before using large (convexity error) or small step sizes (cancellation error).

"Publish"ing in Matlab

I use GNU Octave as the glue to post-process and make preliminary visualizations most of the time. Every once in a while, however, I need to use its commercial cousin, Matlab.

This Fall, I am teaching an undergrad numerical analysis class, and realized that Matlab has this interesting documenting feature called "Publish", which lets you use a double percent sign %% in addition to the regular % sign used to begin a comment. Using these %% signs, you can divide your code into "cells" which is helpful for selectively testing individual cells -- an added benefit.

Check out a video on how to use this feature at Mathworks.

By proper use of %% and % comments, one can actually produce a nice document (in html, LaTeX or some other options), which is ideal, since now students can submit only a single "M-file", and I can check whether the program works, and the accompanying documentation in one shot.

Dispersing points on the surface of a sphere

Dispersing or choosing points, uniformly or randomly, on the surface of a sphere is a task that someone like me who does molecular simulations, runs into, every once in a while. In some form or shape, this problem is encountered, for example, while modeling a 3D random walk, the motion of a Brownian particle in space, decorating a nanoparticle with stabilizers etc.

While the task is not hard, there are some elegant methods and potential pitfalls that a person doing this for the first time should be aware of.

As I mentioned earlier, one may choose points randomly, or uniformly on the surface (of a unit sphere, here, but can be generalized trivially).

1. Random:

The key point here is that it is incorrect to choose points by selecting angles phi and theta corresponding to spherical coordinates from uniform distributions. This incorrectly concentrates points near the poles.

The correct method is to select two random numbers u and v from a uniform distribution, and construct the two angles as follows:

u = rand();            // uniform random number between

v = rand();            // zero and one.

phi = 2 * pi * u;      // pi is 3.141...

theta = acos(2*v - 1)  // acos is inverse cosine

This is explained with figures here.

2. Uniform:

This is trickier than I thought when I first attempted to do it. Here is an old link (circa 1998, but still good!) which touches upon some of the subtleties. The trouble starts with what "uniform" means, for an arbitrary number of points to be distributed.

If we decide that uniform means a distribution which "maximizes the minimum separation" between n points on the sphere, we can start going somewhere.

There are a number of algorithms that can try to solve this problem, as described on this excellent page. Here is C++ code that I wrote while implementing the golden spiral rule from that page.

//
// Spray n points uniformly on the surface of a sphere 

// of radius "radius"
//

// The positions are returned as an array of Vectors

// where, struct Vector { double x; double y; double z }

//

// Based on algorithm from 

// http://www.xsi-blog.com/archives/115

//

void SprayPointsSphere(int n, double radius, Vector * p)
{
  double inc = PI * (3.0 - sqrt(5.0));
  double off = 2.0/n;

  for(int k = 0; k < n; k++) {
    double y   = k * off - 1 + (off/2);
    double r   = sqrt(1 - y*y);
    double phi = k * inc;

    p[k].x = cos(phi)*r * radius;

    p[k].y = y * radius;

    p[k].z = sin(phi)*r * radius);

  }
}