.. !split

Exercises
=========

.. --- begin exercise ---

.. _fem:approx:exer:linalg1:

Problem 1: Linear algebra refresher
-----------------------------------

Look up the topic of *vector space* in your favorite linear algebra
book or search for the term at Wikipedia.

**a)**
Prove that vectors in the plane :math:`(a,b)` form a vector space
by showing that all the axioms of a vector space are satisfied.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

**b)**
Prove that all linear functions of the form :math:`ax+b` constitute a vector space,
:math:`a,b\in\mathbb{R}`.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

**c)**
Show that all quadratic functions of the form :math:`1 + ax^2 + bx` *do not*
constitute a vector space.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

**d)**
Check out
the topic of *inner product spaces*. Suggest a possible inner product
for the space of all linear functions of the form :math:`ax+b`, :math:`a,b\in\mathbb{R}`,
defined on some interval :math:`\Omega=[A,B]`.
Show that this particular inner product satisfies the
general requirements of an inner product in a vector space.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``linalg1``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:exer:vec:3Dby2D:

Problem 2: Approximate a three-dimensional vector in a plane
------------------------------------------------------------

Given :math:`\boldsymbol{f} = (1,1,1)` in :math:`\mathbb{R}^3`, find the best approximation vector
:math:`\boldsymbol{u}` in the plane spanned by the unit vectors :math:`(1,0)` and :math:`(0,1)`.
Repeat the calculations using the vectors :math:`(2,1)` and :math:`(1,2)`.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``vec111_approx``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:exer:parabola_sine:

Problem 3: Approximate a parabola by a sine
-------------------------------------------

Given the function :math:`f(x)=1 + 2x(1-x)` on :math:`\Omega=[0,1]`, we want to
find an approximation in the function space

.. math::
         V = \hbox{span}\{1, \sin(\pi x)\}{\thinspace .}

**a)**
Sketch or plot :math:`f(x)`. Think intuitively how an expansion in terms
of the basis functions of :math:`V`, :math:`{\psi}_0(x)=1`, :math:`{\psi}_1(x)=\sin(\pi x)`,
can be construction to yield a best approximation to :math:`f`. Or
phrased differently, see if you can guess the coefficients :math:`c_0` and :math:`c_1`
in the expansion

.. math::
         u(x) = c_0{\psi}_0 + c_1{\psi}_1 = c_0 + c_1\sin (\pi x){\thinspace .}

Compute the :math:`L^2` error :math:`||f-u||_{L^2}=(\int_0^1(f-u)^2{\, \mathrm{d}x})^{1/2}`.

.. --- begin hint in exercise ---

**Hint.**
If you make a mesh function ``e`` of the error
on some mesh with uniformly spaced coordinates in
the array ``xc``, the integral can be approximated as ``np.sqrt(dx*np.sum(e**2))``,
where ``dx=xc[0]-xc[1]`` is the mesh spacing and ``np`` is an alias for the
``numpy`` module in Python.

.. --- end hint in exercise ---

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

**b)**
Perform the hand calculations for a least squares approximation.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``parabola_sin``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:exer:exp:powers:

Problem 4: Approximate the exponential function by power functions
------------------------------------------------------------------

Let :math:`V` be a function space with basis functions :math:`x^i`,
:math:`i=0,1,\ldots,N`.  Find the best approximation to :math:`f(x)=\exp(-x)` on
:math:`\Omega =[0,8]` among all functions in :math:`V` for :math:`N=2,4,6`. Illustrate
the three approximations in three separate plots.

.. --- begin hint in exercise ---

**Hint.**
Apply the ``lest_squares`` and
``comparison_plot`` functions in the ``approx1D.py`` module as these
make the exercise easier to solve.

.. --- end hint in exercise ---

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``exp_powers``.

.. Taylor: these polynomials go so far off on [0,8] that it is not a

.. good idea to add them.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:exer:sin:powers:

Problem 5: Approximate the sine function by power functions
-----------------------------------------------------------

In this exercise we want to approximate the sine function by polynomials
of order :math:`N+1`. Consider two bases:

.. math::
        
        V_1 &= \{x, x^3, x^5, \ldots, x^{N-2}, x^N \}, \\ 
        V_2 &= \{1,x,x^2,x^3,\ldots, x^N\}{\thinspace .}
        

The basis :math:`V_1` is motivated by the fact that the Taylor polynomial
approximation to the sine function has only odd powers, while :math:`V_2`
is motivated by the assumption that also the even powers could
improve the approximation in a least-squares setting.

