lucasdanb.org

Library Dynamica

#1 - Dynamical Systems

This is a web graphics library I created for analysis of two-dimensional continous dynamical systems. A dynamical system is something that changes state over time, expressed by a set of state variables. The systems that I'm concerned with have two state variables (hence "two-dimensional"). Continous means that each variable has a formula that calculates, based on the values of the current variables, what the instantaneous rate of change of that variable is.

A good example is two objects with different temperatures touching each other. Think of pouring hot coffee into a mug. Over time, the temperature of the coffee decreases, and the temperature of the mug increases. For now, let's ignore the temperature of the atmosphere around this coffee-mug system - we'll get to that later. Think about how fast each temperature changes too: the greater the difference between the two temperatures, the faster they converge toward each other. So, we could write a formula for the rate of change of the mug's temperature \(\dot m\) (pronounced m dot) based on its current temperature \(m\) and the temperature of the coffee \(c\): $$\dot m = c - m$$ As with the temperature of the air, let's ignore the surface area and heat transfer coefficient and all that. Oh yeah, and just like the formula above, the formula for the rate of change of the coffee's temperature \(\dot c\) would be \(\dot c = m - c\).

Now let's put this in dynamical system form. The rate of change of each variable is expressed as a function of the current state of the system. So our full system is: $$\dot c = f(c,\,m) = m - c$$ $$\dot m = g(c,\,m) = c - m$$

#1.1 - Visualization

These rates of change show us where the system is going. Imagine a set of coordinates \((c,\,m)\) at a single point in space, the current state of the system. Then, sticking out from that point, is an arrow that points to where the system is headed to, based on where it is now. In other words, the coordinates of the tip of that arrow are \((c + \dot c,\,m + \dot m)\). Using my library, I can graph the states of this system on a plane (called a vector field), each with one of these arrows. They're each scaled down a bit so the screen isn't cluttered:

Notice how all the possible states of the system converge toward that pink line? That's because the line is a stable equilibrium, an important concept in dynamical systems. An equilibrium occurs anywhere that the system is not changing. Think back to the formulas for our coffee-mug system... the rate of change of each variable is proportional to their difference. So if \(c = m\), their difference is 0 and the system is unchanging. Thus, the formula to graph that line is \(y = x\). And it makes sense when you think about it... the temperatures of the coffee and the mug converge toward each other, they meet in the middle of their initial temperatures.

#1.2 - Nullclines

Now let's add back in the temperate of the surrounding atmosphere that we so rudely ignored earlier. Technically, the atmosphere changes temperate slightly, but so little that we can ignore it again (sorry :/). So instead of adding a third dimension for the atmosphere, we'll just add a constant. Let's call it \(a\) and say it equals about 75 degrees Farenheit. Now our formulas have changed... let's take a look: $$\dot c = f(c,\,m) = \frac{m - c + a - c}{2}$$ $$\dot m = g(c,\,m) = \frac{c - m + a - m}{2}$$ What's happening now is we're taking an average of the differences between the coffee and the mug, and the coffee and the atmosphere. This looks a bit more interesting when graphed:

Now, all states converge to a single point: \((75,\,75)\), the temperature of the surrounding atmosphere. I've added a few purple trajectories to the graph, each a possible path that the system could travel along to eventually reach equilibrium. You can move your cursor around on the graph to see the trajectory from its point. There are a few more linear lines as well: the green line is the c-nullcline, and the blue line is the m-nullcline. Nullclines are sort of like equilibrium lines for only one variable. The c-nullcline is the line along which \(\dot c = 0\), so $$\frac{m - c + a - c}{2} = 0$$ Solved, its formula for graphing is \(m = 2c - 75\). Likewise, the m-nullcline is the line along which \(\dot m = 0\). It's a bit hard to tell with this system, but the arrows along each nullcline are either exactly horizontal or vertical, indicating that one state variable is not changing.

#1.3 - Equilibriums

To introduce the concepts of stable and unstable equilibria, and how to determine which is which, we'll need to take a look at linear transformations. Now, I admit that it is infinitely easier to classify equilibria by just looking at them on a graph. But, because this webpage doubles as my calculus notes, I'm going to explain the mathematical process here. First, let's take a look at types of equilibria:

Each of these four types can be classified as either attractors or repellers. There is also a fifth type called a saddle, which both attracts and repels. It is considered unstable.

#1.3.1 - Linear Transformations

