36.3: Function Approximation

Definition

Let \(C[a,b]\) be a vector space of all possible continuous functions over the interval \([a,b]\) with inner product: \(\langle f,g \rangle = \int_a^b f(x)g(x) dx.\)

Now let \(f\) be an element of \(C[a,b]\), and \(W\) be a subspace of \(C[a,b]\). The function \(g \in W\) such that \(\int_a^b \left[ f(x) - g(x) \right]^2 dx\) is a minimum is called the least-squares approximation to \(f\).

The least-squares approximation to \(f\) in the subspace \(W\) can be calculated as the projection of \(f\) onto \(W\):

\[g = proj_Wf \nonumber \]

If \(\{g_1, \ldots, g_n\}\) is an orthonormal basis for \(W\), we can replace the dot product of \(R^n\) by an inner product of the function space and get:

\[prog_Wf = \langle f,g_1 \rangle g_1 + \ldots + \langle f,g_n \rangle g_n \nonumber \]

Polynomial Approximations

An orthogonal bases for all polynomials of degree less than or equal to \(n\) can be computed using Gram-schmidt orthogonalization process. First we start with the following standard basis vectors in \(W\)

\[ \{ 1, x, \ldots, x^n \} \nonumber \]

The Gram-Schmidt process can be used to make these vectors orthogonal. The resulting polynomials on \([−1,1]\) are called Legendre polynomials. The first six Legendre polynomial basis are:

\[1 \nonumber \]

\[x \nonumber \]

\[x^2 -\frac{1}{3} \nonumber \]

\[x^3 - \frac{3}{5}x \nonumber \]

\[x^4 - \frac{6}{7}x^2 + \frac{3}{35} \nonumber \]

\[x^5 - \frac{10}{9}x^3 + \frac{5}{12}x \nonumber \]

Question

What is the least-squares linear approximations of \(f(x)=e^x\) over the interval \([−1,1]\). In other words, what is the projection of \(f\) onto \(W\), where \(W\) is a first order polynomal with basis vectors \(\{1, x\} (\textit{i.e. } n=1)\).

(Hint: You can give the answer in integrals without computing the integrals. Note the Legendre polynomials are not normalized.)

Here is a plot of the equation \(f(x)=e^x\):

%matplotlib inline
import matplotlib.pylab as plt
import numpy as np

#px = np.linspace(-1,1,100)
#py = np.exp(px)
#plt.plot(px,py, color='red');
import sympy as sym
from sympy.plotting import plot
x = sym.symbols('x')
f = sym.exp(x)
plot(f,(x,-1,1))

We can use sympy to compute the integral. The following code compute the definite integral of \(\int_{-1}^1 e^x dx\). In fact, sympy can also compute the indefinite integral by removing the interval.

sym.init_printing()
x = sym.symbols('x')
sym.integrate('exp(x)',(x, -1, 1))
#sym.integrate('exp(x)',(x))

Use sympy to compute the first order polynomial that approximates the function \(e^x\). The following calculates the above approximation written in sympy:

g_0 = sym.integrate('exp(x)*1',(x, -1, 1))/sym.integrate('1*1',(x,-1,1))*1
g_1 = g_0 + sym.integrate('exp(x)*x',(x,-1,1))/sym.integrate('x*x',(x,-1,1))*x
g_1

Plot the original function \(f(x)=e^x\) and its approximation.

p2 = plot(f, g_1,(x,-1,1))

#For fun, I turned this into a function:
x = sym.symbols('x')

def lsf_poly(f, gb = [1,  x], a =-1, b=1):
    proj = 0
    for g in gb:
#        print(sym.integrate(g*f,(x,a,b)))
        proj = proj + sym.integrate(g*f,(x,a,b))/sym.integrate(g*g,(x,a,b))*g
    return proj

lsf_poly(sym.exp(x))

Question

What would a second order approximation look like for this function? How about a fifth order approximation?

#####Start your code here #####
x = sym.symbols('x')
g_2 = 
g_2
#####End of your code here#####

p2 = plot(f, g_2,(x,-1,1))