Approximation¶
Contents¶
To include the figures in this notebook, instead of in an independent window, we run the following command:
%matplotlib inline
With from __future__ import division we force the result of division with / to be always float. If we do not include this, the division of two integers with / would be an integer.
We also import the modules matplotlib.pyplot and numpy.
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
We set the print options to print matrices neatly
np.set_printoptions(precision = 2) # only 2 fractionary digits
np.set_printoptions(suppress = True) # do not use exponential notation
Discrete polynomial approximation¶
Let us begin with an example. We search for a polynomial of degree $2$ that approximates the function $f(x)=\cos(x)$ but we only can use the values for the $5$ nodes
$$x_0=-1,x_1=-0.5,x_2=0,x_3=0.5,x_4=1$$The polynomial will be
$$P_2(x)=a_0+a_1x+a_2x^2$$and we have three unknowns $a_0,a_1,a_2.$
The equations are
$$ \begin{array}{lll} P_{2}(x_{0})=y_{0} & \quad & a_{0}+a_{1}x_{0}+a_{2}x_{0}^{2}=y_{0}\\ P_{2}(x_{1})=y_{1} & \quad & a_{0}+a_{1}x_{1}+a_{2}x_{1}^{2}=y_{1}\\ P_{2}(x_{2})=y_{2} & \quad & a_{0}+a_{1}x_{2}+a_{2}x_{2}^{2}=y_{2}\\ P_{2}(x_{3})=y_{3} & \quad & a_{3}+a_{1}x_{3}+a_{2}x_{3}^{2}=y_{3}\\ P_{2}(x_{4})=y_{4} & \quad & a_{0}+a_{1}x_{4}+a_{2}x_{4}^{2}=y_{4} \end{array} $$And we must solve the system
$$ \left(\begin{array}{ccccc} 1 & x_{0} & x_{0}^{2} \\ 1 & x_{1} & x_{1}^{2} \\ 1 & x_{2} & x_{2}^{2} \\ 1 & x_{3} & x_{3}^{2} \\ 1 & x_{4} & x_{4}^{2} \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)=\left(\begin{array}{c} y_{0}\\ y_{1}\\ y_{2}\\ y_{3}\\ y_{4} \end{array}\right) $$But we have more equations than unknowns and the system has no solution. But if we multiply the two sides by the transpose of the coefficient matrix:
$$ \left(\begin{array}{ccccc} 1 & 1 & 1 & 1 & 1 \\ x_{0} & x_{1} & x_{2} & x_{3}& x_{4}\\ x_{0}^{2} & x_{1}^{2} & x_{2}^{2} & x_{3}^{2} & x_{4}^{2} \end{array}\right) \left(\begin{array}{ccccc} 1 & x_{0} & x_{0}^{2} \\ 1 & x_{1} & x_{1}^{2} \\ 1 & x_{2} & x_{2}^{2} \\ 1 & x_{3} & x_{3}^{2} \\ 1 & x_{4} & x_{4}^{2} \end{array}\right) \left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)= \left(\begin{array}{ccccc} 1 & 1 & 1 & 1 & 1 \\ x_{0} & x_{1} & x_{2} & x_{3}& x_{4}\\ x_{0}^{2} & x_{1}^{2} & x_{2}^{2} & x_{3}^{2} & x_{4}^{2} \end{array}\right) \left(\begin{array}{c} y_{0}\\ y_{1}\\ y_{2}\\ y_{3}\\ y_{4} \end{array}\right) $$And know we have a consistent system of $3$ unknowns and $3$ equations:
$$ \left(\begin{array}{ccc} n+1 & \sum_{k=0}^{n}x_{k} & \sum_{k=0}^{n}x_{k}^{2}\\ \sum_{k=0}^{n}x_{k} & \sum_{k=0}^{n}x_{k}^{2} & \sum_{k=0}^{n}x_{k}^{3}\\ \sum_{k=0}^{n}x_{k}^{2} & \sum_{k=0}^{n}x_{k}^{3} & \sum_{k=0}^{n}x_{k}^{4} \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)=\left(\begin{array}{c} \sum_{k=0}^{n}y_{k}\\ \sum_{k=0}^{n}x_{k}y_{k}\\ \sum_{k=0}^{n}x_{k}^{2}y_{k} \end{array}\right) \quad \quad (1) $$Write a script that computes the approximating polynomial of degree two for the nodes $x_0=-1,x_1=-0.5,x_2=0,x_3=0.5,x_4=1$ of the function $f(x)=\cos(x)$. Follow the steps:
- Build the coefficients matrix and the right hand side matrix of system (1) and print them.
- Solve the system with the command
numpy.linalg.solveand print the solution. - Find the values of the approximating polynomial for the points in the interval $[-1,1]$ with the command
numpy.polyval. - Plot, in interval $[-1,1]$:
- The nodes.
- The polynomial.
%run Exercise1.py
%run Exercise2.py
We can get the polynomial coefficients with the command polyfit using as input the polynomial degree.
a2 = np.polyfit(x,y,4)
print(a2)
Continuous polynomial approximation¶
If, instead of knowing only the value of the function in some points, we know all the values for the points in an interval, we can perform a continuous fitting. In this case, the additions transform into integrals:
$$ \left(\begin{array}{ccc} \int_{-1}^{1}1 dx & \int_{-1}^{1}x dx & \int_{-1}^{1}x^2 dx\\ \int_{-1}^{1}x dx & \int_{-1}^{1}x^2 dx & \int_{-1}^{1}x^3 dx\\ \int_{-1}^{1}x^2 dx & \int_{-1}^{1}x^3 dx & \int_{-1}^{1}x^4 dx \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)=\left(\begin{array}{c} \int_{-1}^{1}f(x)dx\\ \int_{-1}^{1}x f(x)dx\\ \int_{-1}^{1}x^2 f(x)dx \end{array}\right)\quad \quad (2) $$Write a script that computes the degree two approximating polynomial for the function $f(x)=\cos(x)$ in the interval $[-1,1]$. Follow the steps:
- Build the coefficients matrix and the right hand side matrix of system (2) and print them. Use the function
scipy.integrate.quad. - Solve the system with the command
numpy.linalg.solveand print the solution. - Find the values of the approximating polynomial for the points in the interval $[-1,1]$ with the command
numpy.polyval. - Plot, in interval $[-1,1]$:
- The function $f(x)=\cos(x)$
- The polynomial.
- Compute the relative error $$E_r = \frac{\left\Vert \cos(x)-P(x)\right\Vert}{\left\Vert \cos(x)\right\Vert}$$ with the function
numpy.linalg.norm.
%run Exercise3.py
%run Exercise4.py
Approximation by orthogonal basis¶
We can assume that we have been using the dot product for functions
$$ \left\langle f(x),g(x)\right\rangle=\int_{-1}^ 1 f(x)g(x)dx$$for the polynomial basis $\{P_0,P_1,P_2\}=\{1,x,x^2\}$ we may write the system
$$ \left(\begin{array}{ccc} \left\langle P_0,P_0\right\rangle & \left\langle P_0,P_1\right\rangle & \left\langle P_0,P_2\right\rangle\\ \left\langle P_1,P_0\right\rangle & \left\langle P_1,P_1\right\rangle & \left\langle P_1,P_2\right\rangle\\ \left\langle P_2,P_0\right\rangle & \left\langle P_2,P_1\right\rangle & \left\langle P_2,P_2\right\rangle \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)=\left(\begin{array}{c} \left\langle P_0,f(x)\right\rangle\\ \left\langle P_1,f(x)\right\rangle\\ \left\langle P_2,f(x)\right\rangle \end{array}\right) $$We can improve the previous method making the coefficient matrix diagonal. We can achieve it using an orthogonal polynomial basis $\{L_0,L_1,L_2\}$. Then
$$ \left(\begin{array}{ccc} \left\langle L_0,L_0\right\rangle & 0 & 0\\ 0& \left\langle L_1,L_1\right\rangle & 0\\ 0 & 0 & \left\langle L_2,L_2\right\rangle \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)=\left(\begin{array}{c} \left\langle L_0,f(x)\right\rangle\\ \left\langle L_1,f(x)\right\rangle\\ \left\langle L_2,f(x)\right\rangle \end{array}\right) $$And the solution of the system would be
$$ a_0=\frac{\left\langle L_0,f(x)\right\rangle}{\left\langle L_0,L_0\right\rangle} \quad a_1=\frac{\left\langle L_1,f(x)\right\rangle}{\left\langle L_1,L_1\right\rangle} \quad a_2=\frac{\left\langle L_2,f(x)\right\rangle}{\left\langle L_2,L_2\right\rangle}\quad \quad (3) $$This orthogonal polynomials are Legendre polynomials. We can obtain them from the module scipy.special.
from scipy.special import legendre
def L(x,n):
Leg = legendre(n)
y = Leg(x)
return y
xx = np.linspace(-1,1,100)
plt.plot(xx,0*xx,'k-')
for i in range(6):
plt.plot(xx,L(xx,i),label=r'$L_'+str(i)+'$')
plt.legend(bbox_to_anchor=(1, 1), loc='upper left', ncol=1)
plt.title('Legendre polynomials')
plt.show()
If we solve the latter exercise using this orthogonal polynomials, the result will be the same but it would be expressed as
$$ P_2(x)=a_0\,L_0(x)+a_1\,L_1(x)+a_2\,L_2(x) \quad \quad (4) $$Write a script that computes the degree two approximating polynomial for the function $f(x)=\cos(x)$ in the interval $[-1,1]$. Follow the steps:
- Compute the coefficients in (3). Use the function
scipy.integrate.quadand the Legendre polynomials fromscipy.special. - Find the values of the polynomial (4) in the interval $[-1,1]$
- Plot, in interval $[-1,1]$:
- The function $f(x)=\cos(x)$
- Polynomial (4)
- Compute the relative error $$E_r = \frac{\left\Vert \cos(x)-P(x)\right\Vert}{\left\Vert \cos(x)\right\Vert}$$ with the function
numpy.linalg.norm.
%run Exercise5.py
%run Exercise6.py
Approximation with Fourier series¶
If $f$ is a periodic function with period $T$, the approximation with a basis of trigonometric functions works very well. We use the dot product
$$ \left\langle f(x),g(x)\right\rangle=\int_{\lambda}^{\lambda+T} f(x)g(x)dx\quad\quad (5)$$where the integration interval covers a full period.
The basis of trigonometric functions is
$$ \left\{ 1,\cos\left(\frac{2\pi x}{T}\right),\mathrm{sen} \left(\frac{2\pi x}{T}\right),\cos\left(2\frac{2\pi x}{T}\right),\mathrm{sen} \left(2\frac{2\pi x}{T}\right),\ldots,\cos\left(n\frac{2\pi x}{T}\right),\mathrm{sen} \left(n\frac{2\pi x}{T}\right)\right\} $$that is an orthogonal basis for dot product (5). Then, we may write the system $Ax=b$
$$ A= \left(\begin{array}{cccccc} \int_{\lambda}^{\lambda+T} 1^2 dx & 0 & 0 & \cdots & 0 & 0\\ 0& \int_{\lambda}^{\lambda+T} \cos^2\left(\frac{2\pi x}{T}\right) dx & 0& \cdots & 0 & 0\\ 0 & 0 & \int_{\lambda}^{\lambda+T} \mathrm{sen}^2 \left(\frac{2\pi x}{T}\right) dx & \cdots & 0 & 0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots\\ 0 & 0 & 0 & \cdots & \int_{\lambda}^{\lambda+T} \cos^2 \left(n\frac{2\pi x}{T}\right) dx & 0\\ 0 & 0 & 0 & \cdots & 0 & \int_{\lambda}^{\lambda+T} \mathrm{sen}^2 \left(n\frac{2\pi x}{T}\right) dx\\ \end{array}\right)\left(\begin{array}{c} \frac{a_{0}}{2}\\ a_{1}\\ b_{1}\\ \cdots\\ a_{n}\\ b_{n} \end{array}\right)= $$$$ =\left(\begin{array}{c} \int_{\lambda}^{\lambda+T} f(x) dx\\ \int_{\lambda}^{\lambda+T} f(x)\cos \left(\frac{2\pi x}{T}\right) dx\\ \int_{\lambda}^{\lambda+T} f(x)\mathrm{sen} \left(\frac{2\pi x}{T}\right) dx\\ \cdots\\ \int_{\lambda}^{\lambda+T} f(x)\cos \left(n\frac{2\pi x}{T}\right) dx\\ \int_{\lambda}^{\lambda+T} f(x)\mathrm{sen} \left(n\frac{2\pi x}{T}\right) dx \end{array}\right) $$And the function that approximates $f$ would be of the form:
$$ F(x)=\frac{a_0}{2}+a_1\cos\left(\frac{2\pi x}{T}\right)+b_1\mathrm{sen} \left(\frac{2\pi x}{T}\right)+a_2\cos\left(2\frac{2\pi x}{T}\right)+b_2\mathrm{sen} \left(2\frac{2\pi x}{T}\right)+\cdots+a_n\cos\left(n\frac{2\pi x}{T}\right)+b_n\mathrm{sen} \left(n\frac{2\pi x}{T}\right) $$And now,
$$ \frac{a_0}{2}=\frac{\int_{\lambda}^{\lambda+T} f(x) dx}{\int_{\lambda}^{\lambda+T} 1^2 dx} \quad a_1=\frac{\int_{\lambda}^{\lambda+T} f(x)\cos \left(\frac{2\pi x}{T}\right) dx}{\int_{\lambda}^{\lambda+T} \cos^2\left(\frac{2\pi x}{T}\right)dx} \quad b_1=\frac{\int_{\lambda}^{\lambda+T} f(x)\mathrm{sen} \left(\frac{2\pi x}{T}\right) dx}{\int_{\lambda}^{\lambda+T} \mathrm{sen}^2\left(\frac{2\pi x}{T}\right)dx} ,\ldots, \\ a_n=\frac{\int_{\lambda}^{\lambda+T} f(x)\cos \left(n\frac{2\pi x}{T}\right) dx}{\int_{\lambda}^{\lambda+T} \cos^2\left(n\frac{2\pi x}{T}\right)dx} \quad b_n=\frac{\int_{\lambda}^{\lambda+T} f(x)\mathrm{sen} \left(n\frac{2\pi x}{T}\right) dx}{\int_{\lambda}^{\lambda+T} \mathrm{sen}^2\left(n\frac{2\pi x}{T}\right)dx} $$As
$$ \int_{\lambda}^{\lambda+T} 1^2 dx= T, \quad \int_{\lambda}^{\lambda+T} \cos^2\left(k\frac{2\pi x}{T}\right)dx = \frac{T}{2}, \quad \int_{\lambda}^{\lambda+T} \mathrm{sen}^2\left(k\frac{2\pi x}{T}\right)dx = \frac{T}{2} $$for any positive integer $k$, the coefficients of the basis function are
$$ a_0=\frac{2}{T}\int_{\lambda}^{\lambda+T} f(x) dx, \quad a_1=\frac{2}{T}\int_{\lambda}^{\lambda+T} f(x) \cos\left(\frac{2\pi x}{T}\right) dx, \quad b_1=\frac{2}{T}\int_{\lambda}^{\lambda+T} f(x)\mathrm{sen} \left(\frac{2\pi x}{T}\right) dx, \ldots, a_n=\frac{2}{T}\int_{\lambda}^{\lambda+T} f(x)\cos \left(n\frac{2\pi x}{T}\right) dx, \quad b_n=\frac{2}{T}\int_{\lambda}^{\lambda+T} f(x)\mathrm{sen} \left(n\frac{2\pi x}{T}\right) dx. $$%run Exercise7.py
%run Exercise8.py
