Lecture 8#
Function approximation with global functions#
Consider a generic function in one dimension
We now want to approximate this function
The other functions
A basis is a collection (or set) of basis functions
A functionspace
For example,
If we say that the function
The set
Consider a very simple example. Let
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
x = sp.Symbol('x')
u = 10*(x-1)**2-1
N = 100
xj = np.linspace(1, 2, N)
plt.figure(figsize=(3, 2))
plt.plot(xj, sp.lambdify(x, u)(xj));

Now create a functionspace consisting of all straight lines
and try to find the approximation
How do we find
There are several possible approaches. The three approaches considered in the book are
The least squares method (variational)
The Galerkin method (variational)
The collocation method (interpolation)
where the first two are variational methods and the last is an interpolation method. We will briefly mention the least squares method and focus on the last two. All methods will allow us to find the best possible
For the first two methods we need to define the
for two real functions
Note
Sometimes the inner product is written as
The inner product is used in defining the
then the
For convenience we define a Python function for the inner product of Sympy functions:
def inner(u, v, domain=(-1, 1), x=x):
return sp.integrate(u*v, (x, domain[0], domain[1]))
The least squares method#
The least squares method requires that we solve
This gives us
Find and display the two equations:
from IPython.display import display
u0, u1 = sp.symbols('u0,u1')
err = u-(u0+u1*x)
E = inner(err, err, (1, 2))
eq1 = sp.Eq(sp.diff(E, u0, 1), 0)
eq2 = sp.Eq(sp.diff(E, u1, 1), 0)
display(eq1)
display(eq2)
Solve the two equations for the two unknown u0
and u1
uhat = sp.solve((eq1, eq2), (u0, u1))
uhat
{u0: -38/3, u1: 10}
Finally, plot the function
plt.figure(figsize=(3, 2))
plt.plot(xj, sp.lambdify(x, u)(xj), 'b', xj, uhat[u0]+uhat[u1]*xj, 'r:')
plt.legend(['$u(x)$', '$u_{N}(x)$']);

