Lecture - Implicit MethodsPatrick J. QuirkFebruary 10, 2005Lecture Outline:• Motivation for Implicit Methods: Stiff ODE’s– Stiff ODE Example: y0= −1000y∗ Clearly an analytical solution to this is y = e−1000t. This large negative factor in theexponent is a sign of a stiff ODE. It means this term will drop to zero and becomeinsignficant very quickly. Recalling how Forward Euler’s Method works, solving fory1gives:y1= y0+ h(−1000)y0= y0(1 + h(−1000))Now if |1+h(−1000)| > 1, then |y1| > |y0| and the series will diverge (approximationscontinue to get larger). Thus when solved with Forward Euler, this ODE has astability region of |1 + h(−1000)| < 1. Consequentially as the exponent grows, thetimestep must shrink to give acceptable results. In the worst case, you would haveto make your timestep so small that simulations would appear to be stopped.• Backward Euler’s Method– For a given differential equation system:ddtX(t) = f(X(t))∗ Forward Euler: Xn+1= Xn+ hf(Xn)· Evaluates f at the point we came from.∗ Backward Euler: Xn+1= Xn+ hf(Xn+1)· Evaluates f at the point we’re aiming at.– In general we can’t directly solve for Xn+1unless f happens to be a linear function.To deal with this we replace f(Xn+1) with a linear approximation based on the TaylorExpansion of f.Define ∆X = Xn+1− Xnand rewrite the above Backward Euler Method (BEM):Xn+ ∆X = Xn+ hf(Xn+ ∆X) ⇒ ∆X = hf(Xn+ ∆X)Now approximate f(Xn+ ∆X) by:f(Xn) + f0(Xn)∆X(Note that since f (Xn) is a vector, f0(Xn)). Using this approximation, we can approxi-mate ∆X with:∆X = h(f(Xn) + f0(Xn)∆X) ⇒ ∆X − hf0(X0)∆X = hf(X0)1Multiply through by 1/h and rearrange to get:1hI − f0(X0)∆X = f(X0)where I is the identity matrix. We can then solve for ∆X as:∆X =1hI − f0(X0)−1f(X0)Which allows us to compute:Xn+1= Xn+ ∆X (1)– Clearly this is much more work than using an explicit method. However this is notnecessarily a bad thing for several reasons:∗ The matrix f0will generally be sparse, and there are methods for solving sparsematrix systems quickly. In general equation (1) can be solved in linear time (thedimension of X).∗ Implicit methods allow for much bigger timesteps without losing stability.– As a result, explicit methods are more efficient than implicit methods because theirlarger time step often makes up for the difference in solving a linear system.– Big Picture: Because the code for an implicit solver is much more complicated than thatfor an explicit solver, make your problem un-stiff if possible. If you cannot, c ode up theimplicit method and take advantage of the larger step sizes.• Second-Order Equations– Many dynamics problems are written in terms of second-order ODE’s:¨x(t) = f(x(t),˙x(t))This can be coverted to a first-order ODE by introducing new variables. If we definev = ˙x:ddtx(t)v(t)=v(t)f(x(t), v(t))we then get a first-order system. However, applying BEM to it results in a linear systemsof size 2n × 2n where n is the dimension of x. Through some more substitutions, wecan rewrite the above equation to be an n × n system. Let ∆x = x(t0+ h) − x(t0) and∆v = v(t0+ h) − v(t0), then after applying BEM to the above equation:∆x∆v= hv0+ ∆vf(x0+ ∆x, v0+ ∆v)Applying a Taylor Series expansion to f (remeber is is 2 dimensions, x and v) gives thefirst-order approximation (evaluated at (x0, v0)):f(x0+ ∆x, y0+ ∆y) = f0+∂f∂x∆x +∂f∂v∆vSubstitute this back into the equation:∆x∆v= hv0+ ∆vf0+∂f∂x∆x +∂f∂v∆v2Substitue the first row into the second:∆v = hf0+∂f∂xh(v0+ ∆v) +∂f∂v∆vRegroup terms (I is the identity matrix):I − h∂f∂v− h2∂f∂x∆v = hf0+ h∂f∂xv0Then solve for ∆v and substitute that back into ∆x = h(v0+ ∆v) to get ∆x.This of course assume there is no variance in time. If f describes time-varying externalforces or references moving points or coordinate froms that are not variables of x, thenwe need to add an additional term to the above equation:I − h∂f∂v− h2∂f∂x∆v = hf0+ h∂f∂xv0+∂f∂t3Sources:1. SIGGRAPH Course Notes - http://www.pixar.com/companyinfo/research/pbm2001/notesd.pdf2. Mathworld Notes - http://mathworld.wolfram.com/3. Math 191 Lecture Notes - http://www.amath.unc.edu/Faculty/huang/teaching/math191/Notes/4. Desbrun, Mathieu. ”Interactive Animation of Structured Deformable Objec ts.” 1999 AvailableOnline: