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.