Note
Writing
means the same as writing
It is just two different ways of writing. Using an index set, like
The Galerkin method#
The Galerkin method requires that the error is orthogonal to the basis functions (for the
This is not as much work as the least squares method, and even easier to implement in Sympy:
eq1 = sp.Eq(inner(err, 1, (1, 2)), 0)
eq2 = sp.Eq(inner(err, x, (1, 2)), 0)
display(eq1)
display(eq2)
sp.solve((eq1, eq2), (u0, u1))
{u0: -38/3, u1: 10}
Note
For function approximations the Galerkin method always gives the same result as the least squares method. However, the same does not apply when solving PDEs.
The Galerkin method can be formulated as:
Find
which means exactly the same as:
Find
Note
In order to satisfy
In order for this to always be satisfied, we require that
Note
The Galerkin method described above is sometimes referred to as a projection method and may be formulated as: find
Inserting for
and thus
This is a linear algebra system,
where the matrix
Note
The matrix
Note
In the Galerkin literature, the basis function containing the row index (i.e.,
If the basis functions are chosen as orthonormal to each other, then
In that case we simply get
It is more common that the basis functions are orthogonal. In this case the mass matrix becomes
where
The monomial basis
is a basis for the space of all polynomials. However, it is not a good basis. This is because the basis functions are not orthogonal and the mass matrix
A much better polynomial basis are the Legendre polynomials that are defined on the domain
The first 5 Legendre polynomials are plotted below
Show code cell source
plt.figure(figsize=(6, 4))
xj = np.linspace(-1, 1, 100)
legend = []
p = np.zeros(100)
for n in range(5):
l = sp.legendre(n, x)
p[:] = sp.lambdify(x, l)(xj)
plt.plot(xj, p)
legend.append(f'$P_{n}(x)$')
plt.legend(legend);

The Legendre polynomials are
A slightly complicating factor, though, is that the computational domain needs to be
Note
Mapping from a physical to a computational (reference) domain is very common and it is best to get used to it sooner rather than later. With the finite element method this mapping is used all the time and it is described as one of the first topics for the method in Sec 3.1.8 of Introduction to Numerical Methods for Variational Problems.
Note
The Legendre polynomials are not the only basis functions that require a mapping. For example, the Bernstein polynomials (see, Sec. 2.4.3 in the book) use a computational domain
Let
The reverse mapping is
These mappings are valid for any basis function.
Note
When using a basis function like
For Legendre basis functions we use the computational domain
We need to be careful with defining the
Insert for
This is just the regular Galerkin method. But now we introduce
Note that
where
We can rewrite the inner
function to work for a generic interval
def inner(u, v, domain, ref_domain=(-1, 1)):
A, B = ref_domain
a, b = domain
us = u.subs(x, a + (b-a)*(x-A)/(B-A))
return sp.integrate(us*v, (x, A, B))
For example, for our
a, b = 1, 2
uhat = lambda u, j: (2*j+1) * inner(u, sp.legendre(j, x), (a, b))/2
We can now compute any Legendre coefficient from the above lambda function. For
uhat(10*(x-1)**2-1, 0), uhat(10*(x-1)**2-1, 1)
(7/3, 5)
We can thus plot the function
and compare with the exact
from numpy.polynomial import Legendre
plt.figure(figsize=(4, 3))
xj = np.linspace(1, 2, 100)
Xj = -1 + 2/(b-a)*(xj-a)
plt.plot(xj, sp.lambdify(x, u)(xj))
plt.plot(xj, uhat(u, 0) + uhat(u, 1)*Xj, 'r:')
plt.legend(['$10(x-1)^2-1$', f'{Legendre((uhat(u, 0), uhat(u, 1)))}']);

We see that the result is exactly as before, which is not so strange since the span of this Legendre basis is still the space of all linear functions. But add one more Legendre coefficient and the approximation becomes exact because
l2 = sp.lambdify(x, sp.legendre(2, x))
plt.figure(figsize=(4, 3))
plt.plot(xj, sp.lambdify(x, u)(xj))
plt.plot(xj, uhat(u, 0) + uhat(u, 1)*Xj + uhat(u, 2)*l2(Xj), 'r:')
plt.legend(['$10(x-1)^2-1$', f'{Legendre((uhat(u, 0), uhat(u, 1), uhat(u, 2)))}']);

Note
Adding one more Legendre coefficient does not change the first two.
Note
Note that we have not used any mesh points when computing the Legendre coefficients
In the computations above we have used Sympy and computed the inner product integral exactly. For some functions this may be impossible, or take a very long time. Hence, the inner products are often computed with numerical integration. We can rewrite the inner product above to use numerical integration (from scipy). The results are now almost identical, but no longer exact. See the difference below:
def ninner(u, v, domain, ref_domain=(-1, 1)):
from scipy.integrate import quad
A, B = ref_domain
a, b = domain
us = u.subs(x, a + (b-a)*(x-A)/(B-A))
uv = sp.lambdify(x, us*v)
return quad(uv, A, B)[0]
uhatn = lambda u, j: (2*j+1) * ninner(u, sp.legendre(j, x), (1, 2))/2
display((uhat(10*(x-1)**2-1, 0), uhatn(10*(x-1)**2-1, 0)))
display((uhat(10*(x-1)**2-1, 1), uhatn(10*(x-1)**2-1, 1)))
(7/3, 2.3333333333333335)
(5, 5.0)
The numerical integration routine quad is using adaptive quadrature in order to solve the integral. This adaptive integration routine is using a gradually finer mesh until the results of the integration are no longer improving more than a certain tolerance. The default tolerances (both relative and absolute) are
The collocation method#
The collocation method takes a very different approach, but the objective is still to find
The collocation method requires that for some
Inserting for
which can be written as the linear algebra system
where the matrix components
The Lagrange interpolation method described in lecture 7 is a collocation method. The Lagrange basis functions
are defined such that
Hence for the Lagrange interpolation method the matrix
There is no integration and the method is often favored for its simplicity. There is a problem though. How do you choose the collocation points?!
Lets return to the example with
or more simply
We can choose the end points Lagrangebasis
and Lagrangefunction
from lecture 7. The result is then as shown below.
from lagrange import Lagrangebasis, Lagrangefunction
xj = np.linspace(1, 2, 100)
u = 10*(x-1)**2-1
ell = Lagrangebasis([1, 2])
L = Lagrangefunction([u.subs(x, 1), u.subs(x, 2)], ell)
plt.figure(figsize=(3, 2))
plt.plot(xj, sp.lambdify(x, u)(xj), 'b')
plt.plot(xj, sp.lambdify(x, L)(xj), 'r:')
plt.plot([1, 2], [u.subs(x, 1), u.subs(x, 2)], 'bo')
plt.legend(['$u(x)$', '$(x_0, x_1) = (1, 2)$']);

So this seems to be less accurate than the Galerkin method or the least squares method. However, we have not done any integration, and the lower accuracy has come very easily. We can also very easily use one more point and end up with a perfect approximation:
ell = Lagrangebasis([1, 1.5, 2])
L = Lagrangefunction([u.subs(x, 1), u.subs(x, 1.5), u.subs(x, 2)], ell)
plt.figure(figsize=(3, 2))
plt.plot(xj, sp.lambdify(x, u)(xj), 'b')
plt.plot(xj, sp.lambdify(x, L)(xj), 'r:')
plt.plot([1, 1.5, 2], [u.subs(x, 1), u.subs(x, 1.5), u.subs(x, 2)], 'bo')
plt.legend(['$u(x)$', '$(x_0, x_1, x_2) = (1, 1.5, 2)$']);

Lets consider a more difficult function
and attempt to approximate it with Lagrange polynomials on a uniform grid.
u = 1/(1+16*x**2)
N = 16
xj = np.linspace(-1, 1, N+1)
uj = sp.lambdify(x, u)(xj)
ell = Lagrangebasis(xj)
L = Lagrangefunction(uj, ell)
Plot on a fine grid to plot in between mesh points. We know the Lagrange polynomial is exact on the mesh points.
yj = np.linspace(-1, 1, 1000)
plt.figure(figsize=(6, 4))
plt.plot(xj, uj, 'bo')
plt.plot(yj, sp.lambdify(x, u)(yj))
plt.plot(yj, sp.lambdify(x, L)(yj), 'r:')
ax = plt.gca()
ax.set_ylim(-2, 1.1);

Ok, so that did not work out well! Large over and under-shoots in between the mesh points. So Lagrange polynomials are bad, right?
It turns out that due to something called Runge’s phenomenon it is a very bad idea to interpolate on a uniform mesh.
Lets try something else and use a clustering of points near the two edges. Chebyshev points are defined as
Use these mesh points for the Lagrange polynomials and plot the results:
xj = np.cos(np.arange(N+1)*np.pi/N)
uj = sp.lambdify(x, u)(xj)
ell = Lagrangebasis(xj)
L = Lagrangefunction(uj, ell)
plt.figure(figsize=(6, 4))
plt.plot(xj, uj, 'bo')
plt.plot(yj, sp.lambdify(x, u)(yj))
plt.plot(yj, sp.lambdify(x, L)(yj), 'r:');

The agreement is now much better. And the agreement will become better and better the more points you use. Whereas for a uniform grid the agreement will become worse and worse the more points you use.
How about Legendre polynomials? The domain is
Compute the first 40 Legendre coefficients and plot them
uj = lambda u, j: (2*j+1) * inner(sp.legendre(j, x), u, (-1, 1))/2
ul = []
for n in range(40):
ul.append(uj(u, n).n())
plt.figure(figsize=(3, 2))
plt.semilogy(ul, '+');

Notice that the Legendre coefficients are decreasing. This is because the series is converging.
Runge’s phenomenon#
Where does the large ove- and undershoots come from on the uniform mesh? What is the cause of Runge’s phenomenon?
It can be shown that the error in the Lagrange polynomial
for some
From Eq. eqref
{eq-lagrangeerror}it is evident that large errors occur when
Lets use both a uniform mesh and Chebyshev points:
N = 16
xp = np.linspace(-1, 1, N+1)
xc = np.cos(np.arange(N+1)*np.pi/N)
fig = plt.figure(figsize=(10, 3))
plt.plot(xp, np.ones(N+1), 'bo', xc, 2*np.ones(N+1), 'ro')
plt.gca().set_ylim(0, 2.2)
plt.gca().set_yticks([])
plt.legend(['Uniform grid', 'Chebyshev points']);

We see that Chebyshev points are clustered near the edges. Lets construct
xp = np.linspace(-1, 1, N+1); xc = np.cos(np.arange(N+1)*np.pi/N)
pp = np.poly(xp); pc = np.poly(xc)
xj = np.linspace(-1, 1, 1000)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 3))
ax1.plot(xj, np.polyval(pp, xj)); ax2.plot(xj, np.polyval(pc, xj))
ax1.plot(xp, np.zeros(N+1), 'bo'); ax2.plot(xc, np.zeros(N+1), 'bo')
ax1.set_title('Uniform grid'); ax2.set_title('Chebyshev grid');

