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 \( (a,b) \) form a vector space by showing that all the axioms of a vector space are satisfied.
b) Prove that all linear functions of the form \( ax+b \) constitute a vector space, \( a,b\in\Real \).
c) Show that all quadratic functions of the form \( 1 + ax^2 + bx \) do not constitute a vector space.
d) Check out the topic of inner product spaces. Suggest a possible inner product for the space of all linear functions of the form \( ax+b \), \( a,b\in\Real \), defined on some interval \( \Omega=[A,B] \). Show that this particular inner product satisfies the general requirements of an inner product in a vector space.
Filename: linalg1.
Given \( \f = (1,1,1) \) in \( \Real^3 \), find the best approximation vector \( \u \) in the plane spanned by the unit vectors \( (1,0) \) and \( (0,1) \). Repeat the calculations using the vectors \( (2,1) \) and \( (1,2) \).
Filename: vec111_approx.
Given the function \( f(x)=1 + 2x(1-x) \) on \( \Omega=[0,1] \), we want to find an approximation in the function space $$ V = \hbox{span}\{1, \sin(\pi x)\}\tp$$
a) Sketch or plot \( f(x) \). Think intuitively how an expansion in terms of the basis functions of \( V \), \( \baspsi_0(x)=1 \), \( \baspsi_1(x)=\sin(\pi x) \), can be construction to yield a best approximation to \( f \). Or phrased differently, see if you can guess the coefficients \( c_0 \) and \( c_1 \) in the expansion $$ u(x) = c_0\baspsi_0 + c_1\baspsi_1 = c_0 + c_1\sin (\pi x)\tp$$ Compute the \( L^2 \) error \( ||f-u||_{L^2}=(\int_0^1(f-u)^2\dx)^{1/2} \).
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.
b) Perform the hand calculations for a least squares approximation.
Filename: parabola_sin.
Let \( V \) be a function space with basis functions \( x^i \), \( i=0,1,\ldots,N \). Find the best approximation to \( f(x)=\exp(-x) \) on \( \Omega =[0,8] \) among all functions in \( V \) for \( N=2,4,6 \). Illustrate the three approximations in three separate plots.
Apply the lest_squares and
comparison_plot functions in the approx1D.py module as these
make the exercise easier to solve.
Filename: exp_powers.
In this exercise we want to approximate the sine function by polynomials of order \( N+1 \). Consider two bases: $$ \begin{align*} V_1 &= \{x, x^3, x^5, \ldots, x^{N-2}, x^N \}, \\ V_2 &= \{1,x,x^2,x^3,\ldots, x^N\}\tp \end{align*} $$ The basis \( V_1 \) is motivated by the fact that the Taylor polynomial approximation to the sine function has only odd powers, while \( 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 \( f(x)=\sin(x) \) among all functions in \( V_1 \) and \( V_2 \) on two domains of increasing sizes: \( \Omega_{1,k} = [0, k\pi] \), \( k=2,3\ldots,6 \) and \( \Omega_{2,k} = [-k\pi /2, k\pi/2] \), \( k=2,3,4,5 \). Make plots for all combinations of \( V_1 \), \( V_2 \), \( \Omega_1 \), \( \Omega_2 \), \( k=2,3,\ldots,6 \).
Add a plot of the \( N \)-th degree Taylor polynomial approximation of \( \sin(x) \) around \( x=0 \).
You can make a loop over \( V_1 \) and \( V_2 \), a loop over
\( \Omega_1 \) and \( \Omega_2 \), and a loop over \( k \). Inside the loops,
call the functions least_squares and
comparison_plot from the approx1D module.
\( N=7 \) is a suggested value.
Filename: sin_powers.
Find the best approximation of \( f(x) = \tanh (s(x-\pi)) \) on \( [0, 2\pi] \) in the space \( V \) with basis \( \baspsi_i(x) = \sin((2i+1)x) \), \( i\in\If = \{0,\ldots,N\} \). Make a movie showing how \( u=\sum_{j\in\If}c_j\baspsi_j(x) \) approximates \( f(x) \) as \( N \) grows. Choose \( s \) such that \( f \) is steep (\( s=20 \) is appropriate).
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.
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).
Filename: tanh_sines.
Approximation of a discontinuous (or steep) \( f(x) \) by sines, results in slow convergence and oscillatory behavior of the approximation close to the abrupt changes in \( f \). This is known as the Gibb's phenomenon.
We study the same approximation problem as in Problem 6: Approximate a steep function by sines. Since \( \baspsi_i(0)=\baspsi_i(2\pi)=0 \) for all \( i \), \( u=0 \) at the boundary points \( x=0 \) and \( x=2\pi \), while \( f(0)=-1 \) and \( f(2\pi)=1 \). This discrepancy at the boundary can be removed by adding a boundary function \( B(x) \): $$ u(x) = B(x) + \sum_{j\in\If} c_j\baspsi_j(x), $$ where \( B(x) \) has the right boundary values: \( B(x_L)=f(x_L) \) and \( B(x_R)=f(x_R) \), with \( x_L=0 \) and \( x_R=2\pi \) as the boundary points. A linear choice of \( B(x) \) is $$ B(x) = \frac{(x_R-x)f(x_L) + (x-x_L)f(x_R)}{x_R-x_L}\tp$$
a) Use the basis \( \baspsi_i(x) = \sin((i+1)x) \), \( i\in\If = \{0,\ldots,N\} \) and plot \( u \) and \( f \) for \( N=16 \). (It suffices to make plots for even \( i \).)
b) Use the basis from Problem 6: Approximate a steep function by sines, \( \baspsi_i(x) = \sin((2i+1)x) \), \( i\in\If = \{0,\ldots,N\} \). (It suffices to make plots for even \( i \).) Observe that the approximation converges to a piecewise linear function!
c) Use the basis \( \baspsi_i(x) = \sin(2(i+1)x) \), \( i\in\If = \{0,\ldots,N\} \), and observe that the approximation converges to a piecewise constant function.
Filename: tanh_sines_boundary_term.
The strange results in b) and c) are due to the choice of basis. In b), \( \basphi_i(x) \) is an odd function around \( x=\pi/2 \) and \( x=3\pi/2 \). No combination of basis functions is able to approximate the flat regions of \( f \). All basis functions in c) are even around \( x=\pi/2 \) and \( x=3\pi/2 \), but odd at \( x=0,\pi,2\pi \). With all the sines represented, as in a), the approximation is not constrained with a particular symmetry behavior.
a) Given a function \( f(x) \) on an interval \( [0,L] \), look up the formula for the coefficients \( a_j \) and \( b_j \) in the Fourier series of \( f \): $$ \begin{equation*} 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)\tp \end{equation*} $$
b) Let an infinite-dimensional vector space \( V \) have the basis functions \( \cos j\frac{2\pi x}{L} \) for \( j=0,1,\dots,\infty \) and \( \sin j\frac{2\pi x}{L} \) for \( j=1,\dots,\infty \). Show that the least squares approximation method from the section Approximation of functions leads to a linear system whose solution coincides with the standard formulas for the coefficients in a Fourier series of \( f(x) \) (see also the section Fourier series).
You may choose $$ \begin{equation} \baspsi_{2i} = \cos\left( i\frac{2\pi}{L}x\right),\quad \baspsi_{2i+1} = \sin\left( i\frac{2\pi}{L}x\right), \tag{125} \end{equation} $$ for \( i=0,1,\ldots,N\rightarrow\infty \).
c) Choose \( f(x) = H(x-\half) \) on \( \Omega=[0,1] \), where \( H \) is the Heaviside function: \( H(x)=0 \) for \( x < 0 \), \( H(x)=1 \) for \( x>0 \) and \( H(0)=\half \). Find the coefficients \( a_j \) and \( b_j \) in the Fourier series for \( f(x) \). Plot the sum for \( j=2N+1 \), where \( N=5 \) and \( N=100 \).
Filename: Fourier_ls.
Use interpolation with uniformly distributed points and Chebychev nodes to approximate $$ \begin{equation*} f(x) = -\tanh(s(x-\half)),\quad x\in [0,1],\end{equation*} $$ by Lagrange polynomials for \( s=5 \) and \( s=20 \), and \( N=3,7,11,15 \). Combine \( 2\times 2 \) plots of the approximation for the four \( N \) values, and create such figures for the four combinations of \( s \) values and point types.
Filename: tanh_Lagrange.
Redo Problem 9: Approximate a steep function by Lagrange polynomials, but apply a regression method with \( N \)-degree Lagrange polynomials and \( 2N+1 \) data points. Recall that Problem 9: Approximate a steep function by Lagrange polynomials applies \( N+1 \) points and the resulting approximation interpolates \( f \) at these points, while a regression method with more points does not interpolate \( f \) at the data points. Do more points and a regression method help reduce the oscillatory behavior of Lagrange polynomial approximations?
Filename: tanh_Lagrange_regression.
Consider a domain \( \Omega =[0,2] \) divided into the three elements \( [0,1] \), \( [1,1.2] \), and \( [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.
Filename: fe_numberings1.
Repeat Problem 11: Define nodes and elements, but define the
data structures vertices, cells, and dof_map instead of
nodes and elements.
Filename: fe_numberings2.
Problem 11: Define nodes and elements 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 \( 5\times 5 \) matrix and fill in the entries that will be nonzero.
A matrix entry \( (i,j) \) is nonzero if \( i \) and \( j \) are nodes in the same element.
Filename: fe_sparsity_pattern.
Perform symbolic calculations to find formulas for the coefficient matrix and right-hand side when approximating \( f(x) = \sin (x) \) on \( \Omega=[0, \pi] \) by two P1 elements of size \( \pi/2 \). Solve the system and compare \( u(\pi/2) \) with the exact value 1.
Filename: fe_sin_P1.
Given $$ \begin{equation*} f(x) = \tanh(s(x-\half))\end{equation*} $$ use the Galerkin or least squares method with finite elements to find an approximate function \( u(x) \). Choose \( s=20 \) and try \( N_e=4,8,16 \) P1 elements and \( N_e=2,4,8 \) P2 elements. Integrate \( f\basphi_i \) numerically.
You can automate the computations by calling the approximate method
in the fe_approx1D_numint module.
Filename: fe_tanh_P1P2.
a) Solve Problem 15: Approximate a steep function by P1 and P2 elements using \( N_e=1,2,4 \) P3 and P4 elements.
b) How will an interpolation method work in this case with the same number of nodes?
Filename: fe_tanh_P3P4.
The theory (101) from the section Computing the error of the approximation predicts that the error in the Pd approximation of a function should behave as \( h^{d+1} \), where \( h \) is the length of the element. Use experiments to verify this asymptotic behavior (i.e., for small enough \( h \)). Choose three examples: \( f(x)=Ae^{-\omega x} \) on \( [0,3/\omega] \), \( f(x) = A\sin (\omega x) \) on \( \Omega=[0, 2\pi/\omega] \) for constant \( A \) and \( \omega \), and \( f(x)=\sqrt{x} \) on \( [0,1] \).
Run a series of experiments: \( (h_i,E_i) \), \( i=0,\ldots,m \), where \( E_i \) is the \( L^2 \) norm of the error corresponding to element length \( h_i \). Assume an error model \( E=Ch^r \) and compute \( r \) from two successive experiments: $$ r_i = \ln (E_{i+1}/E_i)/\ln (h_{i+1}/h_i),\quad i=0,\ldots,m-1\tp$$ Hopefully, the sequence \( r_0,\ldots,r_{m-1} \) converges to the true \( r \), and \( r_{m-1} \) can be taken as an approximation to \( r \). Run such experiments for different \( d \) for the different \( f(x) \) functions.
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 \( \sequencei{c} \).
To compute \( u \), use u_glob from the same module.
Use the Trapezoidal rule to integrate the \( L^2 \) error:
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 \( xc \) per cell
such that we can accurately compute the error inside each cell.
Filename: Pd_approx_error.
Approximate the step function $$ \begin{equation*} f(x) = \left\lbrace\begin{array}{ll} 0 & 0\leq x < \halfi,\\ 1 & \halfi \leq x \geq \halfi \end{array}\right. \end{equation*} $$ by 2, 4, 8, and 16 P1, P2, P3, and P4. Compare approximations visually.
This \( f \) can also be expressed in terms of the Heaviside function \( H(x) \): \( f(x) = H(x-\halfi) \). Therefore, \( f \) can be defined by
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.
Filename: fe_Heaviside_P1P2.
a)
Assume we have basis functions \( \basphi_i(x,y) \) in 2D that are
orthogonal such that \( (\basphi_i,\basphi_j)=0 \) when \( i\neq j \). The
function least_squares in the file 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.
b) Apply the function to approximate $$ f(x,y) = x(1-x)y(1-y)e^{-x-y}$$ on \( \Omega = [0,1]\times [0,1] \) via basis functions $$ \basphi_i(x,y) = \sin ((p+1)\pi x)\sin((q+1)\pi y),\quad i=q(N_x+1) + p, $$ where \( p=0,\ldots,N_x \) and \( q=0,\ldots,N_y \).
Get ideas from the function least_squares_orth in
the section Orthogonal basis functions and
file approx1D.py.
c)
Make a unit test for the least_squares_orth function.
Filename: approx2D_ls_orth.
Consider approximation of some \( f(x) \) on an interval \( \Omega \) using the least squares or Galerkin methods with P1 elements. Derive the element matrix and vector using the Trapezoidal rule (109) 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 \( c_i=f(\xno{i}) \) for \( i\in\If \).
Filename: fe_P1_trapez.
We shall approximate the function $$ f(x) = 1 + \epsilon\sin (2\pi nx),\quad x\in \Omega = [0,1],$$ where \( n\in\Integer \) and \( \epsilon \geq 0 \).
a) Plot \( f(x) \) for \( n=1,2,3 \) and find the wave length of the function.
b) We want to use \( N_P \) elements per wave length. Show that the number of elements is then \( 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 \( f \) with just one wave length in \( \Omega = [0,1] \). Set \( \epsilon = 0.5 \).
Run the least squares or projection/Galerkin method for \( N_P=2,4,8,16,32 \). Compute the error \( E=||u-f||_{L^2} \).
Use the fe_approx1D_numint module to compute \( u \) and use
the technique from the section Computing the error of the approximation to
compute the norm of the error.
d) Repeat the set of experiments in the above point, but use interpolation/collocation based on the node points to compute \( u(x) \) (recall that \( c_i \) is now simply \( f(\xno{i}) \)). Compute the error \( E=||u-f||_{L^2} \). Which method seems to be most accurate?
Filename: fe_P1_vs_interp.
Extend the approx2D.py code to 3D applying ideas from the section Extension to 3D. Construct some 3D problem to make a test function for the implementation.
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.
Filename: approx3D.
Redo Exercise 20: Use the Trapezoidal rule and P1 elements, but use P2 elements and Simpson's rule based on sampling the integrands at the nodes in the reference cell.
Filename: fe_P2_simpson.