Software 10.6.1 (Equilibria and their linearisations). SymPy finds the equilibria of the competition model and the eigenvalues, trace and determinant of the Jacobian at each:

>>> import sympy as sp
>>> x, y = sp.symbols('x y', real=True)
>>> F = sp.Matrix([x * (3 - x - 2*y), y * (2 - x - y)])   # competition model
>>> eqs = sp.solve(list(F), [x, y], dict=True); eqs
[{x: 0, y: 0}, {x: 0, y: 2}, {x: 1, y: 1}, {x: 3, y: 0}]
>>> J = F.jacobian([x, y])
>>> for e in eqs:
...     Je = J.subs(e)
...     print((e[x], e[y]), list(Je.eigenvals()), Je.trace(), Je.det())
...
(0, 0) [3, 2] 5 6
(0, 2) [-1, -2] -3 2
(1, 1) [-sqrt(2) - 1, -1 + sqrt(2)] -2 -1
(3, 0) [-1, -3] -4 3

Reading the output. The origin has eigenvalues \(3\) and \(2\), an unstable node; \((0, 2)\) and \((3, 0)\) have negative eigenvalues, stable nodes; \((1, 1)\) has \(\Delta = -1 < 0\), a saddle, with eigenvalues \(-1 \pm \sqrt 2\). In the trace–determinant plane of Figure 8, the saddle lies below the axis and the nodes between the axis and the parabola.

MATLAB/Octave: syms x y; S = solve([x*(3-x-2*y) == 0, y*(2-x-y) == 0]) and eig(subs(jacobian(...))).

Software 10.6.2 (The separatrix and the basins). The stable manifold of the saddle is traced by integrating backwards from points \(10^{-4}\) away along the stable eigenvector, stopping when the orbit leaves the box \([0, 4]^2\); four initial states are followed forward:

>>> import numpy as np
>>> from scipy.integrate import solve_ivp
>>> f = lambda t, z: [z[0] * (3 - z[0] - 2*z[1]), z[1] * (2 - z[0] - z[1])]
>>> Js = np.array(J.subs({x: 1, y: 1}), dtype=float)
>>> lam, V = np.linalg.eig(Js); lam
array([ 0.41421356, -2.41421356])
>>> vs = V[:, np.argmin(lam)]              # stable eigenvector of the saddle
>>> back = lambda t, z: [-c for c in f(t, z)]
>>> out = lambda t, z: max(z) - 4; out.terminal = True   # leave the box
>>> arms = [solve_ivp(back, (0, 12), 1 + s * 1e-4 * vs, events=out,
...                   rtol=1e-10, atol=1e-12).y[:, -1] for s in (1, -1)]
>>> np.round(arms, 4)                      # far ends of the separatrix
array([[4.    , 2.7121],
       [0.    , 0.    ]])
>>> for z0 in ([0.5, 0.6], [0.6, 0.5], [2.5, 2.0], [2.0, 2.5]):
...     print(z0, np.round(solve_ivp(f, (0, 40), z0, rtol=1e-9).y[:, -1], 4))
...
[0.5, 0.6] [3. 0.]
[0.6, 0.5] [3. 0.]
[2.5, 2.0] [0. 2.]
[2.0, 2.5] [0. 2.]

Reading the output. The stable eigenvalue is \(-2.414\). One arm of the separatrix runs back to the unstable node at the origin, the other leaves the box at \((4, 2.71)\). The starting points \((0.5, 0.6)\) and \((0.6, 0.5)\), on the same side of it, both end at \((3, 0)\), although in the first the second species starts larger: the separatrix is not the diagonal. The points \((2.5, 2)\) and \((2, 2.5)\), on the other side, end at \((0, 2)\). A small change in the initial state across the separatrix changes the outcome completely.

MATLAB/Octave: integrate @(t, z) -f(t, z) with ode45 and an Events function for the box.

Software 10.6.3 (The conserved quantity of Lotka–Volterra). For \(a = 1\), \(b = 0.5\), \(c = 0.75\), \(d = 0.25\), the value of \(H\) along a numerical solution over \(50\) time units:

>>> a, b, c, d = 1.0, 0.5, 0.75, 0.25      # x' = ax - bxy, y' = -cy + dxy
>>> lv = lambda t, z: [a*z[0] - b*z[0]*z[1], -c*z[1] + d*z[0]*z[1]]
>>> H = lambda z: d*z[0] - c*np.log(z[0]) + b*z[1] - a*np.log(z[1])
>>> s = solve_ivp(lv, (0, 50), [2.0, 1.0], rtol=1e-10, atol=1e-12,
...               dense_output=True)
>>> Hs = [H(s.sol(t)) for t in np.linspace(0, 50, 6)]; np.round(Hs, 9)
array([0.48013962, 0.48013962, 0.48013962, 0.48013962, 0.48013962,
       0.48013962])

Reading the output. \(H\) is constant to nine decimals over several periods, as Theorem 10.5.1(i) requires: a check on the accuracy of the integration. A general-purpose integrator does not conserve \(H\) exactly; over very long times its small errors accumulate, and the computed orbit drifts slowly across the level sets. Integrators that preserve such invariants by construction, symplectic methods, are discussed in the last section.

MATLAB/Octave: evaluate d*x - c*log(x) + b*y - a*log(y) along the ode45 solution.