.. !split
.. _fem:deq:1D:code:global:
Implementation (4)
===================
It is tempting to create a
program with symbolic calculations to perform all the steps in the
computational machinery,
both for automating the work and for documenting the complete algorithms.
As we have seen, there are quite many details involved with
finite element computations and incorporation of boundary conditions.
An implementation will also act as a structured summary of all these details.
Global basis functions
----------------------
We first consider implementations when :math:`{\psi}_i` are global functions
are hence different from zero on most of :math:`\Omega =[0,L]` so all integrals
need integration over the entire domain. Since the expressions for
the entries in the linear system depend on the differential equation
problem being solved, the user must supply the necessary formulas via
Python functions. The implementations here attempt to perform symbolic
calculations, but fall back on numerical computations if the symbolic
ones fail.
The user must prepare a function
``integrand_lhs(psi, i, j)`` for returning the integrand of the
integral that contributes to matrix entry :math:`(i,j)`.
The ``psi`` variable is a Python dictionary holding the basis
functions and their derivatives in symbolic form. More precisely,
``psi[q]`` is a list of
.. math::
\{\frac{d^q{\psi}_0}{dx^q},\ldots,\frac{d^q{\psi}_N}{dx^q}\}
{\thinspace .}
Similarly, ``integrand_rhs(psi, i)`` returns the integrand
for entry number :math:`i` in the right-hand side vector.
Since we also have contributions to the right-hand side vector,
and potentially also the
matrix, from boundary terms without any integral, we introduce two
additional functions, ``boundary_lhs(psi, i, j)`` and
``boundary_rhs(psi, i)`` for returning terms in the variational
formulation that are not to be integrated over the domain :math:`\Omega`.
Examples shown later will explain in more detail how these
user-supplied function may look like.
The linear system can be computed and solved symbolically by
the following function:
.. code-block:: python
import sympy as sp
def solve(integrand_lhs, integrand_rhs, psi, Omega,
boundary_lhs=None, boundary_rhs=None):
N = len(psi[0]) - 1
A = sp.zeros((N+1, N+1))
b = sp.zeros((N+1, 1))
x = sp.Symbol('x')
for i in range(N+1):
for j in range(i, N+1):
integrand = integrand_lhs(psi, i, j)
I = sp.integrate(integrand, (x, Omega[0], Omega[1]))
if boundary_lhs is not None:
I += boundary_lhs(psi, i, j)
A[i,j] = A[j,i] = I # assume symmetry
integrand = integrand_rhs(psi, i)
I = sp.integrate(integrand, (x, Omega[0], Omega[1]))
if boundary_rhs is not None:
I += boundary_rhs(psi, i)
b[i,0] = I
c = A.LUsolve(b)
u = sum(c[i,0]*psi[0][i] for i in range(len(psi[0])))
return u
Not surprisingly, symbolic solution of differential
equations, discretized by a Galerkin or least squares method
with global basis functions,
is of limited interest beyond the simplest problems, because
symbolic integration might be very time consuming or impossible, not
only in ``sympy`` but also in
`WolframAlpha <http://wolframalpha.com>`_
(which applies the perhaps most powerful symbolic integration
software available today: Mathematica). Numerical integration
as an option is therefore desirable.
The extended ``solve`` function below tries to combine symbolic and
numerical integration. The latter can be enforced by the user, or it
can be invoked after a non-successful symbolic integration (being
detected by an ``Integral`` object as the result of the integration
in ``sympy``).
.. see the section :ref:`fem:approx:global:Lagrange`).
Note that for a
numerical integration, symbolic expressions must be converted to
Python functions (using ``lambdify``), and the expressions cannot contain
other symbols than ``x``. The real ``solve`` routine in the
`varform1D.py <http://tinyurl.com/jvzzcfn/fem/varform1D.py>`_
file has error checking and meaningful error messages in such cases.
The ``solve`` code below is a condensed version of the real one, with
the purpose of showing how to automate the Galerkin or least squares
method for solving differential equations in 1D with global basis functions:
.. code-block:: python
def solve(integrand_lhs, integrand_rhs, psi, Omega,
boundary_lhs=None, boundary_rhs=None, symbolic=True):
N = len(psi[0]) - 1
A = sp.zeros((N+1, N+1))
b = sp.zeros((N+1, 1))
x = sp.Symbol('x')
for i in range(N+1):
for j in range(i, N+1):
integrand = integrand_lhs(psi, i, j)
if symbolic:
I = sp.integrate(integrand, (x, Omega[0], Omega[1]))
if isinstance(I, sp.Integral):
symbolic = False # force num.int. hereafter
if not symbolic:
integrand = sp.lambdify([x], integrand)
I = sp.mpmath.quad(integrand, [Omega[0], Omega[1]])
if boundary_lhs is not None:
I += boundary_lhs(psi, i, j)
A[i,j] = A[j,i] = I
integrand = integrand_rhs(psi, i)
if symbolic:
I = sp.integrate(integrand, (x, Omega[0], Omega[1]))
if isinstance(I, sp.Integral):
symbolic = False
if not symbolic:
integrand = sp.lambdify([x], integrand)
I = sp.mpmath.quad(integrand, [Omega[0], Omega[1]])
if boundary_rhs is not None:
I += boundary_rhs(psi, i)
b[i,0] = I
c = A.LUsolve(b)
u = sum(c[i,0]*psi[0][i] for i in range(len(psi[0])))
return u
Example: constant right-hand side
---------------------------------
To demonstrate the code above, we address
.. math::
-u''(x)=b,\quad x\in\Omega=[0,1],\quad u(0)=1,\ u(1)=0,
with :math:`b` as a (symbolic) constant. A possible basis for the space :math:`V`
is :math:`{\psi}_i(x) = x^{i+1}(1-x)`, :math:`i\in{\mathcal{I}_s}`. Note that
:math:`{\psi}_i(0)={\psi}_i(1)=0` as required by the Dirichlet conditions.
We need a :math:`B(x)` function to take care of the known boundary
values of :math:`u`. Any function :math:`B(x)=1-x^p`, :math:`p\in\mathbb{R}`, is a candidate,
and one arbitrary choice from this family
is :math:`B(x)=1-x^3`. The unknown function is then written as
.. math::
u(x) = B(x) + \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j(x){\thinspace .}
Let us use the Galerkin method to derive the variational formulation.
Multiplying the differential
equation by :math:`v` and integrate by parts yield
.. math::
\int_0^1 u'v' {\, \mathrm{d}x} = \int_0^1 fv {\, \mathrm{d}x}\quad\forall v\in V,
and with :math:`u=B + \sum_jc_j{\psi}_j` we get the linear system
.. math::
\sum_{j\in{\mathcal{I}_s}}\left(\int_0^1{\psi}_i'{\psi}_j' {\, \mathrm{d}x}\right)c_j = \int_0^1(f-B'){\psi}_i {\, \mathrm{d}x},
\quad i\in{\mathcal{I}_s}{\thinspace .}
The application can be coded as follows in ``sympy``:
.. code-block:: python
x, b = sp.symbols('x b')
f = b
B = 1 - x**3
dBdx = sp.diff(B, x)
# Compute basis functions and their derivatives
N = 3
psi = {0: [x**(i+1)*(1-x) for i in range(N+1)]}
psi[1] = [sp.diff(psi_i, x) for psi_i in psi[0]]
def integrand_lhs(psi, i, j):
return psi[1][i]*psi[1][j]
def integrand_rhs(psi, i):
return f*psi[0][i] - dBdx*psi[1][i]
Omega = [0, 1]
u_bar = solve(integrand_lhs, integrand_rhs, psi, Omega,
verbose=True, symbolic=True)
u = B + u_bar
print 'solution u:', sp.simplify(sp.expand(u))
The printout of ``u`` reads ``-b*x**2/2 + b*x/2 - x + 1``.
Note that expanding ``u`` and then simplifying is in the present case
necessary to get a compact, final expression with ``sympy``.
A non-expanded ``u`` might be preferable in other cases - this depends on
the problem in question.
The exact solution :math:`{u_{\small\mbox{e}}}(x)` can be derived by
some ``sympy`` code that closely follows the examples in
the section :ref:`fem:deq:1D:models:simple`. The idea is to integrate
:math:`-u''=b` twice and determine the integration constants from
the boundary conditions:
.. code-block:: python
C1, C2 = sp.symbols('C1 C2') # integration constants
f1 = sp.integrate(f, x) + C1
f2 = sp.integrate(f1, x) + C2
# Find C1 and C2 from the boundary conditions u(0)=0, u(1)=1
s = sp.solve([u_e.subs(x,0) - 1, u_e.subs(x,1) - 0], [C1, C2])
# Form the exact solution
u_e = -f2 + s[C1]*x + s[C2]
print 'analytical solution:', u_e
print 'error:', sp.simplify(sp.expand(u - u_e))
The last line prints ``0``, which is not surprising when
:math:`{u_{\small\mbox{e}}}(x)` is a parabola and our approximate :math:`u` contains polynomials up to
degree 4. It suffices to have :math:`N=1`, i.e., polynomials of degree
2, to recover the exact solution.
We can play around with the code and test that with :math:`f\sim x^p`,
the solution is a polynomial of degree :math:`p+2`, and :math:`N=p+1` guarantees
that the approximate solution is exact.
Although the symbolic code is capable of integrating many choices of :math:`f(x)`,
the symbolic expressions for :math:`u` quickly become lengthy and non-informative,
so numerical integration in the code, and hence numerical answers,
have the greatest application potential.
Finite elements
---------------
Implementation of the finite element algorithms for differential
equations follows closely the algorithm for approximation of functions.
The new additional ingredients are
1. other types of integrands (as implied by the variational formulation)
2. additional boundary terms in the variational formulation for
Neumann boundary conditions
3. modification of element matrices and vectors due to Dirichlet
boundary conditions
Point 1 and 2 can be taken care of by letting the user supply
functions defining the integrands and boundary terms on the
left- and right-hand side of the equation system:
.. code-block:: python
integrand_lhs(phi, r, s, x)
boundary_lhs(phi, r, s, x)
integrand_rhs(phi, r, x)
boundary_rhs(phi, r, x)
Here, ``phi`` is a dictionary where ``phi[q]`` holds a list of
the derivatives of order ``q`` of the basis functions at the
an evaluation point; ``r`` and ``s`` are indices for the corresponding
entries in the element matrix and vector, and ``x`` is the global
coordinate value corresponding to the current evaluation point.
Given a mesh represented by ``vertices``, ``cells``, and ``dof_map`` as
explained before, we can write a pseudo Python code to list all
the steps in the computational algorithm for finite element solution
of a differential equation.
.. code-block:: python
<Declare global matrix and rhs: A, b>
for e in range(len(cells)):
# Compute element matrix and vector
n = len(dof_map[e]) # no of dofs in this element
h = vertices[cells[e][1]] - vertices[cells[e][1]]
<Declare element matrix and vector: A_e, b_e>
# Integrate over the reference cell
points, weights = <numerical integration rule>
for X, w in zip(points, weights):
phi = <basis functions and derivatives at X>
detJ = h/2
x = <affine mapping from X>
for r in range(n):
for s in range(n):
A_e[r,s] += integrand_lhs(phi, r, s, x)*detJ*w
b_e[r] += integrand_rhs(phi, r, x)*detJ*w
# Add boundary terms
for r in range(n):
for s in range(n):
A_e[r,s] += boundary_lhs(phi, r, s, x)*detJ*w
b_e[r] += boundary_rhs(phi, r, x)*detJ*w
# Incorporate essential boundary conditions
for r in range(n):
global_dof = dof_map[e][r]
if global_dof in essbc_dofs:
# dof r is subject to an essential condition
value = essbc_docs[global_dof]
# Symmetric modification
b_e -= value*A_e[:,r]
A_e[r,:] = 0
A_e[:,r] = 0
A_e[r,r] = 1
b_e[r] = value
# Assemble
for r in range(n):
for s in range(n):
A[dof_map[e][r], dof_map[e][r]] += A_e[r,s]
b[dof_map[e][r] += b_e[r]
<solve linear system>