We note that a uniform mesh leads to large oscillations near the edges, whereas Chebyshev points leads to a polynomial with near uniform oscillations.
The oscillations on the uniform mesh grow even larger when increasing
Boundary issues#
Lets consider a slightly different problem
and attempt to find
The sines are orthogonal such that the mass matrix becomes diagonal
Hence we can easily get the coefficients with the Galerkin method
Lets implement and check the results:
u = 10*(x-1)**2-1
uhat = lambda u, i: 2*inner(u, sp.sin((i+1)*sp.pi*x), (0, 1), (0, 1))
ul = []
for i in range(15):
ul.append(uhat(u, i).n())
ul = np.array(ul, dtype=float)
def uN(uh, xj):
N = len(xj)
uj = np.zeros(N)
for i, ui in enumerate(uh):
uj[:] += ui*np.sin((i+1)*np.pi*xj)
return uj
xj = np.linspace(0, 1, 100)
plt.figure(figsize=(6,4))
plt.plot(xj, 10*(xj-1)**2-1, 'b')
plt.plot(xj, uN(ul[:(3+1)], xj), 'g:')
plt.plot(xj, uN(ul[:(12+1)], xj), 'r-')
plt.legend(['$u(x)$', 'N=3', 'N=12']);

Evidently, we get large oscillations and the solution does not seem to converge. The problem is that the basis functions are all such that
since
will always lead to
A possible solution to this problem is to add a function to the right hand side above that satisfies the boundary conditions:
where
Note
The function
We get the slightly untraditional variational problem: find
We get
and using the orthogonality of the sines we get
Compared with Eq. (21) this is only a minor modification, but it makes a world of difference. Lets implement it:
uc = []
b = u.subs(x, 0)*(1-x) + u.subs(x, 1)*x
for i in range(15):
uc.append(uhat(u-b, i).n())
uc = np.array(uc, dtype=float)
def uNc(uh, xj):
N = len(xj)
uj = u.subs(x, 0)*(1-xj) + u.subs(x, 1)*xj
for i, ui in enumerate(uh):
uj[:] += ui*np.sin((i+1)*np.pi*xj)
return uj
plt.figure(figsize=(6,4))
plt.plot(xj, 10*(xj-1)**2-1, 'b')
plt.plot(xj, uNc(uc[:(3+1)], xj), 'g:')
plt.plot(xj, uNc(uc[:(12+1)], xj), 'r-')
plt.legend(['$u(x)$', 'N=3', 'N=12']);

