Dynamical model basics
1. Models and differential equations
In physics, a model is an explicit statement of rules or equations whose outputs or variables correspond to measurable quantities. The significance of models is that they allow us to make explicit predictions about the future behavior of a system.
Models in classical physics are usually fleshed out in the form of differential equations.
Differential equation (DE): any equation referring to the some variable x and its derivatives (dx/dt, d2x/dt2, …). Any equation states an invariant law or relationship between the variable x and the rest of the constants or variables in the equation.
2. Terms: solutions, state (of system), phase space or plane, trajectory
3. First-order ordinary DEs take the general form dx/dt = f(x)
This is a subclass of DEs. First-order means there are no second derivatives and ‘ordinary’ means there are no partial derivatives or that all derivatives are wrt a single variable x. The basic method of solution is via integration.
4. Examples
growth-decay models with constant rate:
dx/dt = w, w 0 x(t) = x(0) + wt
growth-decay models with exponential rate:
dx/dt = ax x(t) = x(0)eat
Verhulst or Logistic growth equation:
dx/dt = kx – ax2x(t) = kx(0)/(k-ax(0))e-kt+ ax(0)
with these examples at hand, we can now define some descriptive terms …
5. Solution = an expression for x(t) that satisfies the equation (leads to identity)
6. State (of system) = value of x at time t, x(t)
7. Phase space = set of all possible states; here the real axis
8. Trajectory = a plot of x(t) in the phase space, a possible evolution of the dynamics
Note that for any given initial condition there is only one trajectory prescribed by the solution to our equation; but the DE tells us about all possible trajectories, e.g., whether these trajectories are lines or exponential curves or S-like curves and so on.
So, it would be good to be able to also represent this more global view of DEs, hence the notion of phase flow.
9. Phase flow, vector field, phase portrait = pattern of trajectories in the phase space
Take a handful of xi values and for each such value draw an arrow of length proportional to |f(xi)| on the x-axis, with its center at xi, and pointing in the direction of increase or decrease depending on the sign of f(xi); if f(xi) is positive the arrow points to the right, if it’s negative it points to the left.
Draw the phase flows for dx/dt = w, dx/dt = ax
dx/dt = w uniform flow with monotonic increase or decrease depending on the sign of w
for w > 0
for w < 0
dx/dt = ax[a1] flow increases monotonically away from 0 for a > 0 (corresponding to an exponential growth solution); flow decreases monotonically toward 0 for a < 0 (… exponential decay)
for a > 0
for a < 0
Intuitive interpretation: the flow of the DE is analogous to the flow of water; a single trajectory is analogous to the path taken by a (massless) particle if it had been placed in the water.
10. Fixed point of the flow (or zero of the vector field)
The points xk where [dx/dt = ] f(xk) = 0 represent states of equilibrium – when our particle is placed initially at xk it remains there for all time (vs. at all other points where the state of the system changes). Such points are called fixed points.
11. Geometric description of fixed points
We can obtain important information about the solutions to a differential equation without actually integrating the equation.
Consider the DE dx/dt = f(x) and the graph of f(x), a continuous function of x. We are interested in the evolution of the system starting at some initial position in the phase space x = a. Since dx/dt = f(x), then if f(a) > 0, x will increase; if f(a) < 0, then when we start at x = a, x will decrease. This increase or decrease will take place until we arrive at a point where x where f(x) (=dx/dt) = 0. This is the fixed-point. In general, determining the fixed points reduces to finding the roots of the equation f(x) = 0.
12. Example: consider the equation dx/dt = f(x) = -x + x3 and plot the graph of f(x) using Matlab [AG2]or your favorite software.
- Determine the fixed points algebraically
- Use the graph of f(x) (first panel below) to draw the flow
13. Intuitive visualization of fixed points – notion of Potential
For first-order systems, we can express f(x) as the negative gradient of a potential function V(x). This means that f(x) = -dV(x)/dx, which implies that for some constant term C we can write V(x) = -f(x)dx + C.
Then, the states of the system can be imagined by placing a ball in the potential V(x). The position of the ball (the state of the system) flows downhill away from the maxima of V(x).
This is illustrated above with the graph of V(x) for our f(x) = -x + x3. This V(x) = (x^2)/2 - (x^4)/4.
14. Two types of fixed points
stable:
around which f(x) is a decreasing function of x, or intuitively, the arrows of flow point towards that point
unstable:
around which f(x) is an increasing function of x, or intuitively, the arrows of flow point away from that point
15. Stability analysis of fixed points
Assume the dynamical system dx/dt = f(x)
Given a fixed point, we perform a Taylor expansion of f(x) in the vicinity of x. That is, we write f(x) = f(x*) + df/dx(x=x*) (x-x*) + ½ d2f/dx2(x=x*) (x-x*)2 + … . Because f(x*) = 0 and because very close to x* the term (x-x*)2 and all the higher-order terms is much smaller than (x-x*), we can linearize by rewriting f(x) = m(x-x*) where m = df/dx(x=x*).
Now, if we perform the minor variable substitution y = x – x*, then dy/dt = dx/dt (as x* is a constant) and we can rewrite f(x) = my. So our original equation takes the simpler form dy/dt = my.
But this is our familiar equation from exponential growth or decay. Its behavior is well understood. There is one fixed point at y = 0 (or x = x* given that y = x-x*). If m > 0, then there is monotonic exponential departure from the fixed point. If m < 0, then there is monotonic exponential approach to the fixed point. In the former case, we have an unstable fixed point and in the latter a stable fixed point.
Indeed, the equation dx/dt = ax for a < 0 is particularly important because it describes the motion in a sufficiently small neighborhood of most fixed points.
16. By way of summary, we have seen two approaches to the study of DEs:
algebraic or analytic solution = find an explicit expression for x(t);
geometric solution = find the fixed points and determine their stability
Simulating DEs
1. Useful guides for solving DEs with simulation
Higham (2001) is a good source and tells you how to solve DEs using Matlab scripts. Brown et al. (2006)is a discussion of DE solvers in the context of cognitive psychology models.
2. Euler’s method of solving DEs
Devised by Euler in 1768. It is simple, well-understood and almost everywhere used in psychology models. Better techniques exist that minimize error (between the true solution and the estimated one) but if you follow some good practice (see articles above) the simple method will do just fine.
Here is the logic of Euler’s method. Take the equation of our model, dx/dt = f(x). Combine that with the well-known formula for a derivative fromCalculus. Then f(x) = dx/dt = x(t+h)-x(t)/h from which we can solve for x(t+h).
x(t+h) = x(t) + h * f(x) OR new = old + change in value
Change in value = step size (=h) * slope of x(t)
3. Euler’s method, step by step
Suppose you have know the initial value of x (say x = c) at some time point t, say t = a, and you want to find out the value of x at some later time point, say t = b. Divide (b-a) in n equal increments, so that h = (b-a)/n. This is the step size. Then by our equation dx/dt = f(x), we can see that the change in x would be dx=dt * f(x). Then add this change to the current value of x (=c) to get the next value. Repeat to get the values of all subsequent time steps.
txchange
______
ach * f(a)
a+hc+ h * f(c)h * f(a+h)
a+2h…
.…
.…
.…
b???xxx
4. Two simple Matlab files implement the simple step-by-step Euler method (copy and paste to Matlab)
% Filename: f.m (name must match with the name of the function)
function y = f(x)
y = -x;
% Filename: ef.m
clear
t0 = 0; % Start time
tfinal = 10; % Final time
x0 = 1; % Initial value (x(t0))
h = 0.05; % Time step
N = ceil( (tfinal-t0)/h ); % Computes the number of loop steps
t(1) = t0; % t( ) will be a vector with the time grid
x(1) = x0; % x( ) will be a vector with the "solution" for each time value
for k=1:N
t(k+1) = h*k; % Fills the time grid vector
x(k+1) = x(k)+h*f(x(k)); % Euler forward iteration
end
plot(t,x)
5. Adding Noise and variability
So far we have looked at deterministic dynamical models. To enter the world of stochastic dynamical models, we add noise on the right-hand side of the deterministic model, that is, dx/dt = f(x) + noise. In Euler’s method terms, the solution will be given by the following equation.
x(k+1) = x(k) + dt * f(x(k)) + dt * noise
Note that dt also multiplies the noise factor (solve for dx in dx/dt = f(x) + noise).Consider the effects of the added noise factor relative to the dt * f(x(k)). Consider in particular regions around the fixed points. Around these, f(x) will be very small. Therefore around these points dt * f(x(k)) is very small compared to dt * noise and as a result these are the regions where we expect noise to have most of an effect greater variability.
6. Scripts: odesim.m and odesimhist.m
odesim.m - runs simulation for any number of different initial values with a given potential function. Outputs graphs of the potential used, the trajectory of each initial value given, and the flow field associated with the given potential.
odesimhist.m - runs simulation a given number of times using either a random initial value each time (the default), or a given specific initial value each time. Output is a graph of the potential used, and a histogram of final particle positions over all simulation runs.
Usage of both files involves changing any desired simulation parameters (in the SETTINGS section at the top of each file) and then just calling either odesim or odesimhist. The two files work for all potentials that can be represented by ax4/4 − bx2/2 + cx. This includes monostable, bistable symmetric, and bistable asymmetric potentials. Simulating different potential types is just done by changing the a,b,c coefficients, and is described in the top comments of the m files.
7. Some specifics for each script
** odesim.m, Noise: d ** Consider different effects of noise strength …
Set d = 0; no bumps in the trajectories -- deterministic behavior
Set d = 10; bumps, stochastic behavior
Set d = 100; particle trajectory is completely erratic
** odesim.m, Time scaling: tau ** dX(j-1) = dX(j-1)*tau; Higher tau values means faster convergence to attractors. Try 1, 2, 4 etc. to see this.
** odesimhist.m ** Convergence to attractor from RANDOM initial poitions.
8. Scripts odesimSTD.m, tts.m
** odesimSTD.m ** Visuals and estimation of std of solution under certain noise values in deep vs. shallow attractor. Specified, rather than random init pos.
** tts.m ** Compute time to settle from different initial values by integrating potential between initial value and attractor.
References
Amari, S. (1977). Dynamics of pattern formation in lateral-inhibition type neural fields. Biological Cybernetics, 27:77–87.
Brown, Scott D., Roger Ratcliff and Philip L. Smith. 2006. Evaluating methods for approximating stochastic differential equations, Journal of Mathematical Psychology, Volume 50, Issue 4, Pages 402-410.
Chomsky, N. (1965). Aspects of the theory of syntax. Cambridge, MA: MIT Press.
Erlhagen, W., & Sch¨oner, G. 2002. Dynamic Field Theory of Movement Preparation. Psychological Review, 109, 545–572.
Fodor, J. A., Bever, T., & Garrett, M. F. (1974). The psychology of language: An introduction to psycholinguistics and generative grammar. New York: McGraw-Hill.
Hawkins, J.A. (2004). Efficiency and Complexity in Grammars. Oxford: OxfordUniversity Press.
Higham, J. D. (2001). An algorithmic introduction to numerical simulation of stochastic differential equations . SIAM Review, 43, 525-546.
Moreno, M. A. (2002). A nonlinear dynamical systems perspective on response time distributions. Unpublished doctoral dissertation, ArizonaStateUniversity, Tempe.
Phillips, C. (1996). Order and structure. Ph.D. thesis, MIT, Cambridge, MA.
[a1]Example of analytic solution in Matlab in analytic.m script.:
subplot(1,2,1); % plot 1 of a 2 by 2 array of plots
t=[0:0.1:4]; % vector of time points
y=exp(-t); % solution vector for initial position 1
axis([0 4 -5 4]);
plot(t,y), hold on % plot this trajectory, stay on this plot
y=2*exp(-t);
plot(t,y), hold on
y=4*exp(-t);
plot(t,y), hold on
y=-3*exp(-t);
plot(t,y), hold on
y=-5*exp(-t);
plot(t,y), hold on
grid on
title('Trajectories');
xlabel('time');
subplot(1,2,2);
[t,y] = meshgrid(0.2:0.4:4,-5:0.4:4); % prepare points on a grid
pt = t .* 0 + 0.1; % t coordinate of vector at each point
py = -0.1 .* y; % y coordinate of vector at each point
quiver(t,y,pt,py) % plot the flow field
title('Flow field');
xlabel('time');
[AG2]1Matlab script to plot f(x) and its potential (in ExampleFlow.m)
x = -3:1/100:3;
y = -x + x.^3;
subplot(2,1,1);
plot(x, y, 'r-');
axis([-2 2 -0.5 0.5]);
title('Velocity function f(x) = -x + x^3');
v = x.^2/2 - x.^4/4;
subplot(2,1,2);
plot(x, v, 'r-');
axis([-2 2 -0.5 0.5]);
title('Potential: f(x) = -dV(x)/x => V(x) = (x^2)/2 - (x^4)/4');
% set graph line properties for displaying
h = findobj('Type', 'line'); % finds the line object and return a handle to it
set(h, 'LineWidth', 2); % using the handle you can change object properties