{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The Interpolation Problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The interpolation problem is illustrated with a simple example that is then solved with the Vandermonde matrix."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider the sine function $f(x) = \\sin(x)$ over the interval $[-\\pi,\\pi]$\n",
"and the points $x_0 = -\\pi$, $x_1 = -\\pi/2$, $x_2 = 0$, $x_3 = \\pi/2$, $x_4 = \\pi$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compute the polynomial $p(x)$ interpolating $f(x)$ at $x_k$: $p(x_k) = f(x_k)$, for $k=0,1,2,3,4$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Solving a Linear System"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To solve this simple problem, a linear system will be solved. The coefficient matrix is the Vandermonde matrix defined by the points $x_k$ and the right hand side vector is defined by the function values $f(x_k)$, for $k=1,2,3,4,5$."
]
},
{
"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": [],
"source": [
"using LinearAlgebra"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Below, the Vandermonde matrix ``V`` is defined for the five points."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5×5 Matrix{Float64}:\n",
" 1.000e+00 -3.142e+00 9.870e+00 -3.101e+01 9.741e+01\n",
" 1.000e+00 -1.571e+00 2.467e+00 -3.876e+00 6.088e+00\n",
" 1.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00\n",
" 1.000e+00 1.571e+00 2.467e+00 3.876e+00 6.088e+00\n",
" 1.000e+00 3.142e+00 9.870e+00 3.101e+01 9.741e+01"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X = [-pi, -pi/2, 0.0, pi/2, pi]\n",
"V = zeros(5,5)\n",
"for i=1:5\n",
" V[i,1] = 1.0\n",
" V[i,2] = X[i]\n",
" for j=3:5\n",
" V[i,j] = V[i,j-1]*X[i]\n",
" end\n",
"end\n",
"V"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us now define the right hand side vector ``F``."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Vector{Float64}:\n",
" -1.225e-16\n",
" -1.000e+00\n",
" 0.000e+00\n",
" 1.000e+00\n",
" 1.225e-16"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"F = [sin(xk) for xk in X]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we solve the system $V c = F$ for the coefficients $c$ of the interpolating polynomial."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Vector{Float64}:\n",
" 0.000e+00\n",
" 8.488e-01\n",
" 5.999e-17\n",
" -8.600e-02\n",
" -6.079e-18"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"c = V\\F"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The interpolating polynomial is defined as a function."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"p (generic function with 1 method)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p(x) = sum(c[i]*x^(i-1) for i=1:5)"
]
},
{
"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": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Vector{Float64}:\n",
" -1.480e-16\n",
" -1.000e+00\n",
" 0.000e+00\n",
" 1.000e+00\n",
" -1.480e-16"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"PX = [p(xk) for xk in X]"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5.206e-16"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"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. A Plot"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"(process:85984): GLib-GIO-WARNING **: 19:12:53.286: Unexpectedly, UWP app `Evernote.Evernote_10.44.8.0_x64__q4d96b2w5wcc2' (AUMId `Evernote.Evernote_q4d96b2w5wcc2!Evernote') supports 1 extensions but has no verbs\n"
]
}
],
"source": [
"using Plots"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(sin, -0.05-pi:0.05:pi+0.1, label=\"sin\", legend=:topleft)\n",
"plot!(p, -0.05-pi:0.05:pi+0.1, label=\"p\")\n",
"plot!(X, F, line=(:scatter), label=\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3. Lagrange polynomials"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Lagrange polynomials interpolate 0/1 functions."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"using Polynomials"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"lagrange"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" lagrange(pts::Array{Float64,1},idx::Int64)\n",
"\n",
"Returns the Lagrange polynomial for the points in pts\n",
"with index idx.\n",
"\n",
"REQUIRED:\n",
"\n",
" All points in pts must be distinct.\n",
" The index idx is between 1 and length(pts).\n",
"\"\"\"\n",
"function lagrange(pts::Array{Float64,1},idx::Int64)\n",
" result = 1.0\n",
" for i=1:idx-1\n",
" result = result*Polynomial([-pts[i], 1.0])/\n",
" (-pts[i]+pts[idx])\n",
" end\n",
" for i=idx+1:length(pts)\n",
" result = result*Polynomial([-pts[i], 1.0])/\n",
" (-pts[i]+pts[idx])\n",
" end\n",
" return result\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us verify that, for a random collection of three points,\n",
"the value of the Lagrange polynomials are 0 or 1,\n",
"and that the sum of the Lagrange polynomials equals one."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" 2.404e-01\n",
" 4.741e-02\n",
" 7.554e-01"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pts = rand(3)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Lagrange polynomial L1:\n",
"-3.603e-01 + 8.077e+00*x - 1.006e+01*x^2\n",
"L1(pts[1]) = 1.000e+00\n",
"L1(pts[2]) = 6.493e-17\n",
"L1(pts[3]) = 1.243e-15\n",
"Lagrange polynomial L2:\n",
"1.329e+00 - 7.287e+00*x + 7.318e+00*x^2\n",
"L2(pts[1]) = -2.160e-16\n",
"L2(pts[2]) = 1.000e+00\n",
"L2(pts[3]) = -3.588e-16\n",
"Lagrange polynomial L3:\n",
"3.126e-02 - 7.894e-01*x + 2.743e+00*x^2\n",
"L3(pts[1]) = 1.068e-17\n",
"L3(pts[2]) = -2.373e-18\n",
"L3(pts[3]) = 1.000e+00\n",
"Sum of Lagrange polynomials :\n",
"1.000e+00 + 8.882e-16*x + 8.882e-16*x^2\n"
]
}
],
"source": [
"sumLagrange = 0.0\n",
"for i=1:3\n",
" Li = lagrange(pts, i)\n",
" println(\"Lagrange polynomial L\", i, \":\")\n",
" println(Li)\n",
" for j=1:3\n",
" print(\"L\", i, \"(pts[\", j, \"]) = \")\n",
" println(Li(pts[j]))\n",
" end\n",
" sumLagrange = sumLagrange + Li\n",
"end\n",
"println(\"Sum of Lagrange polynomials :\")\n",
"println(sumLagrange)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3. Lagrange Interpolation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider the following data in ``x`` and ``f`` below."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Vector{Float64}:\n",
" 3.200e+01\n",
" 2.220e+01\n",
" 4.160e+01\n",
" 1.010e+01\n",
" 5.050e+01"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = [32.0, 22.2, 41.6, 10.1, 50.5]"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Vector{Float64}:\n",
" 5.299e-01\n",
" 3.778e-01\n",
" 6.639e-01\n",
" 1.754e-01\n",
" 6.361e-01"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = [0.52992, 0.37784, 0.66393, 0.17537, 0.63608]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given $(x_i, f_i)$, for $i=0,1,\\ldots, n$, the Lagrange interpolating polynomial is\n",
"\n",
"$$\n",
" P(x) = \\sum_{i=0}^n \\ell_i(x) f_i, \n",
"$$\n",
"\n",
"where $\\ell_i$ is the $i$-th Lagrange polynomial."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"-2.147e-01 + 6.026e-02∙x - 2.829e-03∙x2 + 7.525e-05∙x3 - 7.181e-07∙x4"
],
"text/latex": [
"$-2.147e-01 + 6.026e-02\\cdot x - 2.829e-03\\cdot x^{2} + 7.525e-05\\cdot x^{3} - 7.181e-07\\cdot x^{4}$"
],
"text/plain": [
"Polynomial(-2.147e-01 + 6.026e-02*x - 2.829e-03*x^2 + 7.525e-05*x^3 - 7.181e-07*x^4)"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P = sum([lagrange(x, i)*f[i] for i=1:5])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us verify the interpolation conditions, given by ``x`` and ``f``."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Vector{Float64}:\n",
" 5.299e-01\n",
" 3.778e-01\n",
" 6.639e-01\n",
" 1.754e-01\n",
" 6.361e-01"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Vector{Float64}:\n",
" 5.299e-01\n",
" 3.778e-01\n",
" 6.639e-01\n",
" 1.754e-01\n",
" 6.361e-01"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Px = [P(z) for z in x]"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.722e-14"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"norm(Px - f)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Indeed, at each point ``z`` in the vector ``x``, the polynomial ``P`` takes the same value as the corresponding element in ``f``."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 4. Neville Interpolation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Neville interpolation solves the value problem, computing the value of the interpolating polynomial, given the interpolation data and one point."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"neville"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
" neville(x::Array{Float64,1},f::Array{Float64,1},xx::Float64)\n",
"\n",
"Implements the interpolation algorithm of Neville.\n",
"\n",
"ON ENTRY :\n",
"\n",
" x are the abscisses, given as a column vector\n",
"\n",
" f are the ordinates, given as a column vector\n",
"\n",
" xx is the point where to evaluate the interpolating\n",
" polynomial through (x[i],f[i])\n",
"\n",
"ON RETURN :\n",
" p is the last row of Neville's table where p[1] is\n",
" the value of the interpolator at xx \n",
"\n",
"EXAMPLE :\n",
"\n",
" x = [32.0, 22.2, 41.6, 10.1, 50.5]; \n",
" f = [0.52992, 0.37784, 0.66393, 0.17537, 0.63608];\n",
" xx = 27.5;\n",
" p = neville(x,f,xx)\n",
"\"\"\"\n",
"function neville(x::Array{Float64,1},f::Array{Float64,1},xx::Float64)\n",
" n = length(x)\n",
" p = f\n",
" dx = [0.0 for i=1:n]\n",
" for i=1:n\n",
" dx[i] = xx - x[i]\n",
" end\n",
" for i=2:n\n",
" for j=2:i\n",
" p[i-j+1] = (dx[i]*p[i-j+1] - dx[i-j+1]*p[i-j+2])/\n",
" (x[i-j+1] - x[i])\n",
" end\n",
" end\n",
" return p\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Vector{Float64}:\n",
" 4.575e-01\n",
" 4.790e-01\n",
" 5.584e-01\n",
" 3.738e-01\n",
" 6.361e-01"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = [32.0, 22.2, 41.6, 10.1, 50.5]\n",
"f = [0.52992, 0.37784, 0.66393, 0.17537, 0.63608]\n",
"xx = 27.5\n",
"neville(x,f,xx)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us verify with the Lagrange interpolating polynomial ``P`` computed earlier."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4.575e-01"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P(27.5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Indeed, we see that the function value of ``P`` at ``xx = 27.5`` agrees with the first element returned by the function ``neville``."
]
}
],
"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"
}
},
"nbformat": 4,
"nbformat_minor": 4
}