A linear dynamical system:
$$\dot{x} = Ax + c(t)$$Has the exact solution:
$$x(t) = e^{tA}x(0) + \int_0^t e^{(t-\tau)A}c(\tau) d\tau$$Over a small time interval, $\Delta t$, we can approximate the input, $c(t)$, as a constant. Pulling this term out of the integral we get:
$$x(t + \Delta t) = e^{\Delta t A}x(t) + c(t) \int_0^{\Delta t} e^{(\Delta t - \tau)A} d\tau$$Which evaluates to:
$$x(t + \Delta t) = e^{\Delta t A}x(t) + A^{-1} (e^{\Delta t A} - 1) c(t) $$To simplify notation, we define $\Phi = e^{\Delta t A}$, which exactly describes how the state evolves over a time step $\Delta t$ when the input is zero, $c(t) = 0$.
$$x(t + \Delta t) = \Phi x(t) + A^{-1} (\Phi - 1) c(t)$$
In [1]:
include("integration.jl")
Out[1]:
In [2]:
using PyPlot
In [3]:
# Target linear system behavior
tmax = 1200.0
dt = 0.1
t = [0:dt:tmax]
# Dynamics
nx = 1 # number of states
A = 1e-8*eye(nx) # perfect integrator
x0 = zeros(nx)
ct = [zeros(2000); ones(3000); zeros(2000); -2*ones(3000); zeros(2001)]'
# Solve
x = solve_linsys(A,x0,ct,dt,tmax)
# Plot
figure()
plot(t,x')
show()
In [5]:
# Build a network
N = 400 # number of neurons
G = 0.1*vcat(ones(200),-ones(200))' # readout matrix
lam_d = 10.0
lam_v = 20.0
nu = 1e-5
mu = 1e-6
xhat,spk_ind,V = solve_boerlin(A,G,lam_d,lam_v,nu,mu,ct,dt,tmax);
In [6]:
figure()
plot(t,V[1,:]')
show()
In [ ]: