log log plot and calculation of slope

euler_error/diffusion_fixed_domain_euler_errors.py

@@ -3,6 +3,10 @@ import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib.widgets import Slider
 import progressbar
+from scipy.optimize import curve_fit
+
+def linear_func(x, a ,b):
+    return a*x + b
 
 df.set_log_level(df.LogLevel.ERROR)
 df.parameters['form_compiler']['optimize'] = True
@@ -40,8 +44,8 @@ def advection_diffusion(Nx, L, Nt, tmax, D):
     return c_array
 
 # parameters
-Nx, L, D, tmax = 64, 1.0, 0.1, 1.0
-nt_array = np.array([50, 100, 200, 400, 600, 800, 1600])
+Nx, L, D, tmax = 10, 1, 1, 1
+nt_array = np.array([50, 100, 200, 400, 800, 1600])
 dt_array = tmax/nt_array
 mesh = df.IntervalMesh(Nx, 0, L)
 x = mesh.coordinates()[:, 0]
@@ -59,9 +63,12 @@ for i in range(0, len(nt_array)):
     c = advection_diffusion(Nx, L, nt_array[i], tmax, D)
     error[i] = np.max(np.abs(c - c_exact))
 
+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)
-ax.plot(dt_array, error, 'bo')
-ax.set_xlabel('$dt$')
-ax.set_ylabel('error')
-plt.show()
 
+ax.loglog(dt_array, error, 'bo')
+ax.set_xlabel('log(dt)')
+ax.set_ylabel('log(error)')
+ax.legend()
+plt.show()