Compute the best approximation to :math:`f(x)=\sin(x)` among all functions
in :math:`V_1` and :math:`V_2` on two domains of increasing sizes: :math:`\Omega_{1,k} =
[0, k\pi]`, :math:`k=2,3\ldots,6` and :math:`\Omega_{2,k} = [-k\pi /2, k\pi/2]`,
:math:`k=2,3,4,5`.  Make plots for all combinations of :math:`V_1`, :math:`V_2`,
:math:`\Omega_1`, :math:`\Omega_2`, :math:`k=2,3,\ldots,6`.

Add a plot of the :math:`N`-th degree Taylor polynomial approximation of
:math:`\sin(x)` around :math:`x=0`.

.. --- begin hint in exercise ---

**Hint.**
You can make a loop over :math:`V_1` and :math:`V_2`, a loop over
:math:`\Omega_1` and :math:`\Omega_2`, and a loop over :math:`k`. Inside the loops,
call the functions ``least_squares`` and
``comparison_plot`` from the ``approx1D`` module.
:math:`N=7` is a suggested value.

.. --- end hint in exercise ---

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``sin_powers``.

.. Solveig explanations based on f-B approx

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:exer:tanh:sine1:

Problem 6: Approximate a steep function by sines
------------------------------------------------

Find the best approximation of :math:`f(x) = \tanh (s(x-\pi))` on
:math:`[0, 2\pi]` in the space :math:`V` with basis
:math:`{\psi}_i(x) = \sin((2i+1)x)`, :math:`i\in{\mathcal{I}_s} = \{0,\ldots,N\}`.
Make a movie showing how :math:`u=\sum_{j\in{\mathcal{I}_s}}c_j{\psi}_j(x)`
approximates :math:`f(x)` as :math:`N` grows. Choose :math:`s` such that :math:`f` is
steep (:math:`s=20` is appropriate).

.. --- begin hint in exercise ---

**Hint 1.**
One may naively call the ``least_squares_orth`` and ``comparison_plot``
from the ``approx1D`` module in a loop and extend the basis with
one new element in each pass. This approach
implies a lot of recomputations.
A more efficient strategy is to let ``least_squares_orth``
compute with only one basis function at a time and accumulate
the corresponding ``u`` in the total solution.

.. --- end hint in exercise ---

.. --- begin hint in exercise ---

**Hint 2.**
``ffmpeg`` or ``avconv`` may skip frames when plot files are combined to
a movie. Since there are few files and we want to see each of them,
use ``convert`` to make an animated GIF file (``-delay 200`` is
suitable).

.. --- end hint in exercise ---

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``tanh_sines``.

.. Closing remarks for this Problem

Remarks          (1)
~~~~~~~~~~~~~~~~~~~~

Approximation of a discontinuous (or steep) :math:`f(x)` by sines,
results in slow convergence and oscillatory behavior of the
approximation close to the abrupt changes in :math:`f`.
This is known as the `Gibb's phenomenon <http://en.wikipedia.org/wiki/Gibbs_phenomenon>`__.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:exer:tanh:sine3:

Problem 7: Approximate a steep function by sines with boundary adjustment
-------------------------------------------------------------------------

We study the same approximation problem as in
:ref:`fem:approx:exer:tanh:sine1`. Since :math:`{\psi}_i(0)={\psi}_i(2\pi)=0`
for all :math:`i`, :math:`u=0` at the boundary points :math:`x=0` and :math:`x=2\pi`, while
:math:`f(0)=-1` and :math:`f(2\pi)=1`. This discrepancy at the boundary can be
removed by adding a boundary function :math:`B(x)`:

.. math::
        
        u(x) = B(x) + \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j(x),
        

where :math:`B(x)` has the right boundary values: :math:`B(x_L)=f(x_L)` and
:math:`B(x_R)=f(x_R)`, with :math:`x_L=0` and :math:`x_R=2\pi` as the boundary points.
A linear choice of :math:`B(x)` is

.. math::
         B(x) = \frac{(x_R-x)f(x_L) + (x-x_L)f(x_R)}{x_R-x_L}{\thinspace .}

**a)**
Use the basis
:math:`{\psi}_i(x) = \sin((i+1)x)`, :math:`i\in{\mathcal{I}_s} = \{0,\ldots,N\}`
and plot :math:`u` and :math:`f` for :math:`N=16`. (It suffices to make plots for even :math:`i`.)

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

