Not Honest Abe's hat...
un+1 = un + hf(t0+nh,un)
and u(t0) = u0 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 un in the vector field f with un+1 instead. Simple? Well, now you have un+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 = ut, giving us the system ut = v; vt = -∆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.
No comments:
Post a Comment