solution plot, error plot for fixed dt

moving_domain/moving_heat_equation_analytical.py

@@ -54,30 +54,45 @@ x = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[1]
 # exact solution
 c_exact = np.zeros((Nt+1, Nx+1))
 times = np.linspace(0, tmax, Nt+1)
+
 for j in range(Nt+1):
+    l = L * np.exp(alpha * times[j])
     if alpha == 0:
-        c_exact[j] = 1 + 0.2 * np.cos(np.pi * x[j]/L) * np.exp(-np.pi**2 * D * 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:
-        c_exact[j] = 1 + 0.2 * np.cos(np.pi*x[j]*np.exp(-alpha*times[j])/L) * np.exp(-np.pi**2 * D * (1 - np.exp(-2 * alpha * times[j]))/(2 * alpha * L**2)) * np.exp(-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 = advection_diffusion(Nx, L, Nt, tmax, D, alpha)[0]
 times = np.linspace(0, tmax, Nt+1)
 
-fig, ax = plt.subplots(1,1,figsize=(8,6))
+fig, (ax,ax1) = plt.subplots(2,1,figsize=(8,6))
 ax.set_xlabel(r'$x$')
 ax.set_ylabel(r'$c(x,t)$')
-ax.set_xlim(np.min(x)-2, np.max(x)+2)
-ax.set_ylim(np.min(c)-2, np.max(c)+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)
+
+cplot, = ax.plot(x[0], c[0], 'go', mfc='none', label = 'Numerical solution')
+cexactplot, = ax.plot(x[0], c_exact[0], label = 'Exact solution') 
 
-cplot, = ax.plot(x[0], c[0], 'go', ms=1)
-cexactplot, = ax.plot(x[0], c_exact[0]) 
+error = np.abs(c - c_exact)
+print('dt=%6.3f, max error = %6.5f' % (tmax/Nt, np.max(error)))
+errorplot, = ax1.plot(x[0], error[0], 'bo', mfc = 'none')
+ax1.set_ylabel(r'$|c(x,t) - c_{exact}(x,t)|$')
+ax1.set_xlabel('$x$')
+ax1.set_xlim(np.min(x), np.max(x)+2)
 
+ax1.set_ylim([np.min(error)-1, np.max(error)+1])
 def update(value):
         ti = np.abs(times-value).argmin()
         cplot.set_xdata(x[ti])
         cplot.set_ydata(c[ti])
         cexactplot.set_xdata(x[ti])
         cexactplot.set_ydata(c_exact[ti])
+        errorplot.set_xdata(x[ti])
+
+        errorplot.set_ydata(error[ti])
         plt.draw()
         
 sax = plt.axes([0.1, 0.92, 0.7, 0.02])
@@ -86,5 +101,6 @@ slider = Slider(sax, r'$t/\tau$', min(times), max(times),
         fc='#999999')
 slider.drawon = False
 slider.on_changed(update)
+ax.legend(loc=0)
 plt.show()