# Still Chugging Along

I apologize for not posting anything new in a while.  Most all of my time this summer has went into studying various numerical algorithms for computing differential equations.  I’ve written programs that model all kinds of things, and I’ve been testing how quickly various iterative methods converge and how long it takes my computer to solve matrices of various sizes using say LU Decomposition.   MASSIVE matrices.

It’s been a really productive time for me.  I’ve been practicing breaking partial differential equations into finite elements and solving them numerically, either by putting it all in a giant matrix and letting it crunch away, or using other iterative methods, such as say the Gauss-Seidal method.  There’s different techniques and I’ve been seeing how well they all work.  The current method I’m testing is the Crank-Nicholson method.

What is this sort of stuff used for?  I’ll give you an example.  Say you’re modeling how heat flows through, I dunno, a hubcap, and we want to watch the temperature in the material over time at each location.  I could model all of that using these techniques.

I’d watch that nice 3D model on screen and get a good feeling for what’s going on.  Temperature is important for what I’m modeling because different frequencies of radiation are emitted by objects at different temperatures.   I want to model everyday things like a metal hubcap.  Temperature represents the vibrations and motions of all the tiny atoms making up that hubcap.  The hubcap we actually see with our eyes is based on all the tiny atoms on its surface vibrating, which causes slight distortions in electric and magnetic fields which oscillate their way through space toward your eyeballs in small packets of energy called photons.

I want to be able to zoom in to a section of that hubcap’s surface and watch the individual atoms vibrating, observing the progression and evolution of those electric and magnetic fields as the radiation is emitted.

This may sound like an incredibly boring way to spend your summer, but these are the sorts of things I enjoy doing .  The world we live is so complex, yet most us have no idea what we’re looking at.  What is a moment in time?  As you’re sitting there in your computer chair, what sorts of processes created that moment for you?  I’m curious to find out.  I’m not being paid for this.  I’m not publishing any papers.  I just do it because I want to figure this stuff out.  I like to simulate complex things and render them in 3D on my computer screen and watch how the system behaves over time.  I guess for right now, my primary interest is in optical phenomenon, particularly how light scatters and moves through surfaces.   This ties in with my interest in what space is.  The world I see around me is based on objects scattering light from their surfaces.

Computers are amazing.  If it wasn’t for them, there’d be no way to visualize these sorts of things.  Not in any real complexity.  I’m hoping to get to the place where I can model a handful of atoms on the surface of that hubcap, simulating everything that’s going on in full quantum mechanical detail.  I want to view the wave function of those atoms, watching them get excited and then discharge light.   Computers can translate those mathematical equations, which look a lot like hieroglyphics (even to someone like me who spends most of his time looking at them), and make them into 3D pictures like you see above.  That’s neat.  Takes a long time to master all of this, though.  You have to be a master in computer programming, writing the 3D simulation in say OpenGL, understanding all the mathematics and physics behind what you’re modeling.  Not only that, you have to understand numerical techniques which are different from traditional physics.  You have to be able to break the problem into small pieces and let the computer chug through the problem.  I won’t lie, this stuff isn’t easy, but it’s fulfilling to work on.

With every simulator I write, some complex, ugly equation becomes a vivid picture in my mind.  I deeply understand the physical process each equation represents, which is my ultimate goal.

Just yesterday I was with my mom and dad watching old video tapes.  I saw myself back in 1996, when I was only 12 or 13 years old.   I thought about the kinds of things which excited me back then and compared them to me today.   I sort of laughed to myself reflecting on it.  I’m a totally different person.  To give you an example, the other day I was grinning ear to ear, filled with joy because I learned something I had never even thought of doing.   I came across this in a numerical programming textbook.

“One problem is that many algorithms naturally like to go from 1 to M, not from 0 to M − 1. Sure, you can always convert them, but they then often acquire a baggage of additional arithmetic in array indices that is, at best, distracting.”

I immediately thought, “YES!  I HATE that.  In C and C++, indexed arrays start with 0, not 1, and it’s such a pain dealing with that.  You have say a 3 x 4 matrix.  You want to access the first element.  It’s not [1,1], it’s, [0,0].  *pulls hair out*.

“Consider:

float b[4],*bb;
bb=b-1;

The pointer bb now points one location before b. An immediate consequence is that the array elements bb[1], bb[2], bb[3], and bb[4] all exist. In other words the range of bb is bb[1..4]. We will refer to bb as a unit-offset vector.”

I stared at the book for near twenty minutes in utter disbelief.  It was the most beautiful thing I’d ever seen. How could I have never thought of that?  Countless hours of headaches I’d experienced in the past, gone!  Want to write a for loop from element 1 to 10, you write it from 1 to 10!  So nice.  I’ve been writing code for, I dunno, somewhere around 15 years and I’ve never thought of doing that.

You can get really fancy if you want to.  Using an array of pointers you can create a multi-dimensional matrix and access to elements from any range you like.

#define NR_END 1
double **dmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;

/* allocate pointers to rows */
m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror(“allocation failure 1 in matrix()”);
m += NR_END;
m -= nrl;

/* allocate rows and set pointers to them */
m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror(“allocation failure 2 in matrix()”);
m[nrl] += NR_END;
m[nrl] -= ncl;

for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

/* return pointer to array of pointers to rows */
return m;
}
Whoever wrote that routine, someone should kiss you. Actually, now that I think about it, if you had told the 15 year old Jason how to do that trick, I probably would have kissed you then too. I don’t think I’d have understood it at age 12 though. I had only started programming then.