**b)**
Use the basis from :ref:`fem:approx:exer:tanh:sine1`,
:math:`{\psi}_i(x) = \sin((2i+1)x)`, :math:`i\in{\mathcal{I}_s} = \{0,\ldots,N\}`.
(It suffices to make plots for even :math:`i`.)
Observe that the approximation converges to a piecewise
linear function!

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

**c)**
Use the basis
:math:`{\psi}_i(x) = \sin(2(i+1)x)`, :math:`i\in{\mathcal{I}_s} = \{0,\ldots,N\}`,
and observe that the approximation converges to a piecewise
constant function.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``tanh_sines_boundary_term``.

.. Closing remarks for this Problem

Remarks          (2)
~~~~~~~~~~~~~~~~~~~~

The strange results in b) and c) are due to the choice of
basis. In b), :math:`{\varphi}_i(x)` is an odd function around
:math:`x=\pi/2` and :math:`x=3\pi/2`. No combination of basis functions
is able to approximate the flat regions of :math:`f`.
All basis functions in c) are even around :math:`x=\pi/2` and :math:`x=3\pi/2`,
but odd at :math:`x=0,\pi,2\pi`. With all the sines represented, as in a),
the approximation is not constrained with a particular symmetry
behavior.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:exer:Fourier:

Exercise 8: Fourier series as a least squares approximation
-----------------------------------------------------------

**a)**
Given a function :math:`f(x)` on an interval :math:`[0,L]`, look up the formula
for the coefficients :math:`a_j` and :math:`b_j` in the Fourier series of :math:`f`:

.. math::
        
        f(x) = \frac{1}{2}a_0 +
        \sum_{j=1}^\infty a_j\cos \left(j\frac{2\pi x}{L}\right)
        + \sum_{j=1}^\infty b_j\sin \left(j\frac{2\pi x}{L}\right){\thinspace .}
        

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

**b)**
Let an infinite-dimensional vector space :math:`V` have the basis functions
:math:`\cos j\frac{2\pi x}{L}` for :math:`j=0,1,\dots,\infty` and
:math:`\sin j\frac{2\pi x}{L}` for :math:`j=1,\dots,\infty`.  Show that the least squares
approximation method from the section :ref:`fem:approx:global` leads to a
linear system whose solution coincides with the standard formulas for
the coefficients in a Fourier series of :math:`f(x)` (see also
the section :ref:`fem:approx:global:Fourier`).

.. --- begin hint in exercise ---

**Hint.**
You may choose

.. _Eq:fem:approx:exer:Fourier:basis:

.. math::

    \tag{125}
    {\psi}_{2i} = \cos\left( i\frac{2\pi}{L}x\right),\quad
        {\psi}_{2i+1} = \sin\left( i\frac{2\pi}{L}x\right),
        
        

for :math:`i=0,1,\ldots,N\rightarrow\infty`.

.. --- end hint in exercise ---

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

**c)**
Choose :math:`f(x) = H(x-\frac{1}{2})` on :math:`\Omega=[0,1]`, where :math:`H` is the
Heaviside function: :math:`H(x)=0` for :math:`x<0`, :math:`H(x)=1` for :math:`x>0` and
:math:`H(0)=\frac{1}{2}`. Find the coefficients :math:`a_j` and :math:`b_j` in the
Fourier series for :math:`f(x)`. Plot the sum for :math:`j=2N+1`, where :math:`N=5` and
:math:`N=100`.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``Fourier_ls``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:exer:tanh:Lagrange:

Problem 9: Approximate a steep function by Lagrange polynomials
---------------------------------------------------------------

Use interpolation with uniformly distributed
points and Chebychev nodes to approximate

.. math::
         f(x) = -\tanh(s(x-\frac{1}{2})),\quad x\in [0,1],

by Lagrange polynomials for :math:`s=5` and :math:`s=20`, and :math:`N=3,7,11,15`.
Combine :math:`2\times 2` plots of the approximation for the four
:math:`N` values, and create such figures for the four combinations of
:math:`s` values and point types.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``tanh_Lagrange``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:exer:tanh:Lagrange:regression:

Problem 10: Approximate a steep function by Lagrange polynomials and regression
-------------------------------------------------------------------------------

