CS205b/CME306Lecture 41 Linearized SystemForce is in general a function of both the positions and velocities of the particles in a system.That is, F = F(x, v). In the interests of writing down a linear system to analyse stability, it isconvenient to approximate this force as F (x, v) ≈ F (x0, v0) + Fx(x, v)x + Fv(x, v)v and also ignorethe inhomogeneous term F (x0, v0), since Duhamel’s principle states that the inhomogeneous termdoes not affect stability. Note that this approximation omits gravity, since it is inhomogeneous.We can typically make these simplifications when looking at stability. Forces of the f orm Fx(x, v)xare rather like spring forces, and forces of the form Fv(x, v)v behave like damping forces. UsingFx= maxand Fv= mav, the motion of the particle is described by the first order linear system xv!′= 0 1axav! xv!.The eigenvalues of this system areλ =av±pa2v+ 4ax2and have units of Hertz (s−1). Solutions look eλt, so well posedness requires the real part of theeigenvalues be nonpositive, Re(λ) ≤ 0. This places some restrictions on the way we can modelforces of nature to prevent the system from blow ing up.It is n ecessary for ax≤ 0. If a particle experiences no force at th e origin but experiences astronger force as it moves away, that force should be a restoring force rather than one that push itaway harder as it moves farther away and causes exponential growth.Similarly, it is necessary for av≤ 0. A force that did not satisfy this would tend to apply forcesin the direction of motion that get stronger as the particle moves faster and result in exponentialgrowth.When −av< 2√−ax, we call the system underdamped. The eigenvalues contain imaginarycomponents, and the solution exhibits period behavior. If av< 0, the system has exponentialdamping. If av= 0, the system is undamped. Note that in the undamped case, the eigenvalues arepure imaginary. In the underdamped case |λ| =√−axdoes n ot depend on the amount of dampingapplied.When −av> 2√−ax, we call the system overdamped. The eigenvalues are real and distinct,and the solution exhibits exponential decay only.When −av= 2√−ax, we call the system critically damped. The eigenvalues are equal, and thesolution exhibits exponential decay with at most one overshoot. A critically damped system decaysfaster than it would have with any other amount of damping. If it is u nderdamped, it overshoots1repeatedly and decays slowly. If it is overdamped, the excessive damping causes it to move towardsequilibrium slowly. Severe damping effectively freezes th e system so it can hardly move.2 Time IntegrationWe now consider several popular approaches for integrating an ODE.2.1 Forward EulerForward Euler evolution takes on the form xv!n+1= xv!n+ ∆t va!n.Because forward Euler is unstable when th e system contains pure imaginary eigenvalues, it isunstable without damping. In particular, av< 0. Forward Euler is stable when |∆tλ + 1| < 1. Atime step restriction of ∆ t ≈ |λ|−1is required. As problems get stiffer, m ax |λ| becomes larger, andconsequently the time step that can be taken becomes smaller. The stability region for forwardEuler is shown in Figure 1(a). If the system has eigenvalues n ear the imaginary axis, we wouldbe much better off using using a different method such as one of the higher order Runge Kuttamethods.2.2 Runge Kutta MethodsThere are a number of Runge Kutta methods in common use. T he Runge Kutta methods werediscussed in more detail in CS205A, so they are only highlighted here. The second order RungeKutta method is technically linearly unstable without d ampin g sin ce its stability region does notinclude the imaginary axis. Roundoff with av< 0 may be sufficient, and nonlinear effects cansometimes rectify the p roblem. The third order Runge Kutta method is stable without damping,since its stability region includes part of the imaginary axis. It is suitable for system with av= 0.Stability regions for various Runge Kutta methods are shown in Figure 1.2.3 Backward EulerThe need for a time step restriction is true of all explicit methods. By contract, implicit methodsoften have a much better time step restriction that makes them more suitable even though theyrequire a linear sys tem solve and are thus much more expensive per step.Backward Euler is the simplest of the implicit methods. Its stability region (Figure 2(a)) looksvery different from those of the explicit integrators. It takes the form xv!n+1= xv!n+ ∆t va!n+1.Note that the new values being obtained occur on both sides of th e equation, so a solve is typicallynecessary at each time step. Backward Euler has no time step restriction for stability, but largetime steps do result in poor accuracy.21 2−1−212−1−2(a) Forward Euler (First Order Runge Kutta)1 2−1−212−1−2(b) Second Order Runge Kutta1 2−1−212−1−2(c) Third Order Runge Kutta1 2−1−212−1−2(d) Fourth Order Runge KuttaFigure 1: Stability regions for Runge Ku tta. The scheme is s table with time step size ∆t if ∆tλ, acomplex number in general, is located within the shaded region of the complex plane. Note that thethird and f ourth order s chemes have stability regions that encompass the imaginary axis, makingthem stable on undamped problems. The first and second order schemes do not have this property.31 2−1−212−1−2(a) Backward Euler1 2−1−212−1−2(b) Trapezoid RuleFigure 2: Stability regions for two implicit time integration schemes. The scheme is stable withtime step size ∆t if ∆tλ, a complex number in general, is located within the shaded region of thecomplex plane.0 1 2 301e−t∆t = 4∆t = 2∆t = 1∆t =12∆t =14Figure 3: Effects of large time step sizes on numerical simulation of y′= −y with backward Euler.4The solve required for b ackward Euler can be relatively easy in the case of linear ODEs. Forexample, if the only force is d rag, F = −kv and a = −kv/m. Then, vn+1= vn− ∆tkvn+1/m sothat vn+1= (1 + ∆tk/m)−1vnand xn+1= xn+ ∆tvn+1= xn+ ∆t(1 + ∆tk/m)−1vnare relativelyeasily obtained. Note that this damps velocity each time step, and it damps more with larger timesteps. This is in contrast to forward Euler, where vn+1= (1 − ∆tk/m)vn. Damping here occurs if∆t < 2m/k = (max |λ|)−1but amplification occurs for larger time steps.More generally, backward Euler requires solving a linear system. If the ODE is nonlinear, itwill typically be more difficult and expensive to solve and