paper parameters

moving_domain/moving_heat_equation_analytical.py

@@ -48,7 +48,9 @@ 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
-Nx, L, Nt, tmax, D, alpha = 64, 1, 100, 1, 1, 1
+dx, L, dt, tmax, D, alpha = 0.001, 1, 0.001, 1, 0.01, 0.1
+Nx = int(L/dx)
+Nt = int(tmax/dt)
 x = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[1]
 x = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[1]
 
 
 # exact solution
 # exact solution
@@ -56,16 +58,16 @@ 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)
         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] = 1 + 0.2 * np.cos(np.pi*x[j]/l) * np.exp(-alpha * times[j]) * np.exp(beta)
+        c_exact[j] = np.exp(-alpha * times[j])* (1 + 0.2 * np.cos(np.pi*x[j]/l) * np.exp(beta))
 
 
 c = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[0]
 c = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[0]
-times = np.linspace(0, tmax, Nt+1)
 
 
 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$')
@@ -73,7 +75,7 @@ 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)
 ax.set_ylim(min(np.min(c), np.min(c_exact))-1, max(np.max(c), np.max(c_exact))+1)
 ax.set_ylim(min(np.min(c), np.min(c_exact))-1, max(np.max(c), np.max(c_exact))+1)
 
 
-cplot, = ax.plot(x[0], c[0], 'go', mfc='none', label = 'Numerical solution')
+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)