Function basis methods
======================


```
%matplotlib inline
```


```
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

matplotlib.rcParams.update({'font.size': 14})
```

Introduction to function basis methods
--------------------------------------

### Boundary Value Problems

Considering simple boundary value problem

$$y'' = f(x, y, y'), \quad y(a) = A, \,\, y(b) = B, \quad x \in [a,b].$$

Boundary conditions are only examples here.

Have considered the shooting method (accurate, efficient, may not work)
and the finite difference (relaxation) method (not particularly accurate
or efficient, nearly always works). Both methods require a *grid* with
$n$ points. Accuracy of method depends on $h \propto n^{-1}$.

Can instead use a function basis, independently of a grid.

### Function basis expansion

Aim to solve problem

$${\cal L} y = f,$$

where ${\cal L}$ is a differential operator.

Function basis methods *assume* the solution $y(x)$ can be written

$$y(x) = \sum_{j} c_j u_j (x);$$

constants $c_j$ are *basis coefficients*, $u_j(x)$ are *basis
functions*.

*Choose* “simple” basis functions $u_j$. Get *approximate* solution by
truncating series:

$$y(x) = \sum_{j}^n c_j u_j (x).$$

Additional conditions give linear system for unknowns ($c_j$), defining
solution everywhere.

Collocation methods
-------------------

Simplest function basis method uses a grid. *Collocation* method: fix
$c_j$ by insisting that BVP satisfied exactly at fixed set of points.

1.   Assume approximate solution has form ($n$ fixed)

    $$\sum_{j}^n c_j u_j (x).$$
2.   Boundary conditions $\implies$ two coefficients (e.g. $c_{0, 1}$).
3.  Evaluate ${\cal L} y$ at collocation points $\{x_j\}$ $\implies$ linear system.
4.  Solve the linear system to give the basis coefficients $c_j$.

### Example

We consider the problem (with exact solution $\exp(x)$)

$$y'' - y = 0, \quad y(0) = 1, \,\, y(1) = e.$$

Choose basis $u_j = x^j \implies y = \sum c_j x^j$. Boundary conditions:

$$c_0 = 1, \quad c_1 = e - \sum_{j \ne 1}^n c_j.$$

Differential operator ${\cal L} y = y'' - y$:

$$\sum_{j=2}^n c_j \left( j (j - 1) x^{j-2} \right) - \sum_{k=0}^n
c_k x^k = 0.$$

Linear system $A {\boldsymbol{c}} = {\boldsymbol{b}}$ ($\{ x_l\}$
collocation points)

$$A_{j l} = \left( j (j - 1) x_l^{j-2} \right) - x_l^j , \quad
{\boldsymbol{b}} \text{ from boundary conditions.}$$

When using only three polynomial basis functions we have

$$c_0 = 1, \quad c_1 = e - 1 - c_2.$$

We then evaluate the system at the collocation point $x = 1/2$ to get

$$\begin{aligned}
   && 0 & = 2 \cdot 1 \cdot \left(\tfrac{1}{2}\right)^0 c_2 - 1 -
\tfrac{1}{2} \left( e - 1 - c_2 \right) -
\left(\tfrac{1}{2}\right)^2 c_2 \\ 
   \Rightarrow && c_2 & = \tfrac{2}{9} (1 + e).
  \end{aligned}$$

This, combined with boundary conditions, gives approximate solution

$$y = 1 + \tfrac{1}{9} \left( ( 7e - 11) x + 2(1 + e) x^2 \right).$$

### Example: 3

With only a few points the method is very accurate.

Use Chebyshev collocation points

$$x_k = \tfrac{1}{2} \left( 1 + \cos \left( \frac{(k-1) \pi}{n-1}
\right) \right)$$
to check convergence with $n$.

Convergence is very fast.


```
x_exact = np.linspace(0.0, 1.0, 5000)
y_exact = np.exp(x_exact)
```


```
def CollocationPoly(nbasis):
    """Find the coefficients in the expansion for the above problem."""
    
    # Collocation points are the Chebyshev points
    x = 0.5 * (1.0 + np.cos(np.pi * np.array(range(nbasis-1,-1,-1)) / (nbasis-1)))
        
    b = np.zeros_like(x) # The RHS vector
    A = np.zeros((nbasis, nbasis)) # The matrix
    
    # Use the left boundary to fix the first coefficient
    b[0] = 1.0
    A[0, 0] = 1.0
    # Use the right boundary to fix the last coefficient
    b[-1] = np.exp(1.0)
    for j in range(nbasis):
        A[-1, j] = 1.0
    # Fill the rest of the matrix
    for i in range(1, nbasis-1):
        for j in range(nbasis):
            # The -y' term
            A[i, j] -= x[i]**(j)
        for j in range(2, nbasis):
            # The y'' term
            A[i, j] += j*(j-1)*x[i]**(j-2)
    
    # Solve for coefficients
    c = np.linalg.solve(A, b)
        
    return c
```


```
# Do the simple case with 3 basis functions
c3 = CollocationPoly(3)
y3 = np.zeros_like(x_exact)
for j in range(3):
    y3 += c3[j] * x_exact**j

# Plot result
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111)
ax.plot(x_exact, y_exact, 'k-', label = "Exact solution")
ax.plot(x_exact, y3, 'b--', label = "Using 3 collocation points")
ax.legend(loc=2)
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$");
```


```
# Check the convergence with resolution.
nbasis = range(3, 16)
collocation_error = np.zeros((len(nbasis),))
for i in range(len(nbasis)):
    c = CollocationPoly(nbasis[i])
    yn = np.zeros_like(x_exact)
    for j in range(nbasis[i]):
        yn += c[j] * x_exact**j
    collocation_error[i] = np.linalg.norm(yn - y_exact, 2)
# Plot errors
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111)
ax.semilogy(nbasis, collocation_error, 'kx')
ax.set_xlabel("Number of basis functions")
ax.set_ylabel(r"$\|$"+"Error"+r"$\|_2$");
```

Norm method
-----------

Collocation

1.  not independent of grid;

2.  convergence depends on location of points;

3.  BVP solved only at points, otherwise unconstrained.

Instead try to get best solution “on average”. Minimize average error in
some norm,

$$\| {\cal L} (y) (x) - f(x) \|.$$

Assume a function basis expansion

$$y = \sum_{j}^n c_j u_j (x).$$

Gives a minimization problem for the coefficients $c_j$.

### Example

For the example given above we have

$${\cal L} = y'' - y, \quad f = 0,$$

and hence we want to minimize

$$\| y'' - y \|.$$

Norm is over whole interval – so integrate. Typically use 2-norm:

$$\begin{aligned}
F(c_j) & = \| y'' - y \|_2^2 \\
& = \int_0^1 \left[ y'' - y \right]^2 \, \text{d}x.
  \end{aligned}$$

This is function of $c_j$ using function basis assumption.

### Example: 2

As normal fix two coefficients using boundary conditions. Using three
basis functions

$$c_0 = 1, \quad c_1 = e - 1 - c_2.$$

Explicitly computing the quadratic form gives us the “average error”

$$F(c_j) = \tfrac{47}{10} c_2^2 - \tfrac{13}{6} (1 + e) c_2 +
\tfrac{1+e+e^2}{3}.$$
To minimize, differentiate with respect to $c_2$ and set to zero, giving

$$c_2 = \tfrac{65}{282} (1 + e).$$

### Example: 3

Even with a few coefficients the result is very accurate.

Additional analysis required to set up general problem.

Ritz methods
------------

General framework: working in Hilbert space $L_2$ define *inner product*

$$< u, v > = \int_a^b u \cdot v \, \text{d} x.$$

Need ${\cal L}$ to be symmetric and positive definite.

Inner product can *measure distance* on $L_2$: i.e., measure the
distance to the “exact solution” of the ODE.

So define *energy* of the element $u$ as

$$< {\cal L}(u), u > \,\, \ge 0,$$

where inequality follows by integration by parts. Then to get solution
minimize

$$J(u) = < {\cal L}(u), u > - 2 < f, u >.$$

### Ritz method applied to a function basis

Can minimize the functional

$$J(y) = < {\cal L}(y), y > - 2 < f, y >$$

when $y$ is given with respect to a function basis,

$$y = \sum_j^n c_j u_j(x).$$

Linearity means this is equivalent to minimizing quadratic form

$$J(y) = \sum_{m,k}^n c_m c_k < {\cal L}(u_m), u_k > - 2 \sum_m c_m
< f, u_m >.$$

Know $u_j$, so re-express as a condition on $c_j$.

### Ritz method applied to a function basis

$$y = \sum_j^n c_j u_j(x).$$

$$J(y) = \sum_{m,k}^n c_m c_k < {\cal L}(u_m), u_k > - 2 \sum_m c_m
< f, u_m >.$$
Minimizing this functional requires

$$\begin{align}
&& \frac{\partial{}}{\partial{c_m}} J & = 0, & m = 1, \dots, n, \\
\Rightarrow && \sum_m^n c_m < {\cal L}(u_m), u_k > & = < f, u_k >,
& k = 1, \dots, n.
  \end{align}$$

This is just a linear system to solve.

### Example

For the boundary value problem

$$y'' = -\cos(x), \quad y(0) = 0 = y(\pi)$$

we have the exact solution $y(x) = \cos(x) + 2 x /\pi - 1$. Symmetry of
problem suggests the function basis

$$u_j = \sin(2 j x):$$

satisfies boundary conditions and symmetry.

We need to satisfy

$$\sum_m^n c_m < {\cal L}(u_m), u_k >  = < f, u_k >.$$

Vector of unknowns $c_m$; known matrix has elements $A_{m,k} = <
  {\cal L}(u_m), u_k >$; known vector has elements $< f, u_k >$.

### Example: 2


The *matrix* $A_{m,k} = < {\cal L}(u_m), u_k >$ has elements

$$\begin{aligned}
< {\cal L}(u_m), u_k > & = -\int_0^{\pi} 4 m^2 \sin(2 m x)
\sin(2 k
x) \, \text{d} x \\
& =
\begin{cases}
  -2 \pi k^2, & \text{if $m = k$}, \\
  0, & \text{if $m \ne k$}
\end{cases}.
\end{aligned}$$


The *vector* $< f, u_k >$ has elements

$$\begin{aligned}
< f, u_k > & = \int_0^{\pi} -\cos(x) \sin(2 k x) \, \text{d} x \\
& = \frac{4 k}{4 k^2 - 1}.
\end{aligned}$$


$$\begin{aligned}
< {\cal L}(u_m), u_k > & = \begin{cases}
  -2 \pi k^2, & \text{if $m = k$}, \\
  0, & \text{if $m \ne k$}
\end{cases}, \\
< f, u_k > & = \frac{4 k}{4 k^2 - 1}.
\end{aligned}$$

This gives a *diagonal* system which can be solved to give

$$c_m = \frac{2}{\pi m (4 m^2 - 1)}$$

and hence the approximation

$$y = \frac{2}{\pi} \sum_{m=1}^n \frac{\sin(2 m x)}{m (4 m^2 - 1)}$$

which converges to the Fourier series in the limit $n \rightarrow
\infty$.

### Example: 3

In fact even with very few coefficients the Ritz method is accurate.

It also converges quickly.


```
x_exact = np.linspace(0.0, np.pi, 5000)
y_exact = np.cos(x_exact) + 2.0 * x_exact / np.pi - 1.0
```


```
def RitzCoefficients(nbasis):
    """Returns the Ritz coeffcients for the problem above. Note: cheats and use the form computed analytically above."""
    
    A = np.zeros((nbasis,nbasis))
    b = np.zeros((nbasis,))
    for i in range(nbasis):
        b[i] = 4.0 * (i+1) / (4.0 * (i+1)**2 - 1.0)
        A[i, i] = 2.0 * np.pi * (i+1)**2
    c = np.linalg.solve(A, b)
    
    return c
```


```
# Plot the result with 3 coefficients
c3 = RitzCoefficients(3)
y3 = np.zeros_like(x_exact)
for j in range(3):
    y3 += c3[j] * np.sin(2.0*(j+1)*x_exact)

# Plot result
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111)
ax.plot(x_exact, y_exact, 'k-', label = "Exact solution")
ax.plot(x_exact, y3, 'b--', label = "Using 3 collocation points")
ax.legend(loc=1)
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$");
```


```
# Check the convergence with resolution.
nbasis = range(3, 33, 2)
collocation_error = np.zeros((len(nbasis),))
for i in range(len(nbasis)):
    c = RitzCoefficients(nbasis[i])
    yn = np.zeros_like(x_exact)
    for j in range(nbasis[i]):
        yn += c[j] * np.sin(2.0*(j+1)*x_exact)
    collocation_error[i] = np.linalg.norm(yn - y_exact, 2)
# Plot errors
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111)
ax.loglog(nbasis, collocation_error, 'kx')
ax.set_xlabel("Number of basis functions")
ax.set_ylabel(r"$\|$"+"Error"+r"$\|_2$");
```

Summary
=======


-   When extreme accuracy is essential then function basis methods are
often used.
-   Collocation methods are popular; given the right choice of basis and
collocation points the convergence can be faster than any
polynomial, giving floating point accuracy with a few dozen basis
functions at most.
-   Norm and Ritz type methods do not use a grid at all; this can make
them useful in complex domains, particular as a small part of a
larger method (see finite elements).
-   The complexity of these methods, particularly norm and Ritz, makes
the work involved in setting up the problem quite extreme. The
actual coding of the method is, however, very small for the
impressive accuracy.

2f1c866a9bd75e2ba5ef8dfb5ae6df91c3127f09
ipynb
Jupyter Notebook
Lectures/20 - Function Basis Methods.ipynb
josh-gree/NumericalMethods
03cb91114b3f5eb1b56916920ad180d371fe5283
CC-BY-3.0
2015-02-12T19:51:52.000Z
2022-03-26T15:34:11.000Z
Lectures/20 - Function Basis Methods.ipynb
josh-gree/NumericalMethods
03cb91114b3f5eb1b56916920ad180d371fe5283
CC-BY-3.0
2017-05-24T19:49:52.000Z
2018-01-23T21:40:42.000Z
Lectures/20 - Function Basis Methods.ipynb
josh-gree/NumericalMethods
03cb91114b3f5eb1b56916920ad180d371fe5283
CC-BY-3.0
2015-01-05T13:30:47.000Z
2022-02-15T09:59:39.000Z
Qwen/Qwen-72B
1. YES
2. YES