Differential Equations x i +1 x i x i +1 = x i + x x Physics-based - PowerPoint PPT Presentation

Physics-based simulation x i Differential Equations x i +1 x i x i +1 = x i + ∆ x ∆ x Physics-based simulation Differential equations x i What is a differential equation? Differential equations describe the relation between an Newtonian laws unknown function and its derivatives gravity wind gust Classes of DE x i +1 elastic force Ordinary differential equation (ODE) . . x i . x i +1 = x i + ∆ x Partial differential equation (PDE) ∆ x

Ordinary differential equations Symbolic solutions An ODE is an equality involving a function and its Standard introductory differential equation courses derivatives focus on finding solutions analytically known function x ( t ) = f ( x ( t )) ˙ Linear ODEs can be solved by integral transforms time state of the derivative system Use DSolve[eqn,x,t] in Mathematica Differential equation: x = − kx ˙ What does it mean to “solve” an ODE? x = e − kt Solution: Numerical solutions Initial value problems In this class, we will be concerned with numerical In a canonical initial value problem, the behavior of the solutions system is described by an ODE and its initial condition: For each iteration: x = f ( x, t ) ˙ 1. use f to calculate the derivative as the direction to next state x ( t 0 ) = x 0 2. use this direction to approximate change over a ∆ x time interval ∆ t 3. increment x by to obtain new value ∆ x To solve x ( t ) numerically, we start out from x 0 and follow the changes defined by f thereafter Derivative function f is regarded as a black box: given numerical values for x and t , returns a numerical value for ˙ x

Vector field Integral curves x 2 The differential equation can be Any point on the curve indicates visualized as a vector field the integral of from f t 0 ˙ x = f ( x , t ) f ( x , t ) � t 0 f ( x , t ) dt x How does the vector field look p t 0 like if f depends directly on time? x 1 Numerical methods Numerical methods Get confused yet? Explicit methods What do we know about the system? Explicit Euler’s method The midpoint method What are we trying to solve? Runge-Kutta method Implicit methods Implicit Euler’s method

Explicit Euler ’s method Problems of Euler ’s method Inaccuracy How do we get to the next state from the current state? x ( t 0 + h ) = x 0 + h ˙ x ( t 0 ) The circle turns into a spiral no x ( t 0 + h ) matter how small the step size is Instead of following real integral x ( t 0 ) curve, p follows a polygonal path p Discrete time step h determines the errors Problems of Euler ’s method Performance of Euler ’s method Cost per step is determined by the number of evaluations Instability per step f = − k x At each step, x ( t ) can be written in the form of Taylor Oscillation: series: Divergence: x ( t 0 ) + h 2 x ( t 0 ) + h 3 3! x (3) ( t 0 ) + . . . h n ∂ n x x ( t 0 + h ) = x ( t 0 ) + h ˙ 2! ¨ n ! ∂ t n How small the step size has to be? What is the order of the error term in Euler’s method?

The midpoint method Performance of midpoint method Prove that the midpoint method is correct within O(h 3 ) 1. Compute an Euler’s step ∆ x = hf ( x ( t 0 )) 2. Evaluate f at the midpoint f mid = f ( x ( t 0 ) + ∆ x 2 ) 3. Take a step using f mid x ( t 0 + h ) = x ( t 0 ) + hf mid x ( t + h ) = x 0 + hf ( x 0 + h 2 f ( x 0 )) Runge-Kutta method 4-stage 4 th order Runge-Kutta Runge-Kutta is a numeric method of integrating ODEs k 1 = hf ( x 0 , t 0 ) hf ( x 0 + k 1 2 , t 0 + h k 2 = 2 ) by using a trial step at the midpoint of an interval to hf ( x 0 + k 2 2 , t 0 + h k 3 = 2 ) cancel out lower-order error terms k 4 = hf ( x 0 + k 3 , t 0 + h ) x 0 + 1 6 k 1 + 1 3 k 2 + 1 3 k 3 + 1 x ( t 0 + h ) = 6 k 4 The second order Runge-Kutta is known as the midpoint method x 0 f ( x 0 , t 0 ) 1 . q -stage p -order Runge-Kutta evaluates the derivative f ( x 0 + k 1 2 , t 0 + h 2 . 2 ) function q times in each iteration and its approximation f ( x 0 + k 2 2 , t 0 + h x 3 . 2 ) x ( t 0 + h ) of the next state is correct within O(h p+1 ) f ( x 0 + k 3 , t 0 + h ) 4 . t