Redo :ref:`fem:approx:exer:tanh:Lagrange`, but apply a regression
method with :math:`N`-degree Lagrange polynomials and :math:`2N+1`
data points. Recall that
:ref:`fem:approx:exer:tanh:Lagrange` applies :math:`N+1` points and
the resulting approximation interpolates :math:`f` at these points, while
a regression method with more points does not interpolate :math:`f` at the
data points.
Do more points and a regression method help reduce
the oscillatory behavior of Lagrange polynomial approximations?

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``tanh_Lagrange_regression``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:fe:exer:mesh1:

Problem 11: Define nodes and elements
-------------------------------------

Consider a domain :math:`\Omega =[0,2]` divided into the three elements
:math:`[0,1]`, :math:`[1,1.2]`, and :math:`[1.2,2]`.

For P1 and P2 elements, set up the list of coordinates and nodes
(``nodes``) and the numbers of the nodes that belong to each element
(``elements``) in two cases: 1) nodes and elements numbered from left to
right, and 2) nodes and elements numbered from right to left.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``fe_numberings1``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:fe:exer:mesh2:

Problem 12: Define vertices, cells, and dof maps
------------------------------------------------

Repeat :ref:`fem:approx:fe:exer:mesh1`, but define the
data structures ``vertices``, ``cells``, and ``dof_map`` instead of
``nodes`` and ``elements``.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``fe_numberings2``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:fe:exer:defmesh:sparsity:

Problem 13: Construct matrix sparsity patterns
----------------------------------------------

:ref:`fem:approx:fe:exer:mesh1` describes a element mesh
with a total of five elements, but with two different element and
node orderings. For each of the two orderings,
make a :math:`5\times 5` matrix and fill in the entries that will be nonzero.

.. --- begin hint in exercise ---

**Hint.**
A matrix entry :math:`(i,j)` is nonzero if :math:`i` and :math:`j` are nodes in the
same element.

.. --- end hint in exercise ---

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``fe_sparsity_pattern``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:fe:exer:Asinwt:symbolic:

Problem 14: Perform symbolic finite element computations
--------------------------------------------------------

Perform symbolic calculations to find formulas for the coefficient
matrix and right-hand side when approximating :math:`f(x) = \sin (x)` on
:math:`\Omega=[0, \pi]` by two P1 elements of size :math:`\pi/2`.  Solve the
system and compare :math:`u(\pi/2)` with the exact value 1.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``fe_sin_P1``.

.. Hint: wolframalpha or sympy can help with (1-x)*sin(a*x+b),

.. which is the integral

.. that arises on the right-hand side.

.. solution:

.. from fe_approx1D_numint import *

.. c = approximate(sym.sin(x), symbolic=True, d=1, N_e=2, numint=None,

.. Omega=[0,sym.pi])

.. print sym.simplify(c[1,0].subs('h', sym.pi/2))

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:exer:tanh:P1P2:

Problem 15: Approximate a steep function by P1 and P2 elements
--------------------------------------------------------------

Given

.. math::
         f(x) = \tanh(s(x-\frac{1}{2}))

use the Galerkin or least squares method with finite elements to find
an approximate function :math:`u(x)`. Choose :math:`s=20` and try
:math:`N_e=4,8,16` P1 elements and
:math:`N_e=2,4,8` P2 elements.
Integrate :math:`f{\varphi}_i` numerically.

.. --- begin hint in exercise ---

**Hint.**
You can automate the computations by calling the ``approximate`` method
in the ``fe_approx1D_numint`` module.

.. --- end hint in exercise ---

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``fe_tanh_P1P2``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:exer:tanh:P3P4:

Problem 16: Approximate a steep function by P3 and P4 elements
--------------------------------------------------------------

**a)**
Solve :ref:`fem:approx:exer:tanh:P1P2` using :math:`N_e=1,2,4` P3 and P4
elements.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

**b)**
How will an interpolation method work in
this case with the same number of nodes?

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``fe_tanh_P3P4``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:fe:exer:Asinwt:interpol:error:

Exercise 17: Investigate the approximation error in finite elements
-------------------------------------------------------------------

The theory :ref:`(101) <Eq:fem:approx:fe:error:theorem>` from the section :ref:`fem:approx:fe:error` predicts that the error in the Pd
approximation of a function should behave as :math:`h^{d+1}`, where :math:`h` is
the length of the element. Use experiments to verify this asymptotic
behavior (i.e., for small enough :math:`h`).  Choose three examples:
:math:`f(x)=Ae^{-\omega x}` on :math:`[0,3/\omega]`, :math:`f(x) = A\sin (\omega x)` on
:math:`\Omega=[0, 2\pi/\omega]` for constant :math:`A` and :math:`\omega`, and
:math:`f(x)=\sqrt{x}` on :math:`[0,1]`.

