Physically Based ModelingImplicit Methods for Differential EquationsDavid BaraffPixar Animation StudiosPlease note: This document is 2001 by David Baraff. This chapter may be freelyduplicated and distributed so long as no consideration is received in return, andthis copyright notice remains intact.Implicit Methods for Differential EquationsDavid BaraffPixar Animation Studios1 Implicit MethodsThe methods we have looked at for solving differential equations in the first section of these notes(Euler’s method, the midpoint method) are all called “explicit” methods for solving ODE’s. How-ever, sometimes an ODE can become “stiff,” in which case explicit methods don’t do a very goodjob of solving them. Whenever possible, it is desirable to change your problem formulation so thatyou don’t have to solve a stiff ODE. Sometimes however that’s not possible, and you have just haveto be able to solve stiff ODE’s. If that’s the case, you’ll usually have to use an ODE solution methodwhich is “implicit.”2 Example Stiff ODEFirst, what is the meaning and cause of stiff equations? Lets consider an example that arises fre-quently in dynamics. Suppose that we have a particle, with position (x(t), y(t)), and suppose thatwe want the y-coordinate to always be zero. One way of doing this is to add a component −ky(t)to ˙y(t) where k is a large positive constant. If k is large enough, then the particle will never movetoo far away from y(t) = 0, since the −ky(t) term always brings y(t) back towards zero. However,lets assume that there is no restriction on the x-coordinate, and that we want a user to be able to pullthe particle arbitrarily along the x-axis. So lets assume that over some time interval our differentialequation is simply˙X(t) =ddtx(t)y(t)=−x(t)−ky(t). (2–1)(We’ll also assume that the particle doesn’t start exactly with y0= 0.) What’s happening here is thatthe particle is strongly attracted to the line y = 0, and less strongly towards x = 0. If we solve theODE far enough forward in time, we expect the particle’s location to converge towards (0, 0) andthen stay there once it arrives.Now suppose we use Euler’s method to solve the equation. If we take a step of size h,wegetXnew= X0+ h˙X(t0) =x0y0+ h−x0−ky0.This yieldsXnew=x0− hx0y0− hky0=(1 − h)x0(1 − hk)y0.D1If we look at the y component of this equation, we see that if |1 − hk| > 1 then the ynewwe computewill have an absolute value which is larger than |y0|. In other words, if |1 − hk| > 1, Euler’s methodwill not converge to an answer: each step will result in a value of ynewwhich is larger than the last.Technically, Euler’s method isunstable for |1− hk| > 1. Thus, webetter have 1− hk > −1orhk < 2if we hope to converge. The largest step we can hope to take is less than 2/ k.Now, if k is a large number, we’ll have to take very small steps. This means that the particleslides towards (0, 0) excruciatingly slowly. Even though the particle may nearly satisfy y0= 0, wehave to take such small steps that the particles’ progress along the x-axis is pretty much nonexistent.That’s the embodiment of a stiff ODE. In this case, the stiffness arises from making k very large inorder to keep the particle close to the line y = 0. Later on, when we connect particles with second-order dynamics together with springs, we’ll experience exactly the same effect: stiff ODE’s. Evenif we use a more sophisticated explicit method such as fourth-order Runge-Kutta, we may do a littlebetter in the size of our steps, but we’ll still have major problems.Now as we said above, the name of the game is to pose your dynamics problems so that you don’texperience stiff ODE’s. However, when you can’t, you’ll have to turn towards an implicit solutionmethod. The method we’ll show below is the simplest of the implicit methods, and its based ontaking an Euler step “backwards.”3 Solving Stiff ODE’sGiven a differential equationddtX(t) = f(X(t)),the explicit Euler update would be Xnew= X0+ hf(X(t0)), to advance the system forward h in time.For a stiff problem though, we change the update to instead beXnew= X0+ hf(Xnew) (3–1)That is, we’re going to evaluate f at the point we’re aiming at, rather than where we came from. (Ifyou think about reversing the world and running everything backwards, the above equation makesperfect sense. Then the equation says “if you were at Xnew, and took a step −hf(Xnew), you’d end upat X0.” So if your differential equation represents a system that is reversible in time, this step makessense. It’s just finding a point Xnewsuch that if you ran time backwards, you’d end up at X0.) So,we’re looking for a Xnewsuch that f, evaluated there, times h, points directly back at where we camefrom. Unfortunately, we can’t in general solve for Xnew, unless f happens to be a linear function.To cope with this, we’ll replace f(Xnew) with a linear approximation, again based on f ’s Taylorseries. Lets define 1X by 1X = Xnew− X0. Using this, we rewrite equation (3–1) asX0+ 1X = X0+ hf(X0+ 1X).or just1X = hf(X0+ 1X).Next, lets approximate f(X0+ 1X) byf(X0) + f0(X0)1X.SIGGRAPH 2001 COURSE NOTES D2 PHYSICALLY BASED MODELING(Note that since f(X0) is a vector, the derivative f0(X0) is a matrix.) Using this approximation, wecan approximate 1X with1X = hf(X0) + f0(X0)1X.or1X − hf0(X0)1X = hf(X0)Rewriting this as1hI − f0(X0)1X = f(X0),where I is the identity matrix, we can solve for 1X as1X =1hI − f0(X0)−1f(X0) (3–2)Computing Xnew= X0+ 1X is clearly more work than using an explicit method, since we haveto solve a linear system at each step. While this would seem a serious weakness, computationallyspeaking, don’t despair (yet). For many types of problems, the matrix f0will be sparse—for exam-ple, if we are simulating a spring-lattice, f0will have a structure which matches the connectivityof the particles. (For a discussion of sparsity and solution techniques, see Baraff and Witkin [1].Basic material in Press et al. [2] will also prove useful.) As a result, it is usually possible to solveequation (3–2) in linear time (i.e. time proportional to the dimension of X). In such cases, the payoffis dramatic: we can usually take considerably large timesteps without losing stability (i.e. withoutdivergence, as happens with the explicit case if the stepsize is too large). The time taken in solvingeach
View Full Document
Unlocking...