## Finite Elements Methods in Mathematica

In Mathematica 10 there’s now the ability to solve differential equations using finite element methods (FEM). I’ve spent some time this summer playing around with the tools and I thought I’d write this to help me remember some stuff.

The typical problem to be solved with FEM involves trying to find a function on a 2D surface that obeys a (typically) second order differential equation. A few examples include

• The voltage electric potential inside a region where the voltage potential on the walls is known
• The shape of a rubber sheet hung from a frame
• The possible standing waves on a surface
• In quantum mechanics this could be possible eigenfunctions
• In acoustics this represents the possible resonance modes

Normally you divide the surface into lots of triangles and then turn the continuous differential equation into a finite number of linear equations, solving for the value of the unknown function on the nodes. Here’s a great (though long) description of how it works.

For years I’ve been jealous of my friends who are facile with software that can do this sort of computation, because if you try any of those examples above with Mathematica 9 or below it laughs at you and gives up. But now, in Mathematica 10, it’s built into my favorite command for doing differential equations: NDSolve.

The first thing you have to do is define the surface you’re interested in. There are a ton of built in shapes like Disk[] that you can use, but you can also define them using commands like ImplicitRegion and RegionUnion and RegionDifference. You can also make regions out of graphics and a lot of the built-in data in Mathematica like the shape of countries, states, counties, and cities.

Once you have the region specified, you can do a command like this:

NDSolve[{Laplacian[u[x,y], {x,y}]==0, Boundary conditions ...}, u, {x,y}\Element region]

where I’ll discuss the boundary conditions below (and “region” has already been defined). The Laplacian command is equivalent to:

$\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}$

Alright, here are some examples:

Solves the hanging rubber sheet problem with constantly changing frame.

For this example I set the Lagrangian equal to a nonzero value and set the frame (borders) to be interesting shapes (every frame in the gif is a different solution). You set the border using a Dirichlet boundary condition like this: DirichletCondition[u[x,y]==some cool function, border you’re interested in]. In this example I did a different condition for the inner and outer borders.

Treating Minnesota like a drum head and thumping it.

Here I solved the following equation on the region in the shape of Minnesota (obtained using DiscretizeGraphics[EntityValue[ctrl-=Minnesota, “Polygon”]])

$\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}==\frac{\partial^2 u}{\partial t^2}$

Another way to do the same thing is to find the resonant modes of Minnesota and combining them to start as a thump:

Similar to the last one but using the eigenmodes of Minnesota to approximate a “thump”

By playing with Neumann boundary conditions (where you can set the value of the slope at the boundary) I tried to simulate flapping flags:

The city of St. Paul, MN (colored with the MN flag) “waving”

I didn’t like the edges very much, so I tried again with a square US flag, adjusting the boundary conditions a little more:

US flag flapping

Note that you can get the flag image using flag=EntityValue[ctrl=United States, “Flag”] and you can put it on 3D graphics using PlotStyle->Texture[flag].

Finally I’ve played quite a bit with various drum heads and the sounds they’d make. You solve this equation:

$\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}==\omega^2 u$

But you do it with an awesome brand new command in 10.2: NDEigensystem. It takes just the left hand side of the equation (with an odd negative sign that took me a while to get), along with appropriate boundary conditions. It returns the values of $\omega^2$ and interpolation functions for the resonant functions. Here’s the 25th resonance of a Texas-shaped drum:

25th resonance of a Texas-shaped drum head

You can use the $\omega$ ‘s to figure out what a drum head will sound like, and you can figure out how much of each resonance you’ll hear based on where you thump the drum by expanding your “thump” function into the resonant functions. I cheat and use a delta function thump so I just need to know the value of the resonant function at the location that I want the thump. By the way, the functions returned by NDEigensystem are orthonormal, which makes my last sentence somewhat more palatable to my mathematician friends.

Here’s a movie I made combining Mathematica’s ability to do all the calculations, image generation, Minnesota data look up, Traveling Salesman solutions, and sound generation with FFMpeg’s ability to combine images and sound into a movie. It shows what it would sound like if you thumped a Minnesota-shaped drum head in the center of all 89 counties in an order determined by the shortest tour through those points.

Fun times, huh? I’ve certainly enjoyed learning a ton about FEM and what you can do with it, and I’ve gotten decently facile with Mathematica’s implementation.

Thoughts? Here are some starters for you:

• Why did you strike out voltage in the early examples?
• You said it was typically second order problems. What about the beam equation, huh? That’s fourth order, can you handle that? I bet you’d have to use something like COMSOL to do that, not wimpy Mathematica (you’re right).
• Hey, I’ve done some of this stuff with Mathematica (x<10), you’re a liar.
• What’s up with the negative sign you have to use in NDEigensystem?
• Your description of a delta function thump has tons of problems, not least of which is that you have to calculate all the resonances, and I have a sneaky suspicion you only do, oh, I don’t know, 20.
• How did you decide to put dampening in the sound? It’s clear you’re not just playing those frequencies combined.
• Why didn’t you do this in python, that seems the obvious choice for anything computer related?
• How could you utilize this when teaching Laplacian solutions in a electricity and magnetism class?
• I tried your flapping flag simulation and I got some negative eigenvalues. What the heck? (Yeah, I get that too when I set the Neumann boundary condition equal to the function itself, which is what I thought a flapping flag would look like. My mathematician friend tells me that then I can’t trust any of this.)
• Why do you only do Minnesota examples? What, you think you’re too good for the rest of the country?

Professor of physics at Hamline University in St. Paul, MN
This entry was posted in mathematica. Bookmark the permalink.

### 3 Responses to Finite Elements Methods in Mathematica

1. Zoe Boekelheide says:

Very cool. I have canned questions (1) and (8) for you! I am teaching upper level E&M in the fall. The Laplacian solutions are 3D, not 2D, so harder to visualize, but I guess you could just pick a problem with only two interesting dimensions? Or solve in 3D but visualize only 2 dimensions at a time?

• Andy "SuperFly" Rundquist says:

I struck out voltage because it seems that you should be able to use a volt meter to measure voltage and it doesn’t work in an empty space. What do you think?

Yeah, it works in 3D as well, you just have to choose how to visualize it.