Transcript Implicit Methods

Implicit Methods
Patrick Quirk
February 10, 2005
Introduction, Motivation
• Explicit methods have a difficult time with some
types of ODE’s
• They guess where the function is heading
• Stiff ODE’s
– One component big for a short time
– Does not contribute later
– Must make time step small for explicit to handle it
y '  1000 y  yn 1  yn (1  1000h)
| 1  1000h | 1  0  h  1 / 500
Solution – Implicit Methods
• Backward Euler Method
• Evaluates ƒ at the point we’re aiming at
– Less guessing, more accuracy
• Consequently, the step size can be larger
y '  1000 y  yn 1  yn /(1  1000h)
| 1  1000h | 1  h  0
Implicit Methods
• Last example was simplisitc
• Typically cannot directly solve for yn+1
– Unless ƒ is linear
• So…approximate ƒ with its Taylor
Expansion
f ( xn  x)  f ( xn )  f ' ( xn )x
Implicit Methods
• More equation manipulations are done
• Allows us to approximate Δx
– Albeit with an inverse of a matrix!
• In general this extra computation isn’t bad
– Larger time step offsets this
An Aside - Explicit Method
• Verlet Algorithm
• Does away with the velocity term
– Stores only current, previous position
– Approximates velocity by xn-xn-1
xn 1  2 xn  xn 1  an  dt
2
An Aside - Explicit Method
• But it’s explicit!
– It’s a second-order method, higher accuracy
– Has the computational cost of a first-order
– No numerical drift like in some other explicit
and implicit methods
– Easy to code!
– Widely used in games and physics
– Often chosen over implicit methods for these
reasons
Big Picture - Implicit Methods
• Explicit Method code is easy
• Implicit Method code is harder
• Try to make problem un-stiff if possible
– Then use an explicit method and small Δt
• Otherwise write an implicit solver
– Use the larger time step to your advantage
Second Order ODEs
• Dynamics problems are often second order
x(t )  f (v(t ), a(t ))
• Need to convert it to a first order problem
Second Order ODEs
• Introduce new variables to get first order
v(t )

d  x(t )  

  

dt  v(t )   f ( x(t ), v(t )) 
• But this gives a 2nx2n linear system
– Bad news for BEM, has to solve that system
– Can we reduce the size?
Second Order ODEs
• Can reduce it to nxn
– Introduce some more new variables
x  x(t0  h)  x(t0 )
v  v(t0  h)  v(t0 )
• Taylor Expand ƒ (in two dimensions)
– Put that back into the equation
Second Order ODEs
• Substitute the first row into the second
• Regroup terms and solve for Δv
• If there is a variance in time, an extra term:
f
f
f 


2 f 
I h h
v  h f 0  h v0  
v
x 
x
t 


Sources
•
SIGGRAPH Course Notes
–
•
Mathworld References
–
•
http://mathworld.wolfram.com
“A Practical Dynamics System”. Kačić-Alesić, Nordenstam, Bullock. Eurographics/SIGGRAPH
2003.
–
•
http://www.pixar.com/companyinfo/research/pbm2001/notesd.pdf
http://portal.acm.org/citation.cfm?id=846276.846278
“Advanced Character Physics”. Thomas Jakobsen. Gamasutra. 2003
–
http://www.gamasutra.com/resource_guide/20030121/jacobson_01.shtml
Interactive Animation of
Structured Deformable Objects
Mathieu Desbrun, Peter Schröder, Alan Barr
Caltech
Introduction, Motivation
• Interactive deformation is tough
– Simulation methods are too slow
– Usually used in VR, user has control
• Limited greatly by time step
– Too big means instability in simulation
– Too small means alteration in reality
• Implicit Integration
– Offers low computational time
– Nearly arbitrary time step
Implicit Integration : 1D case
n1
vi
n1
xi
n
n+1
n dt
 v i  Fi
n
n1
 xi  vi
m
dt
• Replace all forces at
time t by t+1
• Now positions are not
blindly reached, but are
logically inferred.
In other words:
• Explicit steps into the unknown with only initial
conditions.
• Implicit tries to hit the next position correctly.
Implicit Integration : 1D case
• Must compute Fn+1 at t+dt (spring forces)
• Use a first order approximation (since we don’t
know what exactly Fn+1(t+dt) is)
F
n1
 F
n
F n1

 x
x
• Here Δn+1x is called the backward difference
operator:
n1

x  x
n1
x
n
Implicit Integration : 1D case
• Add artificial viscosity
– Increases stability
• Force filtering
– Uses matrix multiplication (H=dF/dx)
– Global force effects
– Smoothes large force differences
• Pull it all together…
Implicit Integration : 1D case
• Final Form for 1D:


dt
 v   I 
H 
m 

n 1
2
1


dt
F  dtHv
m
n
n
Damping Force
Force Filter Matrix
• Look closely, is it really implicit?
Implicit Integration : 1D case
• Loss in accuracy
– Force smoothing
– Artificial viscosity
• Force filtering matrix
changes at each time
step (2D/3D only)
– Forces you to solve a
linear system
• Gain in speed
• Gain in speed
• Gain in speed
Implicit Integration: 1D case
• Some of the “Bads” aren’t all that bad:
• Artificial viscosity
– It is frequently added to simulations.
– If it is already taken care of, so much the better.
• Changing filter matrix
– Even if this does take time, the computational
advantages are already so significant it doesn’t
matter too much.
– Approximating it makes it a constant
• Compute it once, use it forever
Implicit Integration: 2D/3D case
• Similar to 1D case, though the smoothing matrix
changes with each step.
– This paper removes the need to solve a linear system
• Splits the force into a linear and non-linear part
– Linear is easily solved, non-linear approximated to
zero
• Magnitude does not change (internal forces), just direction
• Preserve linear/angular momenta
• Includes inverse dynamics to remove stretching
Implicit Integration : 2D/3D case
F
H
x
• In 2D/3D, implicit integration is done by
computing the Hessian matrix H
– It depends on position of particle nodes, changing
• Equation for one element of H is too long to
show here
– Each node is a 3x3 matrix, total of 3nx3n elements
• Too large to solve efficiently
– Approximate it
Implicit Integration : 2D/3D case
• Matrix split into 2 parts
– Linear
 H ij  kij if i  j
H  
 j i kij
 ii
– Non-Linear
• Simple, assume its zero
• It’s constant in magnitude, but will rotate
Implicit Integration : 2D/3D case
• Check if linear, angular momenta are conserved
– Linear is
• Sum of artificial viscosity forces is zero
• Force filtering has no effect
– Angular is not!
• Because of that assumption about the non-linear hessian
• Desbrun corrects this with some inverse dynamics
Results
Sources
•
“Interactive Animation of Structured Deformable Objects”. Desbrun, Schröder, Barr. Graphics
Interface. 1999
–
http://www.multires.caltech.edu/pubs/GI99.pdf