check with error in solution

euler_error/diffusion_fixed_domain_euler_errors.py

@@ -25,7 +25,7 @@ def advection_diffusion(Nx, L, Nt, tmax, D):
 
     # initial condition
     c0 = df.Function(SFS)
-    c0.interpolate(df.Expression('1 + 0.1 * cos(2*pi*x[0]/L)', pi=np.pi,L=L, degree=1))
+    c0.interpolate(df.Expression('1 + 0.1 * cos(2*pi*x[0]/L) + 0.2 * cos(3*pi*x[0]/L)', pi=np.pi,L=L, degree=1))
 
     # arrays
     c_array = np.zeros((Nt+1, Nx+1))
@@ -44,8 +44,8 @@ def advection_diffusion(Nx, L, Nt, tmax, D):
     return c_array
 
 # parameters
-Nx, L, D, tmax = 10, 1, 1, 1
-nt_array = np.array([50, 100, 200, 400, 800, 1600])
+Nx, L, D, tmax = 32, 1, 0.1, 0.25
+nt_array = np.array([50, 100, 200])
 dt_array = tmax/nt_array
 mesh = df.IntervalMesh(Nx, 0, L)
 x = mesh.coordinates()[:, 0]
@@ -59,10 +59,11 @@ for i in range(0, len(nt_array)):
     c_exact = np.zeros((nt_array[i]+1, Nx+1))
     times = np.linspace(0, tmax, nt_array[i]+1)
     for j in range(nt_array[i]+1):
-        c_exact[j] = 1 + 0.1 * np.cos(2*np.pi*x/L) * np.exp(-4*np.pi**2*D*times[j]/L**2)
+        c_exact[j] = 1 + 0.1 * np.cos(2*np.pi*x/L) * np.exp(-4*np.pi**2*D*times[j]/L**2) + 0.2 * np.cos(3*np.pi*x/L) * np.exp(-9*np.pi**2*D*times[j]/L**2)
     c = advection_diffusion(Nx, L, nt_array[i], tmax, D)
-    error[i] = np.max(np.abs(c - c_exact))
+    error[i] = np.max(np.abs(c[-1] - c_exact[-1]))
 
+print('Max error for dt=0.0005 at t=0.25 is',error[0])
 popt, pcov = curve_fit(linear_func, np.log(dt_array), np.log(error))
 print('The slope of the graph of log(error) vs log(dt) is', popt[0])
 fig, ax = plt.subplots(1)
@@ -70,5 +71,4 @@ fig, ax = plt.subplots(1)
 ax.loglog(dt_array, error, 'bo')
 ax.set_xlabel('log(dt)')
 ax.set_ylabel('log(error)')
-ax.legend()
 plt.show()