checked for errors, all clear

no_c_fixed_bdry_fft.py

@@ -3,7 +3,26 @@ from scipy.integrate import odeint
 import matplotlib.pyplot as plt
 import matplotlib.pyplot as plt
 from matplotlib.widgets import Slider
 from matplotlib.widgets import Slider
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 from mpl_toolkits.axes_grid1 import make_axes_locatable
+import time
+import progressbar 
 
 
+#parameters
+Lx = 2*np.pi
+Nx = 200
+T = 100
+dt = 0.1
+Nt = int(T/dt)
+times = np.linspace(0,T,Nt)
+
+K = 1
+eta = 1
+lamda = 1
+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)))
 
 
 def reaction_rho(rho, rho0, krho):
 def reaction_rho(rho, rho0, krho):
     return (-krho * (rho - rho0))
     return (-krho * (rho - rho0))
@@ -14,7 +33,7 @@ def time_derivative(c, t):
 
 
     # compute FFT
     # compute FFT
     uq = np.fft.fftshift(np.fft.fft(u))
     uq = np.fft.fftshift(np.fft.fft(u))
-    rhoq = np.fft.fftshift(np.fft.fft(rho))
+    rhoq = np.fft.fftshift(np.fft.fft(rho - rho0))
     vq = (-K * q**2 * uq + 1j * q * lamda * rhoq)/(1 + eta * q**2)
     vq = (-K * q**2 * uq + 1j * q * lamda * rhoq)/(1 + eta * q**2)
 
 
     # RHS in Fourier-space
     # RHS in Fourier-space
@@ -27,37 +46,19 @@ def time_derivative(c, t):
 
 
     return np.concatenate([RHS_u, RHS_rho])
     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
 # initial conditions
-u0 = np.sin(x)
-rho0 = np.cos(x)
+u00 = np.sin(x)
+# u00 = u00 * np.ones_like(x) + 0.01 * np.random.randn(x.shape[0])
+rho00 = np.cos(x) + np.cos(x/2)
+# rho00 = rho00 * np.ones_like(x) + 0.01 * np.random.randn(x.shape[0])
 
 
-c0 = np.concatenate([u0, rho0])
+c0 = np.concatenate([u00, rho00])
 
 
 #integrate in time
 #integrate in time
 c = odeint(time_derivative, c0, times)
 c = odeint(time_derivative, c0, times)
 
 
 #split and reshape solution arrays
 #split and reshape solution arrays
 u, rho = np.split(c, 2, axis=1)
 u, rho = np.split(c, 2, axis=1)
-u = u.reshape((Nt, Nx))
-rho = rho.reshape((Nt, Nx))
 
 
 #plotting
 #plotting
 fig, (axu, axrho) = plt.subplots(2,1, sharex=True, figsize=(8,6))
 fig, (axu, axrho) = plt.subplots(2,1, sharex=True, figsize=(8,6))