{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Tutorial 05 - Function approximation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Piece-wise linear interpolation, Lagrange interpolation and Neville’s algorithm, Chebyshev polynomials and approximation, polynomial least squares fit."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy import linalg, interpolate, special\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Local interpolation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Linear interpolation\n",
"\n",
"[Linear interpolation](https://en.wikipedia.org/wiki/Linear_interpolation) is defined as a concatenation of straight lines connecting each pair of the set of known points. If the two known points are given by the coordinates $ (x_{i}, y_{i}) $ and $ (x_{i + 1}, y_{i + 1}) $, the linear interpolation in the interval $ (x_{i}, x_{i + 1}) $ is given by the following formula,\n",
"\n",
"$$\n",
"y = y_i + \\alpha (x - x_i), \\qquad x \\in [x_i, x_{i+1}],\n",
"$$\n",
"\n",
"where\n",
"\n",
"$$\n",
"\\alpha = \\frac{y_{i+1} - y_{i}}{x_{i + 1} - x_{i}}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
" \n",
"**Exercise 05.1:** Implement a function that calculates the linear interpolation.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def linear_interpolation(x_p, y_p, x):\n",
" \"\"\"\n",
" Calculates the linear interpolation.\n",
" Args:\n",
" x_p (array_like): X-coordinates of a set of datapoints\n",
" y_p (array_like): Y-coordinates of a set of datapoints\n",
" x (array_like): An array on which the interpolation is calculated\n",
" Returns:\n",
" numpy.ndarray: The linear interpolation\n",
" \"\"\" \n",
"\n",
" # add your code here"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.interpolate.interp1d`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x_p = np.random.rand(10)\n",
"y_p = np.random.rand(10)\n",
"x = np.linspace(np.min(x_p), np.max(x_p), 1000)\n",
"try:\n",
" np.testing.assert_array_almost_equal(linear_interpolation(x_p, y_p, x), interpolate.interp1d(x_p, y_p)(x), decimal=7)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 05.2:** Approximate the $ \\sin x $ function with the linear interpolation using $ 8 $ uniformly spaced points $ x_i $ from the $ [-\\pi, \\ \\pi] $ interval, calculate the absolute error of the approximation, and plot the results. \n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"# add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Global interpolation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Lagrange interpolation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given a set of $ n $ known data points $ (x_{i}, y_{i}) $, $ i = 0, 1, \\dots, n-1 $, where no two $ x_i $ are the same, the [Lagrange interpolation](https://en.wikipedia.org/wiki/Lagrange_polynomial) is the $ (n-1) $-th degree polynomial given by a linear combination\n",
"\n",
"$$\n",
"L(x) = \\sum_{i = 0}^{n-1} y_i F_i(x)\n",
"$$\n",
"\n",
"of Lagrange basis polynomials\n",
"\n",
"$$\n",
"F_i(x) = \\prod_{\\substack{j = 0 \\\\ j \\neq i}}^{n-1} \\frac{x - x_j}{x_i - x_j}.\n",
"$$\n",
"\n",
"The resulting Lagrange interpolating polynomial is unique."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 05.3:** Implement a function that calculates the Lagrange interpolating polynomial.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def lagrange_interpolation(x, y):\n",
" \"\"\"\n",
" Calculates a Lagrange interpolating polynomial.\n",
" Args:\n",
" x (array_like): X-coordinates of a set of datapoints\n",
" y (array_like): Y-coordinates of a set of datapoints\n",
" Returns:\n",
" numpy.poly1d: The Lagrange interpolating polynomial\n",
" \"\"\"\n",
"\n",
" # add your code here"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.interpolate.lagrange`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.lagrange.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x_p = np.random.rand(10)\n",
"y_p = np.random.rand(10)\n",
"try:\n",
" np.testing.assert_array_almost_equal(lagrange_interpolation(x_p, y_p), interpolate.lagrange(x_p, y_p), decimal=7)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 05.4:** Approximate the $ \\sin x $ function with the Lagrange interpolating polynomial using $ 8 $ uniformly spaced points $ x_i $ from the $ [-\\pi, \\ \\pi] $ interval, calculate the absolute error of the approximation, and plot the results. \n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"\n",
"# add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Neville's algorithm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Neville's algorithm](https://en.wikipedia.org/wiki/Neville%27s_algorithm) is used for recursive evaluation of Lagrange interpolating polynomial. The recurence is given by the following relation,\n",
"\n",
"$$\n",
"L_{i, j} = \\frac{(x - x_j) L_{i, j-1} - (x - x_i) L_{i+1, j}}{x_i - x_j}, \\qquad L_{i, i} = y_i, \\qquad i, j = 0, 1, \\dots, n-1.\n",
"$$\n",
"\n",
"The Lagrange interpolating polynomial is then $ L(x) = L_{0, n-1} $."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 05.5:** Implement a function that calculates the Lagrange interpolating polynomial using the Neville's algorithm.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def neville_algorithm(x, y):\n",
" \"\"\"\n",
" Calculates a Lagrange interpolating polynomial using Neville's algorithm.\n",
" Args:\n",
" x (array_like): X-coordinates of a set of datapoints\n",
" y (array_like): Y-coordinates of a set of datapoints\n",
" Returns:\n",
" numpy.poly1d: The Lagrange interpolating polynomial\n",
" \"\"\"\n",
"\n",
" # add your code here"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.interpolate.lagrange`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.lagrange.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x_p = np.random.rand(10)\n",
"y_p = np.random.rand(10)\n",
"try:\n",
" np.testing.assert_array_almost_equal(neville_algorithm(x_p, y_p), interpolate.lagrange(x_p, y_p), decimal=7)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 05.6:** Approximate the Runge's function \n",
" \n",
"$$\n",
"f(x) = \\frac{1}{1 + 25x^2}\n",
"$$\n",
"\n",
"with the Lagrange interpolating polynomial using $ 12 $ uniformly spaced points $ x_i $ from the $ [-1, \\ 1] $ interval, calculate the absolute error of the approximation, and plot the results.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"# add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Note:** The problem of oscillation at the edges of an interval which occurs in Exercise 05.6 is called [Runge's phenomenon](https://en.wikipedia.org/wiki/Runge%27s_phenomenon).\n",
" \n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Chebyshev interpolation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 05.7:** Write a function that calculates a [Chebyshev polynomial](https://en.wikipedia.org/wiki/Chebyshev_polynomials) of degree $ n $ using the following recursive formula,\n",
"\n",
"$$\n",
"T_{n}(x) = 2 x T_{n-1}(x) - T_{n-2}(x), \\quad T_0(x) = 1, \\ T_1(x) = x.\n",
"$$\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def chebyshev_polynomial(n):\n",
" \"\"\"\n",
" Calculates a Chebyshev polynomial of degree n using recursive formula.\n",
" Args:\n",
" n (int): Degree of the polynomial\n",
" Returns:\n",
" numpy.poly1d: The Chebyshev polynomial of degree n\n",
" \"\"\"\n",
"\n",
" # add your code here"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.special.eval_chebyt`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.eval_chebyt.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"n = 10\n",
"x = np.linspace(-1.0, 1.0, 1000)\n",
"try:\n",
" np.testing.assert_array_almost_equal(chebyshev_polynomial(n)(x), special.eval_chebyt(n, x), decimal=7)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 05.8:** Plot the first $ 5 $ Chebyshev polynomials on $ x \\in [-1, 1] $.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"\n",
"# add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 05.9:** Write a function that returns roots of a Chebyshev polynomial of degree $ n $.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def chebyshev_roots(n):\n",
" \"\"\"\n",
" Calculates roots of a Chebyshev polynomial of degree n.\n",
" Args:\n",
" n (int): Degree of the polynomial\n",
" Returns:\n",
" numpy.ndarray: Roots of a Chebyshev polynomial of degree n\n",
" \"\"\"\n",
"\n",
" # add your code here"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`scipy.special.roots_chebyt`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.roots_chebyt.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"n = 10\n",
"try:\n",
" np.testing.assert_array_almost_equal(chebyshev_roots(n), special.roots_chebyt(n)[0], decimal=7)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 05.10:** Approximate the Runge's function \n",
" \n",
"$$\n",
"f(x) = \\frac{1}{1 + 25x^2}\n",
"$$\n",
"\n",
"in the $[ -1, \\ 1] $ interval with the Lagrange interpolating polynomial using the roots of $ 12 $-th degree Chebyshev polynomial, calculate the absolute error of the approximation, and plot the results.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"\n",
"# add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Note:** If the interpolation interval $ [a, b] \\neq [-1, 1] $, then the nodes $ x_i $ of the Chebyshev interpolation can be found using affine transformation\n",
"\n",
"$$\n",
"x_i = \\frac{a + b}{2} + \\frac{b - a}{2} z_i, \\qquad i = 0, 1, \\dots, n - 1.\n",
"$$\n",
"\n",
"where $ z_i $ is the $ i $-th root of the $ n $-th degree Chebyshev polynomial.\n",
" \n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Curve fitting"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Least squares approximation\n",
"\n",
"Given a set of $ n $ known data points $ (x_{i}, y_{i}) $, $ i = 0, 1, \\dots, n-1 $, the [least squares method](https://en.wikipedia.org/wiki/Least_squares) finds an approximation to the data by minimizing the sum of the squares of the residuals,\n",
"\n",
"$$\n",
"S = \\sum_{i = 0}^{n - 1} r_i^2, \\qquad r_i = y_i - f(x_i; \\beta),\n",
"$$\n",
"\n",
"where $ f(x_i; \\beta) $ is the model function and $ \\beta $ is the vector of $ m $ adjustable parameters. The minimum of $ S $ is found by setting its gradient to $ 0 $,\n",
"\n",
"$$\n",
"\\frac{\\partial S}{\\partial \\beta_j} = \\sum_{i = 0}^{n - 1} 2 r_i \\frac{\\partial r_i}{\\partial \\beta_j} = \\sum_{i = 0}^{n - 1} -2 \\left( y_i - f(x_i; \\beta) \\right) \\frac{\\partial f(x_i; \\beta)}{\\partial \\beta_j} = 0, \\qquad j = 0, 1, \\dots, m - 1. \n",
"$$\n",
"\n",
"**Polynomial least squares fit:**\n",
"\n",
"In the case of polynomial least squares fit, the model function can be written in the following form,\n",
"\n",
"$$\n",
"f(x; \\beta) = \\beta_{m-1} x^{m-1} + \\beta_{m-2} x^{m-2} + \\dots + \\beta_1 x + \\beta_0.\n",
"$$\n",
"\n",
"The values of vector $ \\beta $ that minimize $ S $ can be thus obtained by solving a system of $ m $ linear algebraic equations, \n",
"\n",
"$$\n",
"\\sum_{i = 0}^{n - 1} \\left( y_i - \\beta_{m-1} x^{m-1} - \\beta_{m-2} x^{m-2} - \\dots - \\beta_1 x - \\beta_0 \\right) x_i^j = 0, \\qquad j = 0, 1, \\dots, m - 1.\n",
"$$\n",
"\n",
"The system of equations can be expressed as \n",
"\n",
"$$\n",
"\\beta_0 \\sum_{i = 0}^{n - 1} x_i^j + \\beta_1 \\sum_{i = 0}^{n - 1} x_i^{j+1} + \\dots + \\sum_{i = 0}^{n - 1} \\beta_{m-2} x_i^{j + m - 2} + \\beta_{m-1} \\sum_{i = 0}^{n - 1} x_i^{j + m-1} = \\sum_{i = 0}^{n - 1} y_i x_i^j \n",
"$$\n",
"\n",
"and in the matrix form,\n",
"\n",
"$$\n",
"\\begin{pmatrix}\n",
"\\sum 1 & \\sum x_i & \\cdots & \\sum x_i^{m-2} & \\sum x_i^{m-1} \\\\\n",
"\\sum x_i & \\sum x_i^2 & \\cdots & \\sum x_i^{m-1} & \\sum x_i^{m} \\\\\n",
"\\vdots & \\vdots & \\ddots & \\vdots & \\vdots \\\\\n",
"\\sum x_i^{m-2} & \\sum x_i^{m-1} & \\cdots & \\sum x_i^{2m-4} & \\sum x_i^{2m-3} \\\\\n",
"\\sum x_i^{m-1} & \\sum x_i^{m} & \\cdots & \\sum x_i^{2m-3} & \\sum x_i^{2m-2} \\\\\n",
"\\end{pmatrix} \n",
"\\begin{pmatrix}\n",
"\\beta_0 \\\\\n",
"\\beta_1 \\\\\n",
"\\vdots\\\\\n",
"\\beta_{m-2} \\\\\n",
"\\beta_{m-1}\n",
"\\end{pmatrix} \n",
"=\n",
"\\begin{pmatrix}\n",
"\\sum y_i \\\\\n",
"\\sum y_i x_i \\\\\n",
"\\vdots\\\\\n",
"\\sum y_i x_i^{m-2} \\\\\n",
"\\sum y_i x_i^{m-1}\n",
"\\end{pmatrix}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 05.11:** Implement the n-th degree polynomial least squares approximation.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def polynomial_least_squares(x, y, n):\n",
" \"\"\"\n",
" Calculates the n-th degree polynomial least squares approximation.\n",
" Args:\n",
" x (array_like): X-coordinates of a set of datapoints\n",
" y (array_like): Y-coordinates of a set of datapoints\n",
" n (int): degree of the approximating polynomial\n",
" Returns:\n",
" numpy.poly1d: The n-th degree polynomial least squares approximation\n",
" \"\"\"\n",
"\n",
" # add your code here"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may evaluate the correctness of your implementation using the [`numpy.polyfit`](https://numpy.org/doc/stable/reference/generated/numpy.polyfit.html) function."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x_p = np.random.rand(10)\n",
"y_p = np.random.rand(10)\n",
"n = 1\n",
"try:\n",
" np.testing.assert_array_almost_equal(polynomial_least_squares(x_p, y_p, n), np.polyfit(x_p, y_p, n), decimal=7)\n",
"except AssertionError as E:\n",
" print(E)\n",
"else:\n",
" print(\"The implementation is correct.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 05.12:** Consider a set of $ 100 $ random points distributed along a linear function with a Gaussian noise. Fit the data using the 1st degree polynomial least squares approximation.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"\n",
"# add your code here\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"**Exercise 05.13:** Consider a set of $ 100 $ random points distributed along a quadratic function with a Gaussian noise. Fit the data using the 2nd degree polynomial least squares approximation.\n",
" \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"\n",
"# add your code here\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}