Origami Dodecahedron

Physically made this one. 60 folded pieces. Took a couple of months.

Caustic Echoes

Greetings. I know it's been a long while since I last posted. Sorry about that. But rest assured, this has not been abandoned! I come blogging from MSRI in Berkeley, in a conference on General Relativity. Here is the coolest visualization I have seen in the conference today: caustic echoes of waves in a Schwarzschild spacetime (the standard example of a non-rotating black hole and the surrounding space).

A non-cubical example!

Ahh, finally, an example with the appropriate symmetry... Learned how to create some interesting 3D meshes. This example is essentially the same as the previous, but with a different global mesh. I also changed the numerical method to something symplectic, and of higher order: Störmer-Verlet timestepping. The energy that is preserved here (or, only close to being preserved, actually—symplectic integrators do not exactly preserve energy, but do a lot better than more generic methods) is the field energy, 1/2(E, D) + 1/2 (B, H). Or, the energy that would be (close to) preserved, if it weren't for this oscillation being a driven one, caused by currents in the wire, rather than natural oscillations given by some initial, non-equilibrium configuration (as in the vibrating drum post). It's time for a different example, perhaps the analogue of a cavity resonator.

More Vector Finite Element Goodness

Would you like corn on the cob with that?

In case you weren't satisfied with steady-state behavior in the past couple of posts, I finally bring you some dynamics. This is a picture of the magnetic field caused by a central wire with a given alternating current, surrounded by air (although with boundary conditions tangent are to the surrounding cube, so the vector field becomes a bit squashed toward the outside as compared to an actually cylindrically symmetric problem). Current in the wire is modeled by simply taking a current density which is nonzero (and time-varying) within a small radius from the z-axis. The mesh is significantly more refined in the region, to better account for the discontinuity and that I intend the wire to actually be round. It takes a little while for it to reach steady state in terms of amplitude, so don't be disappointed if the vectors aren't very visible initially. If you look closely (and admittedly, you will have to either really look, or play the video in slower motion), you can see that there's a slight delay between what happens at the wire, and what happens some distance away from it. Here c = 3 (so today, E = 9m), that is, it takes 1/3 of a unit of time (about 2 seconds of the video) for the influence of the wave to traverse 1 cube length. The whole video itself is 4 units of time long. The original equations as presented in the paper of Bossavit is in fact more geared toward static situations, in which we can neglect the time derivative term for the electric displacement (the one that Maxwell added), but it turns out adding that term back makes the finite element formulation much more natural—in fact, when reducing it to mass and stiffness matrices and nodal vectors, the time-stepper we use is exactly the same as in the backward Euler discretization described a couple of posts ago. Without that term, I had to use  two different mass matrices. Here's a more interesting example:

And a donut too!

Here, we've modified the setup to include a conductor in the shape of, you guessed it, a torus. The significance of it being "conductive" is that an extra term, Ohm's law, is added in. It affects the dynamics because the field causes currents to flow inside the torus, which generates more field, in different directions. This is expressed in the program as a piecewise function which is zero outside the toroidal region and 1 inside it, and the stiffness matrix is modified to account for this extra term. For similar reasons as the wire above, the mesh is now also significantly more refined in the region of the torus (otherwise it kind of looks like a square block), and in any case, also better for handling the discontinuity, but this time for conductivity. Next up: seeing if there's a way to go about this using symplectic methods!

Happy Thanksgiving from Nested Tori!

It looks kind of turkey-ish, right? Some numerical imperfections give it a green tail, anyway. It's supposed to be the magnetic field of a spinning solid ball of charge. In some sense, this problem started my quest to become more numerical: I remember being frustrated that I could not rattle off what the magnetic field of this configuration of charge was, despite knowing a whole bunch of theory and equations. And attacking the problem analytically was, to use a different species of fowl, a wild goose chase. It turns out I did find out one possible analytic expression for it, using multipole expansions, but that's still a far cry from being able to come up with interesting visualizations (multipole expansions and spherical harmonics DO have cool visualizations which we will definitely post here someday). Gobble gobble! The field is actually similar to that of a dipole, and aside from the boundary condition that it has to be tangent to the big enclosing cube, the vector field is tangent to, you guessed it, nested tori.

Sneak Peek: Vector Finite Elements

