{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Chebyshev Interpolation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The input of the interpolation problem are $n+1$ tuples $(x_k, f_k)$, for $k = 0,1,\\ldots,n$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The output is a polynomial $p$ of degree at most $n$ so that $p(x_k) = f_k$, for $k=0,1,\\ldots,n$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If $f_k = f(x_k)$, for $k=0,1,\\ldots,n$,\n",
"and if we can choose the $x_k$ freely,\n",
"then, for a fixed value of $n$, \n",
"how should we choose $x_k$ to obtain the best approximation for the function $f$?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We consider the sine function $f(x) = \\sin(x)$ over the interval $[-\\pi,\\pi]$ and we apply Chebyshev interpolation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Chebyshev Points"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us take six points in the interval $[-\\pi, +\\pi]$. The points are roots of the Chebyshev polynomials, which are given by explicit formulas."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"using Printf\n",
"Base.show(io::IO, f::Float64) = @printf(io, \"%.3e\", f)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6-element Vector{Float64}:\n",
" 9.659e-01\n",
" 7.071e-01\n",
" 2.588e-01\n",
" -2.588e-01\n",
" -7.071e-01\n",
" -9.659e-01"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 6\n",
"r = [cos((2*i-1)*pi/(2*n)) for i=1:n]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Observe that the roots in ``r`` all lie in the interval $[-1,+1]$, but we want to interpolate the $\\sin(x)$ for $x \\in [-\\pi,\\pi]$. Therefore, we must perform a coordinate transformation and map the points in $[-1,+1]$ into $[-\\pi,\\pi]$. For this example, this coordinate transformation happens with a simple multiplication of all points with $\\pi$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6-element Vector{Float64}:\n",
" 3.035e+00\n",
" 2.221e+00\n",
" 8.131e-01\n",
" -8.131e-01\n",
" -2.221e+00\n",
" -3.035e+00"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X = [pi*z for z in r]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The corresponding function values evaluate the sine function at ``X``."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6-element Vector{Float64}:\n",
" 1.068e-01\n",
" 7.957e-01\n",
" 7.264e-01\n",
" -7.264e-01\n",
" -7.957e-01\n",
" -1.068e-01"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"F = [sin(x) for x in X]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now apply Newton interpolation to construct the interpolating polynomial through those points."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. Newton Interpolation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To compute the divided differences, we apply the function ``divdif``."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"divdif"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" divdif(x::Array{Float64,1},f::Array{Float64,1})\n",
"\n",
"Returns the vector of divided differences for the\n",
"Newton form of the interpolating polynomial.\n",
"\n",
"On entry are x and f;\n",
"x contains the abscisses, given as a column vector;\n",
"f contains the ordinates, given as a column vector.\n",
"\n",
"On return are the divided differences.\n",
"\"\"\"\n",
"function divdif(x::Array{Float64,1},f::Array{Float64,1})\n",
" n = length(x)\n",
" d = deepcopy(f)\n",
" for i=2:n\n",
" for j=1:i-1\n",
" d[i] = (d[j] - d[i])/(x[j] - x[i])\n",
" end\n",
" end\n",
" return d\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6-element Vector{Float64}:\n",
" 1.068e-01\n",
" -8.472e-01\n",
" -4.035e-01\n",
" -3.257e-02\n",
" 1.763e-02\n",
" 5.810e-03"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"c = divdif(X,F)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The interpolating polynomial is defined as a function."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"p"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" p(x)\n",
"\n",
"returns the value of the interpolating polynomial\n",
"through the points X with corresponding function values in F,\n",
"with as its coefficients the divided differences in c.\n",
"\"\"\"\n",
"function p(x)\n",
" n = length(c)\n",
" result = c[n]\n",
" for i=n-1:-1:1\n",
" result = result*(x - X[i]) + c[i]\n",
" end\n",
" return result\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us verify that the interpolation conditions are met by evaluating ``p`` at the points in ``X``."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6-element Vector{Float64}:\n",
" 1.068e-01\n",
" 7.957e-01\n",
" 7.264e-01\n",
" -7.264e-01\n",
" -7.957e-01\n",
" -1.068e-01"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"PX = [p(xk) for xk in X]"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4.119e-16"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using LinearAlgebra\n",
"norm(PX - F)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Indeed, the polynomial $p(x)$ has identical values as $f(x) = \\sin(x)$ at the interpolation points."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. Some Plots"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"(process:42676): GLib-GIO-WARNING **: 20:17:43.172: Unexpectedly, UWP app `Evernote.Evernote_10.45.18.0_x64__q4d96b2w5wcc2' (AUMId `Evernote.Evernote_q4d96b2w5wcc2!Evernote') supports 1 extensions but has no verbs\n"
]
}
],
"source": [
"using Plots"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(sin, -0.1-pi:0.05:pi+0.1, label=\"sin\", legend=:topleft)\n",
"plot!(p, -0.1-pi:0.05:pi+0.1, label=\"p\")\n",
"plot!(X, F, line=(:scatter), label=\"\")"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# savefig(\"plotsinCheb6.png\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us plot the interpolation error."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"e (generic function with 1 method)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"e(x) = sin(x) - p(x)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(e, -0.025-pi:0.05:pi+0.025, label=\"sin-p\", legend=:topleft)\n",
"plot!(X, zeros(6), line=(:scatter), label=\"\")"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# savefig(\"plotsinCheb6error.png\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3. Rational Approximations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The construction of Padé approximations start a Taylor series, for which we use ``SymPy``."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$x$"
],
"text/plain": [
"x"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using SymPy\n",
"x = Sym(\"x\")"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$e^{x}$"
],
"text/plain": [
" x\n",
"e "
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = exp(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"By evaluating the ``f`` at a ``SymPy`` symbol, we obtain a symbolic expression that we then can feed as input to ``series``."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Sym"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"typeof(f)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For a rational approximation with degree 4 in numerator and denominator we need 9 coefficients, 5 for the numerator and 4 for the denominator (because the constant coefficient of the denominator is one)."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$1 + x + \\frac{x^{2}}{2} + \\frac{x^{3}}{6} + \\frac{x^{4}}{24} + \\frac{x^{5}}{120} + \\frac{x^{6}}{720} + \\frac{x^{7}}{5040} + O\\left(x^{8}\\right)$"
],
"text/plain": [
" 2 3 4 5 6 7 \n",
" x x x x x x / 8\\\n",
"1 + x + -- + -- + -- + --- + --- + ---- + O\\x /\n",
" 2 6 24 120 720 5040 "
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ts = series(f, x0=0, n=8)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"9-element Vector{Float64}:\n",
" 1.000e+00\n",
" 1.000e+00\n",
" 5.000e-01\n",
" 1.667e-01\n",
" 4.167e-02\n",
" 8.333e-03\n",
" 1.389e-03\n",
" 1.984e-04\n",
" 0.000e+00"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"c = [Float64(ts.coeff(x,k)) for k=0:8]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The coefficients of the denominator are the solution of a 4-by-4 linear system."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Matrix{Float64}:\n",
" 1.667e-01 5.000e-01 1.000e+00 1.000e+00\n",
" 4.167e-02 1.667e-01 5.000e-01 1.000e+00\n",
" 8.333e-03 4.167e-02 1.667e-01 5.000e-01\n",
" 1.389e-03 8.333e-03 4.167e-02 1.667e-01"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = zeros(4,4)\n",
"A[1,:] = c[4:-1:1]\n",
"A[2,:] = c[5:-1:2]\n",
"A[3,:] = c[6:-1:3]\n",
"A[4,:] = c[7:-1:4]\n",
"A"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Vector{Float64}:\n",
" -8.333e-03\n",
" -1.389e-03\n",
" -1.984e-04\n",
" -0.000e+00"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rhs = [-c[6], -c[7], -c[8], -c[9]]"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Vector{Float64}:\n",
" -3.571e-01\n",
" 1.667e-01\n",
" -3.571e-02\n",
" 3.571e-03"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"b = A\\rhs"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Vector{Float64}:\n",
" 1.000e+00\n",
" 6.429e-01\n",
" 3.095e-01\n",
" 1.190e-01\n",
" 3.333e-02"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = zeros(5)\n",
"a[1] = c[1]\n",
"a[2] = c[2] + b[1]*c[1]\n",
"a[3] = c[3] + b[1]*c[2] + b[2]*c[1]\n",
"a[4] = c[4] + b[1]*c[3] + b[2]*c[2] + b[3]*c[1]\n",
"a[5] = c[5] + b[1]*c[4] + b[2]*c[3] + b[3]*c[2] + b[4]*c[1]\n",
"a"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To check the accuracy of the approximation, we use ``Polynomials`` to define and plot the rational approximation."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"using Polynomials"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"1.000e+00 + 6.429e-01∙x + 3.095e-01∙x2 + 1.190e-01∙x3 + 3.333e-02∙x4"
],
"text/latex": [
"$1.000e+00 + 6.429e-01\\cdot x + 3.095e-01\\cdot x^{2} + 1.190e-01\\cdot x^{3} + 3.333e-02\\cdot x^{4}$"
],
"text/plain": [
"Polynomial(1.000e+00 + 6.429e-01*x + 3.095e-01*x^2 + 1.190e-01*x^3 + 3.333e-02*x^4)"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"numerator = Polynomial(a)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the denominator, we should not forget that the leading coefficient equals one, so we insert ``1.0`` to the coefficient vector ``b``."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Vector{Float64}:\n",
" 1.000e+00\n",
" -3.571e-01\n",
" 1.667e-01\n",
" -3.571e-02\n",
" 3.571e-03"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bp1 = vcat([1.0], b)"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"1.000e+00 - 3.571e-01∙x + 1.667e-01∙x2 - 3.571e-02∙x3 + 3.571e-03∙x4"
],
"text/latex": [
"$1.000e+00 - 3.571e-01\\cdot x + 1.667e-01\\cdot x^{2} - 3.571e-02\\cdot x^{3} + 3.571e-03\\cdot x^{4}$"
],
"text/plain": [
"Polynomial(1.000e+00 - 3.571e-01*x + 1.667e-01*x^2 - 3.571e-02*x^3 + 3.571e-03*x^4)"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"denominator = Polynomial(bp1)"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"q (generic function with 1 method)"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"q(z) = numerator(z)/denominator(z)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"range = -1:0.01:1\n",
"plot(f, range, label=\"f\", legend=:topleft)\n",
"plot!(q, range, label=\"q\", legend=:topleft)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, let us look at the error."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"e(z) = f(z) - q(z)\n",
"plot(e, range, label=\"error\", legend=:topleft)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The shape of the error is such that the error is smallest close to the origin, similar to the Taylor series, which is accurate around ``0.0``. That the error is larger at ``+1.0`` than at ``-1.0`` is because the function grows faster at the right than at the left of the interval ``[-1.0, +1.0]``. "
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.8.0",
"language": "julia",
"name": "julia-1.8"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.0"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"autocomplete": true,
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 1,
"hotkeys": {
"equation": "Ctrl-E",
"itemize": "Ctrl-I"
},
"labels_anchors": false,
"latex_user_defs": false,
"report_style_numbering": false,
"user_envs_cfg": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}