laplace eqn errors for different dx

fem_errors/fem_errors.py

+import dolfin as df
+import matplotlib.pyplot as plt
+import numpy as np
+
+df.set_log_level(df.LogLevel.ERROR)
+df.parameters['form_compiler']['optimize'] = True
+
+def laplace(Nx, L, cL, cR):
+    # mesh, function space, function, test function
+    mesh = df.IntervalMesh(Nx, 0, L)
+    x = mesh.coordinates()[:,0]
+    function_space = df.FunctionSpace(mesh, 'P', 1)
+    c, tc = df.Function(function_space), df.TestFunction(function_space)
+
+    # Dirichlet boundary
+    dbc = df.DirichletBC(function_space, df.Constant(0), 'on_boundary')
+
+    def left(x, on_boundary):                                 # inhomogeneous dirichlet boundary
+        return df.near(x[0], 0.0) and on_boundary
+    def right(x, on_boundary):
+        return df.near(x[0], L) and on_boundary
+
+    dbc_left = df.DirichletBC(function_space, df.Constant(cL), left)
+    dbc_right = df.DirichletBC(function_space, df.Constant(cR), right)
+    dbc = [dbc_left, dbc_right]
+
+    # form
+    form = df.inner(c.dx(0), tc.dx(0)) * df.dx(mesh)
+
+    # solve
+    df.solve(form == 0, c, dbc)
+    c_numerical = c.compute_vertex_values(mesh)
+    # exact solution
+    c_exact = ((cR - cL)/L) * x + cL
+
+    # plot
+    # plt.plot(x, c_numerical, 'bo', label='Numerical solution')
+    # plt.plot(x, c_exact, label='Exact solution')
+    # plt.xlabel('x')
+    # plt.ylabel('$c(x)$')
+    # plt.legend()
+    # plt.show()
+
+    error = np.max(np.abs(c_numerical - c_exact))
+    return(error)
+
+# parameters
+L, cL, cR = 2*np.pi, 1, -1
+Nx_array = np.array([5, 10, 20, 40, 80, 100, 200, 400, 800, 1600])
+dx_array = L/Nx_array
+error_array = np.zeros(len(dx_array))
+
+for i in range(0, len(dx_array)):
+    error_array[i] = laplace(Nx_array[i], L, cL, cR)
+
+plt.plot(dx_array, error_array, 'bo')
+plt.show()
\ No newline at end of file