.. --- begin hint in exercise ---

**Hint 1.**
Run a series of experiments: :math:`(h_i,E_i)`, :math:`i=0,\ldots,m`, where :math:`E_i`
is the :math:`L^2` norm of the error corresponding to element length :math:`h_i`.
Assume an error model :math:`E=Ch^r` and compute :math:`r` from two successive
experiments:

.. math::
         r_i = \ln (E_{i+1}/E_i)/\ln (h_{i+1}/h_i),\quad i=0,\ldots,m-1{\thinspace .}

Hopefully, the sequence :math:`r_0,\ldots,r_{m-1}` converges to the true
:math:`r`, and :math:`r_{m-1}` can be taken as an approximation to :math:`r`.
Run such experiments for different :math:`d` for the different :math:`f(x)` functions.

.. --- end hint in exercise ---

.. --- begin hint in exercise ---

**Hint 2.**
The ``approximate`` function in ``fe_approx1D_numint.py`` is handy for
calculating the numerical solution. This function returns the
finite element solution as the coefficients :math:`\left\{ {c}_i \right\}_{i\in{\mathcal{I}_s}}`.
To compute :math:`u`, use ``u_glob`` from the same module.
Use the Trapezoidal rule to integrate the :math:`L^2` error:

.. code-block:: python

        xc, u = u_glob(c, vertices, cells, dof_map)
        e = f_func(xc) - u
        L2_error = 0
        e2 = e**2
        for i in range(len(xc)-1):
            L2_error += 0.5*(e2[i+1] + e2[i])*(xc[i+1] - xc[i])
        L2_error = np.sqrt(L2_error)

The reason for this Trapezoidal integration is
that ``u_glob`` returns coordinates ``xc`` and corresponding ``u`` values
where some of the coordinates (the cell vertices) coincides, because
the solution is computed in one element at a time, using all local
nodes. Also note that there are many coordinates in :math:`xc` per cell
such that we can accurately compute the error inside each cell.

.. --- end hint in exercise ---

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``Pd_approx_error``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:fe:exer:Heaviside:

Problem 18: Approximate a step function by finite elements
----------------------------------------------------------

Approximate the step function

.. math::
         f(x) = \left\lbrace\begin{array}{ll}
        0 & 0\leq x < {1/2},\\ 
        1 & {1/2} \leq x \geq {1/2}
        \end{array}\right.
        

by 2, 4, 8, and 16 P1, P2, P3, and P4. Compare approximations visually.

.. --- begin hint in exercise ---

**Hint.**
This :math:`f` can also be expressed in terms of the Heaviside function :math:`H(x)`:
:math:`f(x) = H(x-{1/2})`.
Therefore, :math:`f` can be defined by

.. code-block:: python

        f = sym.Heaviside(x - sym.Rational(1,2))

making the ``approximate`` function in the
``fe_approx1D.py`` module an obvious candidate to solve the
problem. However, ``sympy`` does not handle symbolic integration
with this particular integrand, and the ``approximate`` function faces a problem
when converting ``f`` to a Python function (for plotting) since
``Heaviside`` is not an available function in ``numpy``.

An alternative is to perform hand calculations. This is an instructive
task, but in practice only feasible for few elements and P1 and P2 elements.
It is better to copy the functions ``element_matrix``, ``element_vector``,
``assemble``, and ``approximate`` from the ``fe_approx1D_numint.py`` file
and edit these functions such that they can compute approximations
with ``f`` given as a Python function and not a symbolic expression.
Also assume that ``phi`` computed by the ``basis`` function is a Python
callable function. Remove all instances of the ``symbolic`` variable
and associated code.

.. --- end hint in exercise ---

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``fe_Heaviside_P1P2``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:fe:exer:2Dsines:symbolic:

Exercise 19: 2D approximation with orthogonal functions
-------------------------------------------------------

**a)**
Assume we have basis functions :math:`{\varphi}_i(x,y)` in 2D that are
orthogonal such that :math:`({\varphi}_i,{\varphi}_j)=0` when :math:`i\neq j`.  The
function ``least_squares`` in the file `approx2D.py <http://tinyurl.com/nm5587k/approx/fe_approx2D.py>`__ will then spend much time on computing
off-diagonal terms in the coefficient matrix that we know are zero.
To speed up the computations, make a version ``least_squares_orth`` that
utilizes the orthogonality among the basis functions.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

