Notes Monotonicity Test equation with real, negative Assignment 1 - PDF document

Notes Monotonicity � Test equation with real, negative � � Assignment 1 is not out yet :-) • True solution is x(t)=x 0 e � t , which smoothly decays to � http://www.cs.ubc.ca/~rbridson/ zero, doesn’t change sign ( monotone ) courses/533d-winter-2005 � Forward Euler at stability limit: • x=x 0 , -x 0 , x 0 , -x 0 , … � Not smooth, oscillating sign: garbage! � So monotonicity limit stricter than stability � RK3 has the same problem • But the even order RK are fine for linear problems • TVD-RK3 designed so that it’s fine when F.E. is, even for nonlinear problems! cs533d-winter-2005 1 cs533d-winter-2005 2 Monotonicity and Summary 1 Implicit Methods � Particle Systems: useful for lots of stuff � Backward Euler is unconditionally monotone � Need to move particles in velocity field • No problems with oscillation, just too much � Forward Euler damping • Simple, first choice unless problem has oscillation/rotation � Trapezoidal Rule suffers though, because � Runge-Kutta if happy to obey stability limit of that half-step of F.E. • Modified Euler may be cheapest method • Beware: could get ugly oscillation instead of • RK4 general purpose workhorse smooth damping • TVD-RK3 for more robustness with • For nonlinear problems, quite possibly hit nonlinearity (more on this later in the course!) instability cs533d-winter-2005 3 cs533d-winter-2005 4 Summary 2 Second Order Motion � If stability limit is a problem, look at implicit methods • e.g. need to guarantee a frame-rate, or explicit time steps are way too small � Trapezoidal Rule • If monotonicity isn’t a problem � Backward Euler • Almost always works, but may over-damp! cs533d-winter-2005 5 cs533d-winter-2005 6

