851 0585 04 L Modeling and Simulating Social

851 -0585 -04 L – Modeling and Simulating Social Systems with MATLAB Lecture 3 – Dynamical Systems Karsten Donnay and Stefano Balietti Chair of Sociology, in particular of Modeling and Simulation © ETH Zürich | 2011 -03 -07

Schedule of the course Introduction to MATLAB 21. 02. 28. 02. 07. 03. 14. 03. 21. 03. 28. 03. Working on projects (seminar theses) 04. 18. 04. 02. 05. 09. 05. 16. 05. 23. 05. 30. 05. 2011 -03 -07 Introduction to social-science modeling and simulation Handing in seminar thesis and giving a presentation (final deadlines to be communicated) K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 2

Schedule of the course Introduction to MATLAB 21. 02. 28. 02. 07. 03. 14. 03. 21. 03. 28. 03. Working on projects (seminar theses) 04. 18. 04. Introduction to social-science modeling and simulation 02. 05. 09. 05. 16. 05. 23. 05. 30. 05. 2011 -03 -07 Create and Submit a Research Plan Handing in seminar thesis and giving a presentation (final deadlines to be communicated) K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 3

Goals of Lecture 3: students will 1. Consolidate knowledge acquired during lecture 2, through brief repetition of the main points and revision of the exercises. 2. Understand fully the importance and the role of a Research Plan and its main standard components. 3. Learn the basic definition of Dynamical Systems and Differential Equations. 4. Implement Dynamical System in MATLAB a series of examples from Physics and Social Science Literature (Pendulum, Lorenz attractor, Lotka-Volterra equations, Epidemics: Kermack-Mc. Kendrick model) 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 4

Repetition § § § plot(), loglog(), semilogy(), semilogx() § Log-scale plots reduce the complexity of some functional forms and fit well to draw large range of values § The amplitude of the sample matters! law of large numbers – LLN. (Ex. 1) § Outlier observations can significantly skew the distribution, e. g. mean != median. (Ex. 2) Cov() and corrcoef() can measure correlation among variables. (Ex. 3) § hist(), sort() subplot() vs figure() 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 5

Repetition – Ex. 4 >> loglog(germany(: , 2), 'Line. Width', 3); 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 6

Project § Implementation of a model from the Social Science literature in MATLAB § During the next two weeks: § § § 2011 -03 -07 Form a group of two or three persons Choose a topic among the project suggestions (available on line) or propose your own idea Fill in the Research Plan (available on line) and send it to sbalietti@ethz. ch AND kdonnay@ethz. ch by email Subject line: [MATLAB FS 11] name 1, name 2, name 3 E. g. [MATLAB FS 11] donnay, balietti K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 7

Research Plan Structure § 1 -2 (not more!) pages stating: § Brief, general introduction to the problem § Fundamental questions you want to try to answer § Existing literature you will base your model on and possible extensions § Research methods you are planning to use § Use a Word/Openoffice format, to easily track changes § It will become the first page of your report. 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 8

Research Plan § Upon submission, the Research Plan can: be accepted § be accepted under revision § be modified later on only if change is justified (create new version) § § Talk to us if you are not sure § Remember: the sooner the better! § Final deadline for signing-up for a project is: March 14 th 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 9

“Plans are useless, but planning is everything. ” Dwight David Eisenhower (14 Oct 1890 – 28 Mar 1969) 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 10

Dynamical systems § Mathematical description of the time dependence of variables that characterize a given problem/scenario in its state space. 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 11

Dynamical systems § A dynamical system is described by a set of linear/non-linear differential equations. § Even though an analytical treatment of dynamical systems is usually very complicated, obtaining a numerical solution is (often) straight forward. § Differential equation and Difference equation are two different tools for operating with Dynamical Systems 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 12

Differential Equations § A differential equation is a mathematical equation for an unknown function (dependent variable) of one or more (independent) explicative variables that relates the values of the function itself and its derivatives of various orders § Ordinary (ODE) : § Partial (PDE) : 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 13

Numerical Solutions for Differential Equations § Solving differential equations numerically can be done by a number of schemes. The easiest way is by the 1 st order Euler’s Method: x x(t) x(t-Δt) t Δt 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 14

