Scientific computing | Week 9 Linear algebra and dynamical systems Volker Roth, Ivan Dokmani ć 1
Why study differential equations? • But differential equations are so 20th century :-( Life sciences Physics Environment Engineering Finance {Epi, Pan}demics Planetary dynamics CO 2 concentration Heat transfer Black–Scholes Neuroscience Turbulence Ice melting Trajectory design Market crashes 2
Recap: linear ODEs with linear algebra A system of first-order linear ODEs (homogeneous, constant-coefficient) d u u ( t ) ∈ ℝ n A ∈ ℝ n × n dt = A u Entries of model positions, velocities, CO 2 concentrations, … u Scalar case ( ) Vector case n = 1 u ( t ) = e α t u (0) u ( t ) = e At u 0 u ′ ( t ) = α u ( t ) 3
The meaning of e At Defined via the Taylor (power) series ∞ ( At ) k e At := ∑ k ! k =0 Use to check the solution ∞ ∞ ∞ kt k − 1 A k t k − 1 A k − 1 t k A k ∑ ∑ ∑ ( e At ) ′ := = Ae At = A ( k − 1)! = A k ! k ! k =1 k =1 k =0 For diagonalizable matrices , we get a very simple rule λ 1 0 ⋯ 0 e λ 1 t 0 ⋯ 0 e λ 2 t ⋱ 0 λ 2 0 0 ⋱ 0 e At = V ⟹ V − 1 V − 1 A = V 0 ⋱ ⋱ 0 0 ⋱ ⋱ 0 e λ n t 0 ⋯ 0 λ n 0 ⋯ 0 4
Behavior of first-order equations for n = 1 When is a real number, α u ( t ) = e α t u (0) u ′ ( t ) = α u ( t ) In this case the possible dynamics are quite boring… (but they can also be dangerous!) 5
Not so boring: second-order differential equations mx ′ ′ + bx ′ + kx = 0 + x − kx θ i 0 ∆ x u mg 6
Mass on a spring Natural spring Equilibrium position mg = ks ACME ACME ACME ACME F G = mg F S = − k ( s + x ) 7
Mass on a spring Gravitational force Restoring force in the spring F G = mg F S = − k ( s + x ) Newton says F tot = ma = mx ′ ′ F tot = F S + F G ⟹ mx ′ ′ + kx = 0 ′ = − ω 2 x k = ω 2 x ′ A second-order linear ODE! Converting to first order lets us use linear algebra: dt := [ ] = [ u := [ − ω 2 0 ] d u 0 1 x x x ′ [ ] ] x ′ x ′ x ′ ′ ⏟ 8 u A
Writing down the solution By solving we get the eigenvalues of as det( λ I − A ) = 0 A λ 1 = j ω λ 2 = − j ω Solving we further get the eigenvectors Av 1,2 = λ 1,2 v 1,2 v 1 = [ λ 1 ] = [ v 2 = [ λ 2 ] = [ j ω ] − j ω ] 1 1 1 1 Any solution can thus be written as (for some constants and ) c 1 c 2 u ( t ) = c 1 e j ω t v 1 + c 2 e − j ω t v 2 The constants and can be determined from two initial conditions (on and ) c 1 c 2 x x ′ u (0) = c 1 v 1 + c 2 v 2 9
We get the familiar harmonic oscillations x ( t ) = c 1 e j ω t + c 2 e − j ω t ℑ ( x ( t )) = 0 | t =0 ⟹ ℑ ( c 1 ) = − ℑ ( c 2 ) ℑ ( x ( t )) = 0 | t = π 2 ⟹ ℜ ( c 1 ) = ℜ ( c 2 ) ⟹ c 1 = c 2 So finally, for some real and (or and ) A B α ϕ x ( t ) = A cos( ω t ) + B sin( ω t ) = α sin( ω t + ϕ ) 10
Damped oscillations with F D = − bx ′ The force describes damping, friction, proportional to velocity F D = − bx ′ now gives the full second-order equaton F S + F G + F D = mx ′ ′ mx ′ ′ + bx ′ + kx = 0 Rewrite again as a first-order system dt := [ ] = [ − k / m − b / m ] d u x x ′ 0 1 [ ] x ′ x ′ ′ ⏟ u A 11
General solution Solving gives det( λ I − A ) = 0 b 2 − 4 mk b 2 − 4 mk λ 1 = − b + λ 2 = − b − 2 m 2 m General solution x ( t ) = c 1 e λ 1 t + c 2 e λ 2 t Behavior depends on the sign of the discriminant b 2 − 4 mk ≶ 0 12
Different kinds of solutions Overdamped Underdamped b 2 < 4 mk b 2 > 4 mk λ 1,2 < 0 ℜ ( λ 1 ) = ℜ ( λ 2 ) < 0 x ( t ) = c 1 e λ 1 t + c 2 e λ 2 t x ( t ) = e − α t ( A cos( ω t ) + B sin( ω t )) 13
Forced oscillations So far: the right-hand side of the ODE is zero: natural modes of the spring-mass system In most practical applications there is an external forcing • A voltage source in a circuit • Uneven road hits the wheels • Greenhouse gas emissions In our second-order linear case this is modeled as a right-hand side f ( t ) mx ′ ′ ( t ) + bx ′ ( t ) + kx ( t ) = f ( t ) 14
General solution to the non-homogeneous equation A general solution to can be written as mx ′ ′ ( t ) + bx ′ ( t ) + kx ( t ) = f ( t ) x ( t ) = c 1 x 1 ( t ) + c 2 x 2 ( t ) + x p ( t ) Solution to the homogeneous equation A particular solution mx ′ ′ + bx ′ + kx = 0 In a damped system, dictates the long-term behavior—steady-state solution x p 15
Forced oscillations: example x (0) = x ′ x ′ ′ + 8 x ′ + 16 x = 8 sin(4 t ) (0) = 0 Homogeneous e At = [ A = [ − e − 4 t (4 t − 1) ] e − 4 t (4 t + 1) te − 4 t − 16 − 8 ] 0 1 ⟹ − 16 te − 4 t x h ( t ) = c 1 e − 4 t + c 2 te − 4 t 16
Forced oscillations: total solution Particular Method of undetermined coefficients suggests to try x p ( t ) = A cos(4 t ) + B sin(4 t ) x p ( t ) = − 1 4 cos(4 t ) ⟹ Total 4 e − 4 t + te − 4 t − 1 1 4 cos(4 t ) (From Wikipedia) 17
Resonance Systems that are not overdamped have their own natural modes or resonant frequencies Example x ′ ′ ( t ) + x ( t ) = 5 cos( t ) • Homogeneous solution is x h ( t ) = A sin( t + φ ) • Method of undetermined coefficients gives https://sites.lsa.umich.edu/ksmoore/research/tacoma-narrows-bridge/ x p ( t ) = 1 2 t sin(1 t ) • The total solution is then x ( t ) = x h ( t ) + x p ( t ) = A sin( t + φ )+ 1 2 t sin(1 t ) 18
It happens even in real systems with damping! 19
Application 1: Predicting the CO 2 concentration in the atmosphere 20
The DICE model (Therina Groenewald/Shutterstock) (https://www.unenvironment.org/) The Dynamic Integrated Climate-Economy model Economics Carbon cycle Climate science Policy William Nordhaus, 2018 Nobel Prize in Economics Subject of quite a bit of controversy (likely a gross underestimate of the adverse effects) 21
Coupling between the CO 2 containers • An example of a box model : split the total CO 2 into boxes ( atmosphere , upper and lower ocean ) • The boxes exchange CO 2 with certain rates (often determined via experimental fitting) Atmosphere − k a k a Upper ocean k d − k d Lower ocean 22
Coupling between the CO 2 containers dM AT − k a k a = E ( t ) − k a ⋅ ( M AT − A ⋅ B ⋅ M UP ) dt dM UP = k a ⋅ ( M AT − A ⋅ B ⋅ M UP ) − k d ⋅ ( M UP − M LO δ ) − k d k d dt dM LO = k d ⋅ ( M UP − M LO δ ) dt , , model CO 2 mass in atmosphere, upper, and lower ocean (in gigaton) • M AT M UP M LO is the emission rate (gigaton / year) • E ( t ) is the equilibrium ratio of CO 2 between the atmosphere and the upper ocean • AB • is the volume ratio between upper and lower ocean δ , are CO 2 exchange rates between atmosphere/upper ocean and upper/lower ocean • k a k d 23
A linear algebra problem? d m dt = K m + e ( t ) M AT − k a k a AB 0 E ( t ) e ( t ) = M UP k a − k a AB − k d k d / δ m = K = 0 0 M LO 0 k d − k d / δ • Now we have an inhomogeneous system of ODEs • Is there a “principled” way to integrate (solve) such systems? 24
Solving the inhomogeneous equation d m We now know that when is constant in time, the solution to is dt = K m K m ( t ) = e Kt m (0) Duhamel’s principle Massage the inhomogeneous equation into a homogeneous form d m e tK d dt ( e − tK m ( t ) ) = e ( t ) dt = K m + e ( t ) ⟺ It follows that m ( t ) = e tK m (0) + ∫ t e ( t − s ) K e ( s ) ds 0 25
CO emission scenarios 2 • Representative concentration pathways (RCPs) : Emission scenarios from pre-industrial period to year 2050 • Consolidated by the Intergovernmental Panel on Climate Change (IPCC) • RCP 4.5 = intermediate emission; emissions peak in 2040 and then decline • RCP 8.5 = business-as-usual emissions; worst case 26
Approximating the integral by the trapezoidal rule t N 2 ( e ( t − j Δ t ) K + e ( t − ( j +1) Δ t ) K ) ∫ Δ t ∑ e ( t − s ) K ds ≃ 0 j =0 (Wikipedia) 27
CO2 concentration and the surface temperature • The temperature change with respect to a pre-industrial reference is estimated as α = 3.8 W/m 2 λ log 2 ( M AT , ref ) Δ T = α M AT λ = 1.3 W/m 2 / ∘ C M AT , ref = 596.4 GtC 28
Limitations of the model • A major limitation of the model is that is a constant matrix independent of time and K the current concentrations • In reality the carbonate chemistry dictates that the absorption capacity of the ocean drops after initial absorption, resulting in huge errors over longer timescales • One remedy is to allow the coefficients to depend on time and the k a , k d , AB , … current concentrations M AT , M UP , M LO • NB: Nordhaus’s work and models have been even more heavily criticized for how they measure economic utility, in that they “overemphasize growth as the ultimate measure of economic success” (https://www.sciencemag.org/news/2018/10/roles-ideas-and-climate-growth-earn-duo-economics-nobel-prize) 29
Recommend
More recommend