Notice the very good agreement. The expansion coefficients are now converging, which we can see by printing out the 15 computed values
print(uc)
[-2.58012275e+00 0.00000000e+00 -9.55601020e-02 0.00000000e+00
-2.06409820e-02 0.00000000e+00 -7.52222377e-03 0.00000000e+00
-3.53926304e-03 0.00000000e+00 -1.93848441e-03 0.00000000e+00
-1.17438450e-03 0.00000000e+00 -7.64480816e-04]
These should be compared with the expansion coefficients without the boundary adjustment:
print(ul)
[2.51283542 3.18309886 1.60209262 1.59154943 0.99795065 1.06103295
0.72004323 0.79577472 0.56234498 0.63661977 0.46105771 0.53051648
0.39059163 0.45472841 0.33876606]
The sine basis functions need modification because they are are all zero on the boundary. Legendre polynomials, on the other hand, can be used without any boundary adjustments since
We get a minor modification to the Legendre coefficients because of the new physical domain, but otherwise there is no issue and again we get a perfect approximation using three Legendre polynomials:
a, b = 0, 1
uhat = lambda u, j: (2*j+1) * inner(u, sp.legendre(j, x), (0, 1))/2
uL = []
for i in range(6):
uL.append(uhat(u, i))
print(uL)
[7/3, -5, 5/3, 0, 0, 0]
Weekly assignments#
Experiment with the Galerkin and collocation methods and approximate the global functions
where
Experiment using either Legendre polynomials or sines as basis functions for the Galerkin method and Lagrange polynomials for the collocation method. Use both exact and numerical integration for the Galerkin method. Measure the error as a function of
Hint
Number 5 is difficult because of several factors. i) The Bessel function is expensive to integrate exactly and you need to use numerical integration. ii) It requires approximately 80 Legendre coefficients to completely converge. The Legendre polynomials returned by Sympy are such that when simply lambdified and evaluated, they will lead to severe roundoff errors. As such, it is better to use Numpy’s Legendre class for evaluation in the numerical integral.
The approximation of the Bessel function with Legendre polynomials is shown below for