only exponential growth

moving_domain/moving_heat_equation_analytical.py

@@ -48,7 +48,7 @@ def advection_diffusion(Nx, L, Nt, tmax, D, alpha):
     return c_array, x_array
     return c_array, x_array
 
 
 # plot c(x,t) numerical and analytical for given dt
 # plot c(x,t) numerical and analytical for given dt
-dx, L, dt, tmax, D, alpha = 0.001, 1, 0.001, 1, 0.01, 0.1
+dx, L, dt, tmax, D, alpha = 0.01, 1, 0.01, 20, 0.01, 0.1
 Nx = int(L/dx)
 Nx = int(L/dx)
 Nt = int(tmax/dt)
 Nt = int(tmax/dt)
 x = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[1]
 x = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[1]
@@ -58,18 +58,17 @@ c_exact = np.zeros((Nt+1, Nx+1))
 times = np.linspace(0, tmax, Nt+1)
 times = np.linspace(0, tmax, Nt+1)
 
 
 for j in range(Nt+1):
 for j in range(Nt+1):
-    
+    l = L * np.exp(alpha * times[j])
     if alpha == 0:
     if alpha == 0:
         beta = -D * np.pi**2 * times[j]/L**2
         beta = -D * np.pi**2 * times[j]/L**2
-        c_exact[j] = 1 + 0.2 * np.cos(np.pi * x[j]/L) * np.exp(beta)
     else:
     else:
-        l = L * np.exp(alpha * times[j])
         beta = (-D * np.pi**2 /(2 *alpha * L**2)) * (1 - np.exp(-2 * alpha * times[j]))
         beta = (-D * np.pi**2 /(2 *alpha * L**2)) * (1 - np.exp(-2 * alpha * times[j]))
     c_exact[j] = np.exp(-alpha * times[j])* (1 + 0.2 * np.cos(np.pi*x[j]/l) * np.exp(beta))
     c_exact[j] = np.exp(-alpha * times[j])* (1 + 0.2 * np.cos(np.pi*x[j]/l) * np.exp(beta))
 
 
+# numerical solution
 c = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[0]
 c = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[0]
 
 
-fig, (ax,ax1) = plt.subplots(2,1,figsize=(8,6))
+fig, (ax, ax1) = plt.subplots(2,1,figsize=(8,6))
 ax.set_xlabel(r'$x$')
 ax.set_xlabel(r'$x$')
 ax.set_ylabel(r'$c(x,t)$')
 ax.set_ylabel(r'$c(x,t)$')
 ax.set_xlim(np.min(x), np.max(x)+2)
 ax.set_xlim(np.min(x), np.max(x)+2)
@@ -79,13 +78,13 @@ cplot, = ax.plot(x[0], c[0], '--',label = 'Numerical solution')
 cexactplot, = ax.plot(x[0], c_exact[0], label = 'Exact solution') 
 cexactplot, = ax.plot(x[0], c_exact[0], label = 'Exact solution') 
 
 
 error = np.abs(c - c_exact)
 error = np.abs(c - c_exact)
-print('dt=%6.3f, max error = %6.5f' % (tmax/Nt, np.max(error)))
+print('dt = %6.3f, max error = %6.5f' % (tmax/Nt, np.max(error)))
 errorplot, = ax1.plot(x[0], error[0], 'bo', mfc = 'none')
 errorplot, = ax1.plot(x[0], error[0], 'bo', mfc = 'none')
 ax1.set_ylabel(r'$|c(x,t) - c_{exact}(x,t)|$')
 ax1.set_ylabel(r'$|c(x,t) - c_{exact}(x,t)|$')
 ax1.set_xlabel('$x$')
 ax1.set_xlabel('$x$')
 ax1.set_xlim(np.min(x), np.max(x)+2)
 ax1.set_xlim(np.min(x), np.max(x)+2)
-
 ax1.set_ylim([np.min(error)-1, np.max(error)+1])
 ax1.set_ylim([np.min(error)-1, np.max(error)+1])
+
 def update(value):
 def update(value):
         ti = np.abs(times-value).argmin()
         ti = np.abs(times-value).argmin()
         cplot.set_xdata(x[ti])
         cplot.set_xdata(x[ti])
@@ -104,5 +103,13 @@ slider = Slider(sax, r'$t/\tau$', min(times), max(times),
 slider.drawon = False
 slider.drawon = False
 slider.on_changed(update)
 slider.on_changed(update)
 ax.legend(loc=0)
 ax.legend(loc=0)
+
+fig, ax2 = plt.subplots(1,1,figsize=(8,6))
+ax2.semilogy(times, c[:,0], label='$c_{num}(0,t)$')
+if alpha == 0:
+    ax2.semilogy(times, 1 + 0.2 * np.exp(-D *np.pi**2 * times/L**2), label = '$c_{an}(0,t)$')
+else:
+    ax2.semilogy(times, np.exp(-alpha * times) * (1 + 0.2 * np.exp(-D * (np.pi**2/(2*alpha*L**2)) * (1 - np.exp(-2 * alpha * times)))), label = '$c_{an}(0,t)$' )
+ax2.legend()
 plt.show()
 plt.show()