no signaling -- fixed bdry, fft

no_c_fixed_bdry_fft.py

+import numpy as np
+from scipy.integrate import odeint
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+
+def reaction_rho(rho, rho0, krho):
+    return (-krho * (rho - rho0))
+
+def time_derivative(c, t):
+    # split
+    u, rho = np.split(c, 2)
+
+    # compute FFT
+    uq = np.fft.fftshift(np.fft.fft(u))
+    rhoq = np.fft.fftshift(np.fft.fft(rho))
+    vq = (-K * q**2 * uq + 1j * q * lamda * rhoq)/(1 + eta * q**2)
+
+    # RHS in Fourier-space
+    RHS_u_q = vq
+    RHS_rho_q = -rho0 * 1j * q * vq - Drho * q**2 * rhoq
+
+    # RHS in real space
+    RHS_u = np.real(np.fft.ifft(np.fft.ifftshift(RHS_u_q)))
+    RHS_rho = np.real(np.fft.ifft(np.fft.ifftshift(RHS_rho_q))) + reaction_rho(rho, rho0, krho)
+
+    return np.concatenate([RHS_u, RHS_rho])
+
+#parameters
+Lx = 2*np.pi
+Nx = 100
+T = 100
+dt = 0.01
+Nt = int(T/dt)
+times = np.linspace(0,T,Nt)
+
+K = 1
+eta = 1
+lamda = 8
+krho = 1
+rho0 = 1
+Drho = 1
+
+x = np.linspace(0, Lx, Nx)
+q = np.fft.fftshift(np.fft.fftfreq(len(x), d=Lx/(2*np.pi*Nx)))
+
+# initial conditions
+u0 = np.sin(x)
+rho0 = np.cos(x)
+
+c0 = np.concatenate([u0, rho0])
+
+#integrate in time
+c = odeint(time_derivative, c0, times)
+
+#split and reshape solution arrays
+u, rho = np.split(c, 2, axis=1)
+u = u.reshape((Nt, Nx))
+rho = rho.reshape((Nt, Nx))
+
+#plotting
+fig, (axu, axrho) = plt.subplots(2,1, sharex=True, figsize=(8,6))
+axu.set_xlabel(r'$x$')
+axu.set_ylabel(r'$u(x,t)$')
+axrho.set_ylabel(r'$\rho(x,t)$')
+axu.set_ylim(np.min(u), np.max(u))
+axrho.set_ylim(np.min(rho), np.max(rho))
+
+uplot, = axu.plot(x, u[0])
+rhoplot, = axrho.plot(x, rho[0])
+
+def update(value):
+    ti = np.abs(times-value).argmin()
+    uplot.set_ydata(u[ti])
+    rhoplot.set_ydata(rho[ti])
+    plt.draw()
+    
+sax = plt.axes([0.1, 0.92, 0.7, 0.02])
+slider = Slider(sax, r'$t/\tau$', min(times), max(times),
+        valinit=min(times), valfmt='%3.1f',
+        fc='#999999')
+slider.drawon = False
+slider.on_changed(update)
+plt.show()
+