A famous example: Pendulum § A pendulum is a simple dynamical system: L = length of pendulum (m) ϴ = angle of pendulum g = acceleration due to gravity (m/s 2) The motion is described by: 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 15

Pendulum: MATLAB code 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 16

Set time step 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 17

Set constants 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 18

Set starting point of pendulum 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 19

Time loop: Simulate the pendulum 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 20

Perform 1 st order Euler’s method 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 21

Plot pendulum 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 22

Set limits of window 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 23

Make a 10 ms pause 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 24

Pendulum: Executing MATLAB code 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 25

Lorenz attractor § The Lorenz attractor defines a 3 -dimensional trajectory by the differential equations: § σ, r, b are parameters. 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 26

Lorenz attractor: MATLAB code 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 27

Set time step 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 28

Set number of iterations 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 29

Set initial values 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 30

Set parameters 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 31

Solve the Lorenz-attractor equations 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 32

Compute gradient 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 33

Perform 1 st order Euler’s method 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 34

Update time 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 35

Plot the results 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 36

Animation 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 37

2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 38

Food chain 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 39

Lotka-Volterra equations § The Lotka-Volterra equations describe the interaction between two species, prey vs. predators, e. g. rabbits vs. foxes. x: number of prey y: number of predators α, β, γ, δ: parameters 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 40

Lotka-Volterra equations § The Lotka-Volterra equations describe the interaction between two species, prey vs. predators, e. g. rabbits vs. foxes. 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 41

Lotka-Volterra equations § The Lotka-Volterra equations describe the interaction between two species, prey vs. predators, e. g. rabbits vs. foxes. 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 42

Epidemics Source: Balcan, et al. 2009 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 43

SIR model § A general model for epidemics is the SIR model, which describes the interaction between Susceptible, Infected and Removed (immune) persons, for a given disease. 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 44

Kermack-Mc. Kendrick model § Spread of diseases like the plague and cholera? A popular SIR model is the Kermack-Mc. Kendrick model. § The model assumes: A constant population size. § A zero incubation period. § The duration of infectivity is as long as the duration of the clinical disease. § 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 45

Kermack-Mc. Kendrick model § The Kermack-Mc. Kendrick model is specified as: S: Susceptible persons I: Infected persons R: Removed (immune) persons β: Infection rate γ: Immunity rate 2011 -03 -07 S β transmission R γ recovery K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch I 46

Kermack-Mc. Kendrick model § The Kermack-Mc. Kendrick model is specified as: S: Susceptible persons I: Infected persons R: Removed (immune) persons β: Infection rate γ: Immunity rate 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 47

Kermack-Mc. Kendrick model § The Kermack-Mc. Kendrick model is specified as: S: Susceptible persons I: Infected persons R: Removed (immune) persons β: Infection rate γ: Immunity rate 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 48

Exercise 1 § Implement and simulate the Kermack. Mc. Kendrick model in MATLAB. Use the starting values: S=I=500, R=0, β=0. 0001, γ =0. 01 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 49

Exercise 2 § A key parameter for the Kermack-Mc. Kendrick model is the epidemiological threshold, βS/γ. Plot the time evolution of the model and investigate the influence of the epidemiological threshold, in particular the cases: 1. 2. βS/γ < 1 βS/γ > 1 Starting values: S=I=500, R=0, β=0. 0001 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 50

Exercise 3 - optional § Implement the Lotka-Volterra model and investigate the influence of the timestep, dt. § How small must the timestep be in order for the 1 st order Euler‘s method to give reasonable accuracy? § Check in the MATLAB help how the functions ode 23, ode 45 etc, can be used for solving differential equations. 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 51

References § Matlab and Art: http: //vlab. ethz. ch/ROM/DBGT/MATLAB_and_art. html § Kermack, W. O. and Mc. Kendrick, A. G. "A Contribution to the Mathematical Theory of Epidemics. " Proc. Roy. Soc. Lond. A 115, 700 -721, 1927. § "Lotka-Volterra Equations”. http: //mathworld. wolfram. com/Lotka. Volterra. Equations. html § § http: //en. wikipedia. org/wiki/Numerical_ordinary_differential_equations "Lorenz Attractor”. http: //mathworld. wolfram. com/Lorenz. Attractor. html 2011 -03 -07 K. Donnay & S. Balietti / kdonnay@ethz. ch sbalietti@ethz. ch 52