Lotka-Volterra automatically differentiated using tensorflow.

lv.py

import numpy as np
import tensorflow as tf


def LV(N1, N2, eps1, eps2, gam1, gam2):
    dt = tf.constant(1.)
    states = [(N1, N2)]
    for i in range(1, 13):
        states.append((states[i-1][0] + (states[i-1][0] * (eps1-gam1*states[i-1][1])) * dt, states[i-1][1] - states[i-1][1] * (eps2-gam2*states[i-1][0])) * dt)
    return states[-1]

with tf.GradientTape(persistent=True) as tape:
    (N1, N2, eps1, eps2, gam1, gam2) = arg = [tf.Variable(x) for x in [120., 60., 7e-3, 4e-2, 5e-4, 5e-4]]
    (fN1, fN2) = LV(N1, N2, eps1, eps2, gam1, gam2)

print(tape.jacobian(fN1, arg))
print(tape.jacobian(fN2, arg))

# [<tf.Tensor: shape=(), dtype=float32, numpy=0.6423073>, <tf.Tensor: shape=(), dtype=float32, numpy=-0.570499>, <tf.Tensor: shape=(), dtype=float32, numpy=1046.2244>, <tf.Tensor: shape=(), dtype=float32, numpy=194.50531>, <tf.Tensor: shape=(), dtype=float32, numpy=-68459.875>, <tf.Tensor: shape=(), dtype=float32, numpy=-21517.64>]
# [<tf.Tensor: shape=(), dtype=float32, numpy=0.34937298>, <tf.Tensor: shape=(), dtype=float32, numpy=1.0397532>, <tf.Tensor: shape=(), dtype=float32, numpy=228.17729>, <tf.Tensor: shape=(), dtype=float32, numpy=-801.1748>, <tf.Tensor: shape=(), dtype=float32, numpy=-14474.355>, <tf.Tensor: shape=(), dtype=float32, numpy=83849.51>]