Notes Trapezoidal Rule Again � Most of assignment 1 hasn � t been covered � The method: in class yet, but after today you should be ( ) x n + 1 = x n + � t 2 v ( x n , t n ) + 1 1 2 v ( x n + 1 , t n + 1 ) able to do a lot of it � Forgot to include instructions about � Let � s work out stability: view_obj: ( ) x n + 1 = x n + � t 2 � x n + 1 2 � x n + 1 1 • To navigate, hold down shift and click/drag ( ) x n + 1 = 1 + 1 ( ) x n with left, right, or middle mouse buttons (same 1 � 1 2 � � t 2 � � t navigation model as Maya) x n + 1 = 1 + 1 2 � � 2 � � t x n 1 � 1 cs533d-term1-2005 1 cs533d-term1-2005 2 Monotonicity and Monotonicity Implicit Methods � Test equation with real, negative � � Backward Euler is unconditionally • True solution is x(t)=x 0 e � t , which smoothly decays to monotone zero, doesn � t change sign ( monotone ) • No problems with oscillation, just too much � Forward Euler at stability limit: damping • x=x 0 , -x 0 , x 0 , -x 0 , … � Trapezoidal Rule suffers though, because � Not smooth, oscillating sign: garbage! of that half-step of F.E. � So monotonicity limit stricter than stability in this case • Beware: could get ugly oscillation instead of � RK3 has the same problem smooth damping • 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-term1-2005 3 cs533d-term1-2005 4
Summary 1 Summary 2 � Particle Systems: useful for lots of stuff � If stability limit is a problem, look at implicit methods � Need to move particles in velocity field • e.g. need to guarantee a frame-rate, or � Forward Euler • Simple, first choice unless problem has explicit time steps are way too small oscillation/rotation � Trapezoidal Rule � Runge-Kutta if happy to obey stability limit • If monotonicity isn � t a problem • Modified Euler may be cheapest method � Backward Euler • RK4 general purpose workhorse • Almost always works, but may over-damp! • TVD-RK3 for more robustness with nonlinearity (more on this later in the course!) cs533d-term1-2005 5 cs533d-term1-2005 6 Second Order Motion Second Order Motion � If particle state is just position (and colour, size, …) then 1st order motion • No inertia • Good for very light particles that stay suspended: smoke, dust… • Good for some special cases (hacks) � But most often, want inertia • State includes velocity, specify acceleration • Can then do parabolic arcs due to gravity, etc. � This puts us in the realm of standard Newtonian physics • F=ma � 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-term1-2005 7 cs533d-term1-2005 8
What’s New? Linear Analysis � If x =(x,v) this is just a special form of 1st � Approximate acceleration: order: d x /dt= v (x,t) ) � a 0 + � a � x x + � a � But since we know the special structure, ( a x , v � v v can we take advantage of it? (i.e. better time integration algorithms) • More stability for less cost? � Split up analysis into different cases • Handle position and velocity differently to � Begin with first term dominating: better control error? constant acceleration • e.g. gravity is most important cs533d-term1-2005 9 cs533d-term1-2005 10 Constant Acceleration Linear Acceleration � Solution is � Dependence on x and v dominates: v ( t ) = v 0 + a 0 t a(x,v)=-Kx-Dv x ( t ) = x 0 + v 0 t + 1 2 a 0 t 2 � Do the analysis as before: � No problem to get v(t) right: � � � � � � just need 1st order accuracy x 0 I x d � = � � � � � � But x(t) demands 2nd order accuracy � K � D � � � � � � dt v v � So we can look at mixed methods: • 1st order in v � Eigenvalues of this matrix? • 2nd order in x cs533d-term1-2005 11 cs533d-term1-2005 12
More Approximations… Simplification � � is the eigenvalue of the Jacobian, and � Typically K and D are symmetric semi-definite 2 � � K (there are good reasons) ( ) � = � 1 2 � D ± 2 � D 1 • What does this mean about their eigenvalues? � � � Often, D is a linear combination of K and I � Same as eigenvalues of 0 1 (“Rayleigh damping”), or at least close to it � � � � K � � D � � • Then K and D have the same eigenvectors (but different eigenvalues) � Can replace K and D (matrices) with • Then the eigenvectors of the Jacobian are of the form corresponding eigenvalues (scalars) (u, � u) T • Just have to analyze 2x2 system • [work out what � is in terms of � K and � D ] cs533d-term1-2005 13 cs533d-term1-2005 14 Split Into More Cases Three Test Equations � Constant acceleration (e.g. gravity) � Still messy! Simplify further • a(x,v,t)=g � If D dominates (e.g. air drag, damping) • Want exact (2nd order accurate) position { } � � � � D , 0 � Position dependence (e.g. spring force) • a(x,v,t)=-Kx • Want stability but low or zero damping • Exponential decay and constant • Look at imaginary axis � If K dominates (e.g. spring force) � Velocity dependence (e.g. damping) • a(x,v,t)=-Dv • Want stability, monotone decay � � ± � 1 � K • Look at negative real axis cs533d-term1-2005 15 cs533d-term1-2005 16
Explicit methods from before Implicit methods from before � Backward Euler � Forward Euler • Constant acceleration: bad (1st order) • Constant acceleration: bad (1st order) • Position dependence: very bad (unstable) • Position dependence: ok (stable, but damps) • Velocity dependence: ok (conditionally • Velocity dependence: great (monotone) monotone/stable) � Trapezoidal Rule � RK3 and RK4 • Constant acceleration: great (2nd order) • Constant acceleration: great (high order) • Position dependence: great (stable, no • Position dependence: ok (conditionally stable, but damping) damps out oscillation) • Velocity dependence: good (stable but only • Velocity dependence: ok (conditionally conditionally monotone) monotone/stable) cs533d-term1-2005 17 cs533d-term1-2005 18 Setting Up Implicit Solves Symmetry � Let � s take a look at actually using Backwards � Why multiply by M? � F position � F Euler, for example � Physics often demands that and velocity x n + 1 = x n + � t v n + 1 � x � v are symmetric ( ) • And M is symmetric, so this means matrix is v n + 1 = v n + � t M � 1 F x n + 1 , v n + 1 symmetric, hence easier to solve � Eliminate position, solve for velocity: • (physics generally says matrix is SPD - even better) ( ) v n + 1 = v n + � t M � 1 F x n + � t v n + 1 , v n + 1 • If the masses are not equal, the acceleration form of � Linearize at guess v k , solving for v n+1 � v k + � v the equations results in an unsymmetric matrix - bad. � F � � ) + � t � F � x � v + � F � Unfortunately the matrix is usually v k + � v = v n + � t M � 1 F x n + � tv k , v k ( velocity � � v � v � unsymmetric � x � � • Makes solving with it considerably less efficient � Collect terms, multiply by M • See Baraff & Witkin, “Large steps in cloth simulation”, � � M � � t � F � v � � t 2 � F ( ) + � t F x n + � tv k , v k ( ) � � v = M v n � v k � SIGGRAPH � 98 for one solution: throw out bad part � � x � cs533d-term1-2005 19 cs533d-term1-2005 20
Specialized 2nd Order Methods Symplectic Euler � Like Forward Euler, but updated velocity used � This is again a big subject for position � Again look at explicit methods, implicit ( ) v n + 1 = v n + � ta x n , v n methods x n + 1 = x n + � tv n + 1 � Also can treat position and velocity dependence differently: � Some people flip the steps (= relabel v n ) mixed implicit-explicit methods � Symplectic means certain qualities of the underlying physics are preserved in discretization - quite desirable visually! � [work out test cases] cs533d-term1-2005 21 cs533d-term1-2005 22 Symplectic Euler performance Tweaking Symplectic Euler � [sketch algorithms] � Constant acceleration: bad • Velocity right, position off by O( � t) � Stagger the velocity to improve x � Start off with � Position dependence: good ( ) v 12 = v 0 + 1 2 � ta x 0 , v 0 � t < 2 • Stability limit � Then proceed with K ( ) v n + 12 = v n � 12 + 1 2 ( t n + 1 � t n � 1 ) a x n , v n � 12 • No damping! (symplectic) x n + 1 = x n + � tv n + 12 � Velocity dependence: ok • Monotone limit � t < 1 D � Finish off with ( ) • Stability limit � t < 2 D v N = v N � 12 + 1 2 � ta x N , v N � 12 cs533d-term1-2005 23 cs533d-term1-2005 24
Recommend
More recommend