**b)**
Apply the function to approximate

.. math::
         f(x,y) = x(1-x)y(1-y)e^{-x-y}

on :math:`\Omega = [0,1]\times [0,1]` via basis functions

.. math::
         {\varphi}_i(x,y) = \sin ((p+1)\pi x)\sin((q+1)\pi y),\quad i=q(N_x+1) + p,
        

where :math:`p=0,\ldots,N_x` and :math:`q=0,\ldots,N_y`.

.. --- begin hint in exercise ---

**Hint.**
Get ideas from the function ``least_squares_orth`` in
the section :ref:`fem:approx:global:orth` and
file `approx1D.py <http://tinyurl.com/nm5587k/approx/fe_approx1D.py>`__.

.. --- end hint in exercise ---

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

**c)**
Make a unit test for the ``least_squares_orth`` function.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``approx2D_ls_orth``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:fe:exer:1D:trapez:

Exercise 20: Use the Trapezoidal rule and P1 elements
-----------------------------------------------------

Consider approximation of some :math:`f(x)` on an interval :math:`\Omega` using
the least squares or Galerkin methods with P1 elements. Derive
the element matrix and vector using the
Trapezoidal rule :ref:`(109) <Eq:fem:approx:fe:numint1:trapez>` for calculating
integrals on the reference element. Assemble the contributions, assuming
a uniform cell partitioning, and show that the resulting linear system
has the form :math:`c_i=f(x_{i})` for :math:`i\in{\mathcal{I}_s}`.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``fe_P1_trapez``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:fe:exer:1D:P1:vs:interp:

Exercise 21: Compare P1 elements and interpolation
--------------------------------------------------

We shall approximate the function

.. math::
         f(x) = 1 + \epsilon\sin (2\pi nx),\quad x\in \Omega = [0,1],

where :math:`n\in\mathbb{Z}` and :math:`\epsilon \geq 0`.

**a)**
Plot :math:`f(x)` for :math:`n=1,2,3` and find the wave length of the function.

**b)**
We want to use :math:`N_P` elements per wave length. Show that the number
of elements is then :math:`nN_P`.

**c)**
The critical quantity for accuracy is the number of elements per
wave length, not the element size in itself. It therefore suffices
to study an :math:`f` with just one wave length in :math:`\Omega = [0,1]`.
Set :math:`\epsilon = 0.5`.

Run the least squares or projection/Galerkin method for
:math:`N_P=2,4,8,16,32`. Compute the error :math:`E=||u-f||_{L^2}`.

.. --- begin hint in exercise ---

**Hint.**
Use the ``fe_approx1D_numint`` module to compute :math:`u` and use
the technique from the section :ref:`fem:approx:fe:error` to
compute the norm of the error.

.. --- end hint in exercise ---

**d)**
Repeat the set of experiments in the above point, but
use interpolation/collocation based on the node points to
compute :math:`u(x)` (recall that :math:`c_i` is now simply :math:`f(x_{i})`).
Compute the error :math:`E=||u-f||_{L^2}`.
Which method seems to be most accurate?

Filename: ``fe_P1_vs_interp``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:fe:exer:3D:approx3D:

Exercise 22: Implement 3D computations with global basis functions
------------------------------------------------------------------

Extend the `approx2D.py <http://tinyurl.com/nm5587k/approx/approx2D.py>`__ code to 3D
applying ideas from the section :ref:`fem:approx:3D:global`.
Construct some 3D problem to make a test function for the
implementation.

.. --- begin hint in exercise ---

**Hint.**
Drop symbolic integration since it is in general too slow for 3D problems.
Also use ``scipy.integrate.nquad`` instead of ``sympy.mpmath.quad``
for numerical integration, since it is much faster.

.. --- end hint in exercise ---

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``approx3D``.

.. --- end exercise ---

.. --- begin exercise ---

.. _fem:approx:fe:exer:1D:simpson:

Exercise 23: Use Simpson's rule and P2 elements
-----------------------------------------------

Redo :ref:`fem:approx:fe:exer:1D:trapez`, but use P2
elements and Simpson's rule based on sampling the integrands at
the nodes in the reference cell.

.. removed !bsol ... !esol environment (because of the command-line option --without_solutions)

Filename: ``fe_P2_simpson``.

.. --- end exercise ---