Approximation¶
Contents¶
We import the modules matplotlib.pyplot, numpy and the numpy functions for polynomials.
import numpy as np
import matplotlib.pyplot as plt
import numpy.polynomial.polynomial as pol
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)=\sin(x)$ and we have values for $5$ nodes
$$x_1=0,x_2=0.5,x_3=1,x_4=1.5,x_5=2$$The polynomial would 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_{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_{0}+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}\\ P_{2}(x_{5})=y_{5} & \quad & a_{0}+a_{1}x_{5}+a_{2}x_{5}^{2}=y_{5} \end{array} $$And we must solve the system
$$ \left(\begin{array}{ccccc} 1 & x_{1} & x_{1}^{2} \\ 1 & x_{2} & x_{2}^{2} \\ 1 & x_{3} & x_{3}^{2} \\ 1 & x_{4} & x_{4}^{2} \\ 1 & x_{5} & x_{5}^{2} \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)=\left(\begin{array}{c} y_{1}\\ y_{2}\\ y_{3}\\ y_{4}\\ y_{5} \end{array}\right)\quad \quad (1) $$That is,
$$ Vp=y\quad \quad (1) $$But we have more equations than unknowns and the system has no solution. If we multiply the two sides by the transpose of the coefficient matrix:
$$ \left(\begin{array}{ccccc} 1 & 1 & 1 & 1 & 1 \\ x_{1} & x_{2} & x_{3}& x_{4} &x_{5}\\ x_{1}^{2} & x_{2}^{2} & x_{3}^{2} & x_{4}^{2} & x_{5}^{2} \end{array}\right) \left(\begin{array}{ccccc} 1 & x_{1} & x_{1}^{2} \\ 1 & x_{2} & x_{2}^{2} \\ 1 & x_{3} & x_{3}^{2} \\ 1 & x_{4} & x_{4}^{2} \\ 1 & x_{5} & x_{5}^{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_{1} & x_{2} & x_{3}& x_{4} &x_{5}\\ x_{1}^{2} & x_{2}^{2} & x_{3}^{2} & x_{4}^{2} & x_{5}^{2} \end{array}\right) \left(\begin{array}{c} y_{1}\\ y_{2}\\ y_{3}\\ y_{4}\\ y_{5} \end{array}\right) $$And now we have a consistent system of $3$ unknowns and $3$ equations:
$$ \left(\begin{array}{ccc} n & \sum\limits_{k=1}^{n}x_{k} & \sum\limits_{k=1}^{n}x_{k}^{2}\\ \sum\limits_{k=1}^{n}x_{k} & \sum\limits_{k=1}^{n}x_{k}^{2} & \sum\limits_{k=1}^{n}x_{k}^{3}\\ \sum_\limits{k=1}^{n}x_{k}^{2} & \sum\limits_{k=1}^{n}x_{k}^{3} & \sum\limits_{k=1}^{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\limits_{k=1}^{n}y_{k}\\ \sum\limits_{k=1}^{n}x_{k}y_{k}\\ \sum\limits_{k=1}^{n}x_{k}^{2}y_{k} \end{array}\right) \quad \quad (2) $$That is, we have a system
$$Cp=d \quad \quad (2)$$where the data are $C$ y $d$ and we want to calculate the coefficients of the approximation polynomial $p.$
Write a function approx1(f,deg,a,b,n) that computes and plots the approximating polynomial of degree deg for the n equispaced points of the function f in the interval [a,b] (the border points are nodes). The function will not return anything: it will do the necessary calculations to obtain the approximation polynomial and it will plot the polynomial and the points. Follow the steps:
- Build the nodes
xwithnp.linspaceand obtainy = f(x) - Build the coefficients matrix $C$ and the right hand side matrix $d$ of system (2).
- Build the matrix $V$ from system (1) and the obtain
C = np.dot(V.T,V)andd = np.dot(V.T,y), whereV.Tis the transposed matrix of $V$ or - Build $C$ directly using
np.sumthat sums the elements of a vector. Also, take into account thatx**2is a vector whose elements arexelements squared.
- Build the matrix $V$ from system (1) and the obtain
- Solve the system $Cp=d$ with the command
np.linalg.solve(C,d). - Obtain 50 points in the interval $[a,b]$ and the corresponding values of the polynomial
pwithpol.polyval. - Plot, in interval $[a,b]$:
- The points.
- The polynomial.
Use this function to obtain and plot:
- The degree $\mathrm{deg} = 2$ approximation polynomial for the function $f_1(x)=\sin(x)$ using $n = 5$ equispaced points in the interval $[0,2].$
- The degree $\mathrm{deg} = 4$ approximation polynomial for the function $f_2(x) =\cos(\arctan(x))-\log(x+5)$ using $n = 10$ equispaced points in the interval $[-2,0].$
%run Exercise1.py
We can get the polynomial coefficients with the command polyfit using as input the polynomial degree.
x = np.linspace(0,2,5)
y = np.sin(x)
p = pol.polyfit(x,y,2)
print(p)
Using the dataset Auto MPG from UC Irvine Machine Learning Repository We will estimate the power (horsepower) needed for a car of a given weight (weight) and the miles per gallon (mpg) that it would be able to travel. The car models in this database are from the year 1970 to 1982. The database is also in kaggle.
Using the pandas library, we read the data file cars.csv and store it as a dataframe, which is the equivalent of an Excel sheet: the data is organized in a matrix with one sample per row and one variable per column. The data can be non homogeneous. The dataframe can contain variables of different types (categorical, ordinal, numerical,...)
import pandas as pd
df = pd.read_csv('http://www.unioviedo.es/compnum/labs/new/cars.csv',sep=',')
Another way to read the file would be to download it to the working directory from cars and read it with
df = pd.read_csv('cars.csv',sep=',') # It reads data separated by commas
Let us have a look at the dataframe
df
For intance, to extract the column horsepower
x = df['horsepower']
Using the file cars that contains the dataset Auto MPG from UC Irvine Machine Learning Repository, estimate the necessary power (horsepower) for a car of a given weight(weight) and the miles per gallon (mpg) that would be able to travel in a city. Follow the steps:
- Extract the columns weight y horsepower.
- Using
pol.polyfitcalculate and plot the degree 1 fitting polynomial using as variable $x$ the weight and as variable $y$ the horsepower. - Using
pol.polyvalestimate the necessary power (horsepower) of a $3000$ lbs car. - Extract the column mpg and using as variable $x$ the horsepower and as variable $y$ the mpg calculate and plot the degree 2 fitting polynomial.
- With the horsepower obtained for the $3000$ lbs car, estimate the miles per gallon that this car could travel within a city using the approximation polynomial.
- Plot the data, the corresponding fitting polynomial and the point over the polynomial used to estimate the corresponding mpg.
%run Exercise2.py
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_{a}^{b}1 dx & \int_{a}^{b}x dx & \int_{a}^{b}x^2 dx\\ \int_{a}^{b}x dx & \int_{a}^{b}x^2 dx & \int_{a}^{b}x^3 dx\\ \int_{a}^{b}x^2 dx & \int_{a}^{b}x^3 dx & \int_{a}^{b}x^4 dx \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2} \end{array}\right)=\left(\begin{array}{c} \int_{a}^{b}f(x)dx\\ \int_{a}^{b}x f(x)dx\\ \int_{a}^{b}x^2 f(x)dx \end{array}\right)\quad \quad (3) $$Write a function approx2(f,deg,a,b) that computes and plots the degree deg fitting polynomial p for the function f in the interval [a,b]. Follow the steps:
- Build the coefficients matrix and the right hand side matrix of system (3) and print them. Use the function
scipy.integrate.quad. For instance, to compute the value of an integral $$I = \int_{a}^{b}g(x)dx$$wherefrom scipy.integrate import quad Ia = quad(g,a,b)[0]
gis a lambda function. We can define a lambda function from other one. For exampleg = lambda x: x * f(x)
- Solve the system with the command
np.linalg.solveand print the solution. - Find the values of the approximating polynomial for the points in the interval $[a,b]$ with the command
pol.polyval. - Plot, in interval $[a,b]$:
- The function.
- The polynomial.
Use this function to obtain and plot:
- The degree $\mathrm{deg} = 2$ fitting polynomial for the function
in the interval $[a,b]=[0,2].$
- The degree $\mathrm{deg}=4$ fitting polynomial for the function
in the interval $[a,b]=[-2,0].$
%run Exercise3.py
We can assume that we have been using the dot product for functions
$$ \left\langle f(x),g(x)\right\rangle=\int_{a}^ b 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 the Legendre polynomials. We can obtain them from the module scipy.special.
from scipy.special import eval_legendre
xp = np.linspace(-1,1,100)
plt.figure(figsize=(8,6))
plt.plot(xp,0*xp,'k-')
for i in range(6):
plt.plot(xp,eval_legendre(i,xp),label=r'$L_'+str(i)+'$')
plt.legend(bbox_to_anchor=(1, 1), loc='upper left', ncol=1,fontsize=16)
plt.title('Legendre polynomials',fontsize=16)
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) $$But this Legendre polynomials are defined for $[-1,1].$
We can apply this to any $[a,b]$ interval with a linear transformation.
We want to change variable that takes us from the interval $[-1,1]$ to the interval $[a,b].$
If we perform a linear change of variable
$$x = m\,t +n,$$that if $x = -1$, then $t = a$ then
$$-1 = a\,m+n.$$And if $x = 1$, then $t = b$ then
$$1 = b\,m+n.$$We multiply the first equation by $-1$ and we add these two equations
$$2 = (b-a)\,m \qquad m = \dfrac{2}{b-a}$$Now we substitute this value in the first equation
$$-1 = a\,\dfrac{2}{b-a} +n$$And solving this equation for $n$ we get
$$n =-\dfrac{a+b}{b-a}$$Thus, the change of variable is
$$\fbox{$x=\dfrac{2t-(a+b)}{b-a}$}$$And if we want to use the Legendre polynomials in an interval that it is not $[-1,1]$
a = -1.; b = 3.
xp = np.linspace(a,b,100)
plt.figure(figsize=(8,6))
plt.plot(xp,0*xp,'k-')
for i in range(6):
plt.plot(xp,eval_legendre(i,(2*xp-(a+b))/(b-a)),label=r'$L_'+str(i)+'$')
plt.legend(bbox_to_anchor=(1, 1), loc='upper left', ncol=1,fontsize=16)
plt.title('Legendre polynomials',fontsize=16)
plt.show()
Write a function approx3(f,deg,a,b) that computes and plots the degree deg fitting polynomial p for the function f in the interval [a,b]. Follow the steps:
- Compute the coefficients in (3). Use the function
scipy.integrate.quadand the Legendre polynomials fromscipy.specialtransformed for the interval $[a,b]$ - Find the values of the polynomial (4) in the interval $[a,b]$
- Plot, in the interval $[0,2]$:
- The function $f(x)=\cos(x)$
- The degree $2$ approximation polynomial (4)
Also compute the degree $4$ polynomial approximation of $f_2(x)=\cos(\arctan(x))-\log(x+5)$ in the interval $[-2,0]$. Plot the function and the polynomial approximation.
%run Exercise4.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 Exercise5.py
