In bmarkslash7/QPot: Quasi-Potential Analysis for Stochastic Differential Equations

Calculating the local quasi-potential

Example 1 from Moore et al. (2016)

Conceptualizing the local quasi-potential

The next step is to compute a local quasi-potential for each basin of attractor (sys. with stable equilibrium). In the example above we have 2. For each stable equilibrium.

To understand the local quasi-potential, it is useful consider the analogy of a particle traveling according to the system of equations in 2. Stochastic simulation. In the context of example 1, the coordinates of the particle correspond to population densities, and the particle’s path corresponds to how those population densities change over time. The deterministic skeleton can be visualized as a force field influencing the particle’s trajectory. Suppose that the particle moves along a path from a stable equilibrium to another point $(x, y)$. If this path does not coincide with a solution of the deterministic skeleton, then the stochastic terms must be doing some ``work'' to move the particle along the path. The more work is required, the less likely it is for the path to be a realization of the system. $\Phi (x, y)$ is the amount of work required to traverse the easiest path from the stable equilibrium to $(x, y)$. Note that $\Phi (x, y)$ is non-negative, and it is zero at the stable equilibrium.

In the basin of attraction for esi, $\Phi (x, y)$ has many properties analogous to the potential function for gradient systems. Key among these properties is that the quasi-potential is non-increasing along deterministic trajectories. This means that the quasi-potential can be interpreted as a type of energy surface, and the rolling ball metaphor is still valid. The difference is that, in non-gradient systems, there is an additional component to the vector field that causes trajectories to circulate around level sets of the energy surface. This is discussed in more detail in [vignette on vector field decomposition], below.

QPot calculates quasi-potentials using an adjustment developed by Cameron (2012) to the ordered upwind algorithm (Sethian and Vladimirsky, 2001, 2003). The idea behind the algorithm is to calculate $\Phi (x, y)$ in ascending order, starting with the known the stable equilibrium. The result is an expanding area where the solution is known.

Calculating $\Phi (x, y)$ with the function QPotential requires:

1. a text string of the equations and parameter values,
2. the stable equilibrium points,
3. the computational domain, and
4. the mesh size.

The coordinates of the stable equilibrium points, which were determined in 1. Analyzing the deterministic skeleton, are $e_{s1} = (1.4049, 2.8081)$ and $e_{s2} = (4.9040, 4.0619)$.

    eq1.x <- eqs[1, 1]; eq1.y <- eqs[1, 2] # stable focus
    eq2.x <- eqs[3, 1]; eq2.y <- eqs[3, 2] # stable node

Determining the computational domain

Next, the boundaries of the computational domain need to be entered. The ordered-upwind method terminates when the solved area encounters a boundary of this domain, so it's important to choose boundaries carefully. For example, if a stable equilibrium lies on one of the coordinate axes, one should not use that axis as a boundary because the algorithm will immediately terminate. Instead, one should add padding space. For this example, a good choice of boundaries is $x = -0.5-20$ and $y = -0.5-20$. This choice of domain was obtained by examining both the direction field and stochastic realizations.

    bounds.x <- c(-0.5, 20.0)
    bounds.y <- c(-0.5, 20.0)

Determining the mesh size

Finally, the mesh size for the discretization of the domain needs to be specified. In general, the best choice of mesh size will be a compromise between resolution and computational time. The mesh size must be fine enough to precisely track how information moves outward along characteristics from the initial point, but too fine of a mesh size can lead to very long computational times.

    step.number.x <- 1000
    step.number.y <- 1000

Calculating local quasi-potentials

For each stable equilibrium, calculate the local quasi-potential:

    eq1.local <- QPotential(x.rhs = parms.eqn.x, x.start = eq1.x, x.bound = bounds.x, x.num.steps = step.number.x, y.rhs = parms.eqn.y, y.start = eq1.y,  y.bound = bounds.y, y.num.steps = step.number.y)

    eq2.local <- QPotential(x.rhs = parms.eqn.x, x.start = eq2.x, x.bound = bounds.x, x.num.steps = step.number.x, y.rhs = parms.eqn.y, y.start = eq2.y, y.bound = bounds.y, y.num.steps = step.number.y)

You did it! QPotential returns a matrix that has step.number.x by step.number.y rows and columns. We will visualize these local quasi-potential surfaces later in vignette [5. Global quasi-potential visualization].

bmarkslash7/QPot documentation built on Jan. 11, 2020, 11:11 a.m.