Stage vs. order Adaptive step size p 1 2 3 4 5 6 7 8 9 10 Ideally, we want to choose h as large as possible, but q min (p) 1 2 3 4 6 7 9 11 12-17 13-17 not so large as to give us big error or instability We can vary h as we march forward in time The minimum number of stages necessary for an explicit method to attain order p is still an open problem Step doubling Why is fourth order the most popular Runge Kutta Embedding estimate method? Variable step, variable order Step doubling Embedding estimate Estimate by taking a full Euler step x a x a = x 0 + hf ( x 0 , t 0 ) Also called Runge-Kutta-Fehlberg e Estimate by taking two half Euler steps x b Compare two estimates of x ( t 0 + h ) x temp = x 0 + h 2 f ( x 0 , t 0 ) Fifth order Runge-Kutta with 6 stages x b = x temp + h 2 f ( x temp , t 0 + h 2 ) Forth order Runge-Kutta with 6 stages e = | x a − x b | is bound by O ( h 2 ) � � � 1 2 h Given error tolerance , what is the optimal step size? � e

Problems of explicit methods Variable step, variable order Change between methods of different order as well Do not work well with stiff ODEs as step based on obtained error estimates Simulation explodes if the step size is too big These methods are currently the last work in Simulation progresses slowly if the step size is too small numerical integration Example: a bead on the wire Stiff equations y ( t ) Y ( t ) = ( x ( t ) , y ( t )) � � � � Y = d x ( t ) − x ( t ) ˙ = x ( t ) Stiffness constant: k y ( t ) − ky ( t ) dt Step size is limited by the largest k Explicit Euler’s method: Systems that has some big k’s mixed in are called � � � � x ( t ) − x ( t ) Y new = Y 0 + h ˙ Y ( t 0 ) = + h “stiff system” y ( t ) − ky ( t ) � � (1 − h ) x ( t ) Y new = (1 − kh ) y ( t )

Implicit methods Implicit method Our goal is to solve for Y new such that Explicit Euler: Y new = Y 0 + hf ( Y 0 ) Y new = Y 0 + hf ( Y new ) Implicit Euler: Approximating f ( Y new ) by linearizing f ( Y ) Y new = Y 0 + hf ( Y new ) , where f ( Y new ) = f ( Y 0 ) + ∆ Y f � ( Y 0 ) ∆ Y = Y new − Y 0 Solving for such that � , at time , points t 0 + h Y new f Y new = Y 0 + hf ( Y 0 ) + h ∆ Y f � ( Y 0 ) directly back at � � � Y 0 f ( Y , t ) = ˙ Y ( t ) � 1 � − 1 f � ( Y , t ) = ∂ Y ( t ) ∆ Y = h I − f � ( Y 0 ) f ( Y 0 ) ∂ Y Example: A bead on the wire Example: A bead on the wire Apply the implicit Euler’s method to the bead-on-wire What is the largest step size the implicit Euler’s method example can take? � 1 � − 1 h I − f � ( Y 0 ) ∆ Y = f ( Y 0 ) h � � h +1 x 0 h →∞ ∆ Y = lim lim � − x ( t ) h →∞ − h � 1+ kh ky 0 f ( Y ( t )) = − ky ( t ) � − 1 f � ( Y ( t )) = ∂ f ( Y ( t )) � � � � � x 0 0 x 0 = = − = − ∂ Y 0 − k 1 k ky 0 y 0 � − 1 � � 1+ h � − x 0 0 ∆ Y = h 1+ kh 0 − ky 0 h � � − x 0 Y new = Y 0 + ( − Y 0 ) = 0 � h � 0 h +1 = h 0 − ky 0 1+ kh h � � h +1 x 0 = − h 1+ kh ky 0

Implicit vs. explicit Second-order implicit Euler ¨ x = f ( x ( t ) , ˙ x ( t )) x d � � � � x ( t ) v ( t ) x ( h ) = − kx ( h ) ˙ = v ( t ) f ( x, ˙ x ) dt x (0) = 1 correct solution: x ( h ) = e − hk explicit Euler: x ( h ) = 1 − hk 1 . implicit Euler: x ( h ) = . 1 + hk . h � � � � I − h ∂ f ∂ v − h 2 ∂ f x 0 ) + h ∂ f ∆ v = h f ( x 0 , ˙ ∂ xv 0 ∂ x Helpful hints Modular implementation Write solvers in terms of Reusable solver code Simple model implementation Don’t use explicit Euler’s method Generic operations: Do use adaptive step size Get dim(x) Get/Set x and t Derivative evaluation at current (x, t)

Solver interface Summary How do we solve an ODE numerically? GetDim What are the drawbacks of Explicit Euler’s method? Get/Set What is a stiff system? System Solver State Why does implicit Euler allow us to take bigger Deriv Eval steps than explicit methods? What’s next? How do we use Euler’s method to simulate the motion of a mass point? What are the equations of motion when forces are applied to the mass point?