The mathematical process for determining the stability of equilibria involves examining the eigenvalues of its Jacobian matrix. Up until now we've been thinking about these vector fields quite literally, with coordinates \((x,\,y)\) at the base of each arrow, and coordinates \((x + f(x,\,y),\,y + g(x,\,y))\) at the tip. And indeed, this is exactly how I've written the code to draw them. But we can also think of them as linear transformations (at least some of them). I'll show you what I mean.

This visual explanation of linear transformations is adopted from Grant Sandersons' YouTube video on the subject, and if this part is confusing, I recommend you watch his series on linear algebra. In essence, any transformation (spinning, stretching, scaling) of 2D space can be represented by a 2 x 2 matrix. Take a look at this one: $$\begin{bmatrix} 1 & k \\ k & 1 \\ \end{bmatrix}$$ Here, I'll let you play around with it in the graph below by scaling \(k\) and watching the rectangle compress and stretch.

\(k =\,\) $$\begin{bmatrix} 1 & k \\ k & 1 \\ \end{bmatrix}$$

Can you tell what kind of vector field this transformation would draw? Think about arrows going from the original coordinates of each point in the rectangle to the new coordinates after transformation. This will help:

\(k =\,\) $$\begin{bmatrix} 1 & k \\ k & 1 \\ \end{bmatrix}$$ $$\dot x = ky$$ $$\dot y = kx$$

This linear transformation creates a saddle. The most simple form of a saddle equilibrium can be created with the formulas \(\dot x = y,\,\dot y = x\), and in the graph above I've added the \(k\) term to each. Likewise, each type of equilibria can be described with a simple transformation matrix in this way. You can probably picture how stable or unstable nodes could be created by a rectangle that shrinks or grows. This interesting ability to turn an equilibrium of a dynamical system into a linear transformation allows us to analyze the properties of the system by analyzing the properties of the transformation. Right now, we're only concerned with the eigenvalues of the transformation, but there are likely countless additional interesting facts to be discovered about any dynamical system through analyzing it as a matrix.

#1.3.2 - Jacobian Matrices

Remember how I said that only some dynamical systems can be fully represented with a matrix? I'll show you what I mean. First, let's go through the process to transform the system $$\dot x = ky$$ $$\dot y = kx$$ into a linear transformation. Keep in mind the properties of matrix multiplication. We are trying to find a matrix \(M\) such that the \(x\) and \(y\) coordinate vector \(\vec v\) multiplied by \(M\) equals a vector of \(\dot x\) and \(\dot y\). It's easier to read if I make it bigger: $$\begin{bmatrix} \dot x \\ \dot y \end{bmatrix} = \begin{bmatrix} a & b \\ c & d \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\dot v = M\vec v$$ As we know from the matrix multiplication rules that I'm sure we all perfectly memorized from high school: $$\begin{bmatrix} x \\ y \end{bmatrix} \begin{bmatrix} a & b \\ c & d \end{bmatrix} = \begin{bmatrix} ax + by \\ cx + dy \end{bmatrix}$$ Maybe you see where this is going? What values could we plug in for \(a\) through \(d\) to create our system formulas in that final matrix? These values give us the solution: $$M = \begin{bmatrix} 0 & k \\ k & 0 \end{bmatrix}$$ Wait, that doesn't look right. The transformation matrix in the graphs above had 1s instead of 0s. Remember: this matrix gives us \(\dot v = M\vec v\). This is the rate of change of the coordinate vector; the actual coordinates of the tip of each arrow are \(\vec v + \dot v\). Obviously, for formulas \(x_1 = x + ky,\,y_1 = kx + y\), we get the actual visual transformation matrix.

The reason why this system was so easy to describe as a matrix was because it is linear. All this means is that, quite literally, the formulas that describe the rate of change of the system vector are straight lines. Remember the nullclines from earlier in our coffee-mug system? Those lines are straight, and therefore that system is linear. It's also worth noting that linear systems only have one equilibrium - two non-parallel straight nullclines can only intersect once. Therefore, the stability of the single equilibrium can be analyzed by analyzing the matrix of the entire system. For nonlinear systems with exponential nullclines and multiple equilibria, we instead look at the transformation matrix of each equilibrium. This is done by finding linear systems that approximate the behavior of each equilibrium, and then finding their transformation matrices through the process above. These matrices are called Jacobian matrices.

Let's take a look at the system graphed on the left, which has the formulas: $$\dot x = f(x,\,y) = -x^2 - y + 5$$ $$\dot y = g(x,\,y) = \frac{x}{2} - y - 3$$ and equilibriums at \((-3.09,\,-4.54)\) and \((2.59,\,-1.71)\). Starting with the first one, we translate it into a linear system $$\dot x = ax + by$$ $$\dot y = cx + dy$$ where \(a\) through \(d\) are the partial derivates of the system functions evaluated with the equilibrium coordinates \((x_0,\,y_0)\) like so: $$a = \frac{\partial f}{\partial x}(x_0,\,y_0) = -2x_0 = 6.18 \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, b = \frac{\partial f}{\partial y}(x_0,\,y_0) = -1$$ $$c = \frac{\partial g}{\partial x}(x_0,\,y_0) = \frac{1}{2} \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, d = \frac{\partial g}{\partial y}(x_0,\,y_0) = -1$$ which gives us the linear transformation matrix $$L = \begin{bmatrix} 6.18 & -1 \\ 0.5 & -1 \end{bmatrix}$$ I promise I'll get back to visual explanations soon, just bear with me. From \(L\) we can compute eigenvectors and eigenvalues. I'm not gonna bother to show you how to do that because, honestly, I just used Symbolab. All you need to know is what these concepts represent visually. The eigenvalues and eigenvectors of this transformation are: $$\lambda_1 = 6.11,\,\lambda_2 = -0.93 \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \vec v_1 = \begin{bmatrix} 14.22 \\ 1 \end{bmatrix},\, \vec v_2 = \begin{bmatrix} 0.14 \\ 1 \end{bmatrix}$$ These eigenvectors fall upon two lines on the graph of this system. Any vector on one of these lines will stay directly on it as the transformation takes place. Again, I turn to Grant Sanderson to explain this better. Let's switch back to the first linear transformation that we graphed, and I'll add in the eigenvectors to show you what I mean.

\(k =\,\) $$\begin{bmatrix} 1 & k \\ k & 1 \\ \end{bmatrix}$$ $$\dot x = ky$$ $$\dot y = kx$$

See how any arrows who's bases sit on either of the eigenvector lines only extend further along those lines as the transformation takes place? Those two lines are the only lines that pass through the origin with those properties. Try to draw another one, you'll notice that any arrows originating from the line will point further and further off of it. Now we finally have the foundation to understand what the eigenvalues mean. Each eigenvector line has a corresponding eigenvalue and it represents the factor by which any arrows on it are scaled. The orange line has an eigenvalue of \(1 + k\). The actual value doesn't really matter, only the fact that it is positive. This means that any arrow along that line will extend outward as the transformation takes place. Likewise, the eigenvalue of the blue line is \(1 - k\), and arrows along it extend inward. This is the conclusion of this section on eigenvalues - by analyzing the signs of the eigenvalues of the Jacobian matrix of the equilibrium of a dynamical system, you can determine its stability. Negative signs mean all states extend away from the equilibrium, and it is unstable. Positive signs mean all states converge toward it.

Now that we understand this property, we know that the leftmost equilibrium in our example nonlinear system is a saddle. If we follow the process outlined to determine the eigenvalues of the rightmost equilibrium, we'd find that it's a stable node.

#1.4 - Solutions

Apologies for the vertical learning curve with that previous section. We can move on to some cool examples now. The last capability of my library that I want to show off is the ability to solve these dynamical systems. To solve a system of differential equations would be to take, for example, our coffee-mug system $$\dot c = f(c,\,m) = m - c$$ $$\dot m = g(c,\,m) = c - m$$ and integrate \(f\) and \(g\) to get \(c\) and \(m\) as functions of time, based on their inital values. The actual mathematical process of doing this is incredibly complicated, involving lots of powers of \(e\), and even if I knew how to do it, I don't think I'd want to explain it. My library takes care of this by simply taking the current system state, moving forward a tiny amount, recalculating where to go next, and moving again. That's how you're able to create those purple trajectories from your cursor on the graphs. My library can also take the values of those movements and print them onto another graph over time, like below. Click somewhere on the large graph to launch a trajectory and watch each value graph itself over time on the two smaller graphs, until the system has stopped changing.

#2 - Code

The source for this library is on Github. If you got to this page from a link in my resume, check it out for proof that I can actually write good code. The actual bundled distribution is on NPM, and I've linked it in the head of this document through a CDN. The README on Github and NPM explains how to use the library in your website if you're interested.

#3 - Cool Examples

#3.1 - Van der Pol Oscillator

#3.2 - Predator and Prey Populations (Lotka - Volterra)

#3.3 - Neuron Action Potentials

\(I =\,\)\(\,pA\)

#3.4 - Simple Pendulum

Updated March 3, 2023