Vector-valued differential equations are often a topic felt to be too advanced to be mentioned in elementary differential equations courses. Part of the problem, of course, is that many of the equations actually involve things like differential forms and other vector-like quantities that are best explained in terms of differential geometry, which takes an introductory course too far afield. The following picture is the magnetic field of current flowing through a solid cube (it is a constant current coming approximately out of the plane of the picture) requiring it to be tangent to the boundary. The method of solution involves special kinds of finite elements specifically for vector-valued functions (really, differential forms), and as talked about in Anil Hirani's talk. I'll be working with these types of elements for a while, so expect to see more posts on this!

Hyperbolic Equations and Symplectic Integrators

One well-known analytic fact about the wave equation (hyperbolic equations) vs. the heat equation (parabolic equations) is that the latter always results in smooth solutions, for any time t > 0, for even very nonsmooth initial data; whereas the former has no smoothing properties at all. This is usually summed up as "Wave singularities propagate." Here is the video with the same (very nonsmooth) initial data that was given for the heat equation last time:

Not Honest Abe's hat...

Pretty ridiculously jagged, right? However, it turns out that in order to actually (reasonably) faithfully preserve the jaggedness, you have to tweak your numerical methods a bit to achieve this. (We'll give a few words on why one would want to preserve jaggedness at the end). The most basic method to discretize time evolution problems is simply to replace time derivatives with their difference quotient approximations—the derivative without a limit: u'(t) ≈ 1/∆t (u(t+∆t) - u(t)). This enables us to solve for the future time in terms of the past: u(t+∆t) ≈ u(t) + u'(t)∆t. Now the essence of differential equations is that the u' is also given by some other function (vector field) f(t,u). Fixing the ∆t to some standard "mesh size" h, and having the initial value, this gives us an iteration process (recurrence relation):

u_n+1 = u_n + hf(t₀+nh,u_n)

and u(t₀) = u₀ is the given initial condition. This is called Euler's Method. It is well-known, however, that this method is sensitive to the size of the time step (it is unstable). It turns out that one can get a stable method by simply considering backward differences, that is, replacing the u_n in the vector field f with u_n+1 instead. Simple? Well, now you have u_n+1 on both sides of the equation, and for all but the simplest fs, it can't be easily solved algebraically (so we have to do some more numerical computations). For linear ODEs, like the heat and wave equations, this turns out to simply be inverting a different matrix.

For wave equations, with two time derivatives, one usually does a classic (but as far as I know, nameless) trick, setting v = u_t, giving us the system u_t = v; v_t = -∆u, then treating it as first-order. It's a bit more complex than that, because we also have to discretize the spatial part (with finite elements—more fun with FEniCS). The backward Euler discretization gives: 

Hey, why is it a lot smoother?! As explained to me by Prof. Melvin Leok in his class on geometric integration, it is essentially because stability theory is based on exponential decay (the model problem), and so in the goal to achieve stability, the methods actually introduce artificial dissipation (smoothing) terms. This means the solution loses energy. Energy conservation (which preserves the jaggedness) in the wave equation, however, is important for many applications. For example, if you simulate the solar system, you'll get some wildly off predictions, since energy loss there corresponds to orbits spiraling in. And even if we have a dissipative system, we probably don't want the wrong amount (however, there's one case I can think of in which it might not matter: solving an elliptic problem by running a parabolic problem to equilibrium).

The fix in this particular case (Euler-type discretizations) is interesting and almost seems like a hack at first (but really, there's lots of good theory backing this!) is that we use forward Euler in one variable (say, u) and the backward Euler in the other variable (v). Either way works (one is called Euler-A, the other Euler-B). This is called a symplectic integrator. In some sense, the dissipative behavior of the backward method, combined with the instability (where things exponentially grow) of the forward way, cancel out. Not completely, as it does have timestep sensitivity (too large a timestep causes our solution to simply blow up, even if is smooth like last time). Things can be written down a lot more explicitly with simpler, 1-particle systems such as a pendulum, or simulation of one planet in the solar system, but that is a story for another entry.

For fun, I leave you one last video, which is the symplectic case, but now with only the color-coding and not the actual displacement in 3D, which might be more informative. Note how you can really see the waves "reflecting" at the boundary—that's the Dirichlet boundary condition.

This is your brain. On math.

And here's a fun still: