Note: This tutorial was generated from an IPython notebook that can be downloaded here.

Error calculations

FEMpy provides the ability to compute the \(L^\infty\) and \(L^2\) norm errors as well as the \(H^1\) semi-norm error.

This tutorial was made with the following version of FEMpy:

import FEMpy
FEMpy.__version__

'1.0'

Let us examine the error of

\[-\frac{{\rm d}}{{\rm d} x} \left( e^x \frac{{\rm d}}{{\rm d} x} u(x)\right) = -e^x [\cos x -2\sin x - x \cos x - x \sin x]; x \in [0, 1]\]

\[u(0) = 0, u'(1) = \cos(1) - \sin(1)\]

as we vary the mesh step size \(h\).

import numpy as np


def coefficient_function(x):
    return  np.exp(x)


def source_function(x):
    return -np.exp(x) * (np.cos(x) - 2*np.sin(x) - x*np.cos(x) - x*np.sin(x))


def dirichlet_function(x):
    if x == 0:
        return 0


def neumann_function(x):
    if x == 1:
        return np.cos(1) - np.sin(1)

We will need the analytical solution to our problem for the \(L^\infty\) and \(L^2\) norm errors and the derivative of the solution for the \(H^1\) semi-norm.

def analytical_sol(x):
    return x * np.cos(x)


def dx_analytical_sol(x):
    return np.cos(x) - x * np.sin(x)

We will vary our mesh size for \(h \in \left\{ \frac{1}{4}, \frac{1}{8}, \frac{1}{16}, \frac{1}{32}, \frac{1}{128}, \frac{1}{256} \right\}.\)

h_list = [1/(2**n) for n in np.arange(2, 9)]

For our case we will use quadratic finite elements

basis = FEMpy.IntervalBasis1D('quadratic')

Now we can iterate through our mesh sizes and store our errors.

L_inf_err = []
L2_err = []
H1_err = []
for h in h_list:
    mesh = FEMpy.Interval1D(0, 1, h, 'quadratic')
    bcs = FEMpy.BoundaryConditions(mesh, ('dirichlet', 'neumann'), dirichlet_fun=dirichlet_function, neumann_fun=neumann_function, coeff_fun=coefficient_function)

    poisson_eq = FEMpy.Poisson1D(mesh, fe_trial_basis=basis, fe_test_basis=basis, boundary_conditions=bcs)
    poisson_eq.solve(coefficient_function, source_function)

    L_inf_err.append(poisson_eq.l_inf_error(analytical_sol))
    L2_err.append(poisson_eq.l2_error(analytical_sol))
    H1_err.append(poisson_eq.h1_seminorm_error(dx_analytical_sol))

To display our results we can use a pandas dataframe or an astropy table.

from astropy.table import Table

error_table = Table([h_list, L_inf_err, L2_err, H1_err], names=['h', 'L_inf Norm Error', 'L_2 Norm Error', 'H_1 Semi-norm Error'],)
error_table['h'].format = '.4f'; error_table['L_inf Norm Error'].format = '.4e'; error_table['L_2 Norm Error'].format = '.4e'; error_table['H_1 Semi-norm Error'].format = '.4e'
error_table

Table length=7

h	L_inf Norm Error	L_2 Norm Error	H_1 Semi-norm Error
float64	float64	float64	float64
0.2500	3.3279e-04	2.1050e-04	5.4213e-03
0.1250	3.9288e-05	2.6147e-05	1.3534e-03
0.0625	4.7533e-06	3.2632e-06	3.3823e-04
0.0312	5.8395e-07	4.0774e-07	8.4550e-05
0.0156	7.2344e-08	5.0962e-08	2.1137e-05
0.0078	9.0022e-09	6.3701e-09	5.2842e-06
0.0039	1.1227e-09	7.9626e-10	1.3211e-06