Today’s earthquakes near Sumatra fortunately didn’t lead to a major tsunami. But figuring out when tsunamis will develop is a difficult matter–and an interesting exercise in applied mathematics.

The main phenomenon is the propagation of so-called shallow water waves–water waves whose wavelength is large compared to the depth of the ocean. Those waves satisfy a partial differential equation (PDE) that was figured out in the 1800s. The equation is a nasty nonlinear one–that can’t be solved exactly.

I’ve been working on the numerical differential equation capabilities of *Mathematica* for more than a decade. Our goal is to automate the solutions of all types of equations–so users just have to enter their equation, and *Mathematica* then does all the analysis to select and apply the best algorithm.

The shallow-water equations are a good test–that I’m happy to say *Mathematica* passes with flying colors. One essentially just has to type the equations in, and get the solution, which is then easy to visualize–especially using the new visualization capabilities of *Mathematica* 6. (Click the image below to see the *Mathematica* animation.)

Let me explain a little about what’s involved in getting this.

Here’s an image of a *Mathematica* notebook that produces the main picture (you can download the notebook here):

The first thing is just the equations, given in *Mathematica* `StandardForm` (one could use `TraditionalForm` so it looks exactly like a traditional textbook, but it’s slightly harder to understand that way). Then we just use `NDSolve` to solve the equations. Then we use `Plot3D` to create a 3D visualization–complete with specular surfaces and everything. And finally, we use `Animate` to create an animation. (We can immediately export for the web using `Export`.)

What’s going on inside `NDSolve`? That one function is doing some pretty complicated things–it’s almost a microcosm of a century or two of applied mathematics.

`NDSolve` is doing a lot of things that can really only be done in *Mathematica*–by combining *Mathematica*‘s strength not only in numerical computation, but also in symbolic computation and discrete mathematics.

That you can just type an equation into `NDSolve` relies on *Mathematica*‘s general symbolic architecture. Once `NDSolve` gets an equation, it immediately determines the structure. In this case, it finds out that it’s been given a (2+1)-dimensional initial-value PDE.

For that type of equation, it forms a discrete grid in space (with the grid structure determined automatically to meet a certain accuracy criterion), then generates a large system of first-order ODEs on that grid, automatically incorporating all the necessary boundary conditions.

`NDSolve` has about 20 different families of methods for solving ODEs. In this case, it actually switches automatically between explicit and implicit methods depending on local stiffness. (The implicit methods get to use some of our very fast sparse matrix solving capabilities.)

In the end, `NDSolve` packages its solution as an `InterpolatingFunction`–a very convenient *Mathematica* construct that directly represents an approximate function, that in this case we chose to plot using `Plot3D`.

(Even though `NDSolve` can figure out what to do automatically, you can always give it hints, or even tell it exactly what method to use. In this case, you can improve the quality of the solution by choosing to use a pseudospectral method with a specified number of grid points.)

As one of the people responsible for `NDSolve`, I know how complicated all the things it does inside are. But when you use `NDSolve`, all you have to do is type your equation in, and let *Mathematica* do the rest. We’ve done the work (and it’s been a lot of work, I might add) to have everything run efficiently and automatically from there.

It might take a while to work out the physics of the shallow water wave approximation to a tsunami. But I think I can say that *Mathematica*‘s part of solving the equations could be accomplished before a tsunami has propagated very far.