Second Order Motion What’s New? � If particle state is just position (and colour, size, …) then � If x =(x,v) this is just a special form of 1st 1st order motion order: d x /dt= v (x,t) • No inertia • Good for very light particles that stay suspended : smoke, dust… � But since we know the special structure, • Good for some special cases (hacks) can we take advantage of it? � But most often, want inertia (i.e. better time integration algorithms) • State includes velocity, specify acceleration • More stability for less cost? • Can then do parabolic arcs due to gravity, etc. • Handle position and velocity differently to � This puts us in the realm of standard Newtonian physics • F=ma better control error? � Alternatively put: • dx/dt=v • dv/dt=F(x,v,t)/m (i.e. a(x,v,t) ) � For systems (with many masses) say dv/dt=M -1 F(x,v,t) where M is the “mass matrix” - masses on the diagonal cs533d-winter-2005 7 cs533d-winter-2005 8 Linear Analysis Constant Acceleration � Solution is v ( t ) = v 0 + a 0 t � Approximate acceleration: 2 a 0 t 2 x ( t ) = x 0 + v 0 t + 1 ) � a 0 + � a � x x + � a a x , v ( � v v � No problem to get v(t) right: just need 1st order accuracy � Split up analysis into different cases � But x(t) demands 2nd order accuracy � So we can look at mixed methods: � Begin with first term dominating: • 1st order in v constant acceleration • 2nd order in x • e.g. gravity is most important cs533d-winter-2005 9 cs533d-winter-2005 10 Linear Acceleration More Approximations… � Dependence on x and v dominates: � Typically K and D are symmetric semi-definite (there are good reasons) a(x,v)=-Kx-Dv • What does this mean about their eigenvalues? � Do the analysis as before: � Often, D is a linear combination of K and I (“Rayleigh damping”), or at least close to it � � � � � � x 0 I x d • Then K and D have the same eigenvectors � = � � � � � (but different eigenvalues) dt v � K � D v � � � � � � • Then the eigenvectors of the Jacobian are of the form (u, � u) T • [work out what � is in terms of � K and � D ] � Eigenvalues of this matrix? cs533d-winter-2005 11 cs533d-winter-2005 12

Simplification Split Into More Cases � � is the eigenvalue of the Jacobian, and � Still messy! Simplify further 2 � � K � If D dominates (e.g. air drag, damping) � = � 1 ( 1 ) 2 � D ± 2 � D � � � � D , 0 { } � Same as eigenvalues of � 0 1 � � � • Exponential decay and constant � � K � � D � � � If K dominates (e.g. spring force) � Can replace K and D (matrices) with corresponding eigenvalues (scalars) • Just have to analyze 2x2 system � � ± i � K cs533d-winter-2005 13 cs533d-winter-2005 14 Three Test Equations Explicit methods from before � Constant acceleration (e.g. gravity) � Forward Euler • a(x,v,t)=g • Constant acceleration: bad (1st order) • Want exact (2nd order accurate) position • Position dependence: very bad (unstable) � Position dependence (e.g. spring force) • Velocity dependence: ok (conditionally • a(x,v,t)=-Kx monotone/stable) • Want stability but low damping � RK3 and RK4 • Look at imaginary axis • Constant acceleration: great (high order) � Velocity dependence (e.g. damping) • Position dependence: ok (conditionally stable, but • a(x,v,t)=-Dv damps out oscillation) • Want stability, smooth decay • Velocity dependence: ok (conditionally • Look at negative real axis monotone/stable) cs533d-winter-2005 15 cs533d-winter-2005 16 Implicit methods from before Setting Up Implicit Solves � Backward Euler � Let’s take a look at actually using Backwards Euler, for example • Constant acceleration: bad (1st order) x n + 1 = x n + � t v n + 1 • Position dependence: ok (stable, but damps) v n + 1 = v n + � t M � 1 F x n + 1 , v n + 1 ( ) • Velocity dependence: great (monotone) � Eliminate position, solve for velocity: � Trapezoidal Rule v n + 1 = v n + � t M � 1 F x n + � t v n + 1 , v n + 1 ( ) • Constant acceleration: great (2nd order) � Linearize at guess v k , solving for v n+1 � v k + � v • Position dependence: great (stable, no � � ) + � t � F � x � v + � F v k + � v = v n + � t M � 1 F x n + � tv k , v k damping) ( � v � v � � • Velocity dependence: good (stable but only � � � Collect terms, multiply by M conditionally monotone) � M � � t � F � v � � t 2 � F � � v = M v n � v k ( ) + � t F x n + � tv k , v k ( ) � � � � x � cs533d-winter-2005 17 cs533d-winter-2005 18

Symmetry Specialized 2nd Order Methods � Why multiply by M? � This is again a big subject � F position � F � Physics often demands that and velocity � Again look at explicit methods, implicit � x are symmetric � v methods • And M is symmetric, so this means matrix is symmetric, hence easier to solve � Also can treat position and velocity • (Actually, physics generally says matrix is SPD - even dependence differently: better) • If the masses are not equal, the acceleration form of mixed implicit-explicit methods the equations results in an unsymmetric matrix - bad. � F � Unfortunately the matrix is usually velocity unsymmetric � x • Makes solving with it considerably less efficient • See Baraff & Witkin, “Large steps in cloth simulation”, SIGGRAPH ‘98 for one solution: throw out bad part cs533d-winter-2005 19 cs533d-winter-2005 20 Symplectic Euler Symplectic Euler performance � Like Forward Euler, but updated velocity used � Constant acceleration: bad for position • Velocity right, position off by O( � t) ( ) v n + 1 = v n + � ta x n , v n � Position dependence: good � t < 2 • Stability limit x n + 1 = x n + � tv n + 1 K � Some people flip the steps (= relabel v n ) • No damping! � (Symplectic means certain qualities are � Velocity dependence: ok preserved in discretization; useful in science, not • Monotone limit necessarily in graphics) � t < 1 D • Stability limit � [work out test cases] � t < 2 D cs533d-winter-2005 21 cs533d-winter-2005 22 Tweaking Symplectic Euler Staggered Symplectic Euler � [sketch algorithms] � Constant acceleration: great! • Position is exact now � Stagger the velocity to improve x � Other cases not effected � Start off with v 12 = v 0 + 1 2 � ta x 0 , v 0 ( ) • Was that magic? Main part of algorithm unchanged (apart from relabeling) yet now it’s more accurate! � Then proceed with � Only downside to staggering ( ) v n + 12 = v n � 12 + 1 2 ( t n + 1 � t n � 1 ) a x n , v n � 12 • At intermediate times, position and velocity not known x n + 1 = x n + � tv n + 12 together • May need to think a bit more about collisions and � Finish off with other interactions with outside algorithms… ( ) v N = v N � 12 + 1 2 � ta x N , v N � 12 cs533d-winter-2005 23 cs533d-winter-2005 24