1

I'm trying to come up with a python implementation for the Fitzhugh-Nagumo model. 

V_t = V_xx + V(V - a)(1 - V) - W + I
W_t = eps(beta*V - W)

Using the really basic code for eps = 0.05, a = 0.2, beta = 5, I = .1 I can numerically solve the system(with out the V_xx), but I can't seem to figure out how to implement the spacial diffusion.

def func_v(v, w):
    return v * (1 - v) * (v - .2) - w + .1


def func_w(v, w):
    return .05 * (5 * v - w)


def get_yn(t0, v, w, h, t):
    while t0 < t:
        w += h * func_w(v, w)
        v += h * func_v(v, w)
        t0 += h
    return v, w

I know the centered difference formula for second order derivatives is

V_xx(x_i, t) = (V(x_i+1, t) - 2*V(x_i, t) + V(x_i-1, t)) / dx^2

but how would I implement the different values for x_i(let's say from x=0 to 10) in order to get the wave to propagate along the x-axis?

enter image description here

The results should give a wave that propagates something like this.

vlovero
  • 190
  • 1
  • 11
  • Do you have any test data or test samples we can use? – artemis Dec 16 '19 at 20:53
  • @wundermahn what do you mean test data? Like the time interval and initial conditions I use to solve it with out the `V_xx`? – vlovero Dec 16 '19 at 20:57
  • What do you think of previous questions on the topic like https://stackoverflow.com/q/14915398/3088138 or https://stackoverflow.com/q/38116507/3088138 – Lutz Lehmann Dec 17 '19 at 10:14
  • @Dr.LutzLehmann I saw the question in MATLAB and when I translated the code into python, it didn't work for me. – vlovero Dec 17 '19 at 12:42

1 Answers1

1

An ODE solver (any computer program really) only can treat problems that have a finite dimensional state. The state in this PDE is a pair of functions v,w of x. These can not, in the necessary generality, be represented in a computer. Thus you need to work with finite approximations. A first one that is deemed sufficient in many contexts is the simple function table. Then the x derivatives are computed using finite difference formulas.

x = np.linspace(0,L,N+1);
dx = x[1]-x[0];
v0,w0 = initial_functions(x);

def func_v(v, w):
    d2v = -2*v;
    d2v[0] += v[-1]+v[1];
    d2v[1:-1] += v[:-2] + v[2:]
    d2v[-1] += v[-2]+v[0];
    return d2v/dx**2 + v * (1 - v) * (v - .2) - w + .1

etc.

For a proof-of-concept the Euler method may be sufficient, but the values obtained will be questionable. Use a higher order method to get usable results without employing ridiculously small time steps.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • To use the new `func_v` with the initial `v0` and `w0`, would I iterate like `v0 += func_v(v0, w0) ; w0 += func_w(v0, w0)`? Because when I used it, the new plot just got translated up the v-axis instead of traversing the x-axis. – vlovero Dec 17 '19 at 16:45
  • You can call your unchanged Euler method `vn, wn = get_yn(t0, v0, w0, h, t)`. Then use `plot(x,vn,x,wn)`. What you get depends on the initial condition (for instance, a loop around the circle of radius 4, `x=linspace(0,2*pi,N+1)[:-1]; v0=4*np.cos(x); w0=4*sin(x)` or similar, and most critically, on the step size `h`. You will need some stabilization time where the solution converges towards the limit cycle of the VanderPol equation. – Lutz Lehmann Dec 17 '19 at 16:53
  • I'm not sure why it's still not working. Could it be because my `get_yn` was originally supposed to take floats instead of arrays? Whats happening now is I'm getting something similar to what I want, but it's shifted left instead of to the right. – vlovero Dec 17 '19 at 17:37
  • What do you mean with "not working"? If I read right, then the code runs but produces unexpected results? Note that the ODE for a traveling wave solution has order 2+1, so that its first order system has dimension 3. The gif-animation shows a homoclinic orbit in that 3D space, which is quite a special situation, I'm not sure if it is stable under perturbation. It is easy to determine the location of the stationary point, but everything after that is a mathematical problem, off-topic here. – Lutz Lehmann Dec 19 '19 at 00:25
  • The [plots](https://imgur.com/KUxSG9y) are "diffusing" but instead of moving along the x-axis, they're behaving more like what you would get from the head equation instead of a traveling wave. – vlovero Dec 19 '19 at 01:09
  • Yes, that's what I'm seeing too using odeint. You will need a carefully constructed initial condition and the right parameters. Also check the coefficient for the dissipation term, too strong mixing will equalize the solution over space. – Lutz Lehmann Dec 19 '19 at 07:15