{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Cubic spline interpolation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will approximate the function\n", "$$\n", "f(x) = (1-x^2)^2 \\sin(4\\pi x) \\exp(\\sin(2\\pi x)), \\qquad x \\in [-1,+1]\n", "$$\n", "We have taken an example for which \n", "$$\n", "f'(-1) = f'(+1) = 0\n", "$$\n", "which will be used in constructing the cubic spline below." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "%config InlineBackend.figure_format = 'svg'\n", "\n", "from numpy import linspace,sin,cos,exp,pi,polyfit,polyval,zeros\n", "from numpy.linalg import solve\n", "from matplotlib.pyplot import plot,legend\n", "\n", "xmin, xmax = -1.0, +1.0\n", "f = lambda x: (1-x**2)**2*sin(4*pi*x)*exp(sin(2*pi*x))\n", "\n", "N = 20\n", "x = cos(linspace(0,pi,N+1))\n", "y = f(x)\n", "xe = linspace(xmin,xmax,200)\n", "ye = f(xe)\n", "plot(xe,ye,x,y,'o')\n", "legend(('f(x)','Data'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us first try polynomial interpolation using Chebyshev points." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "p = polyfit(x,y,N)\n", "yp = polyval(p,xe)\n", "plot(xe,ye,x,y,'o',xe,yp,'r-')\n", "legend(('Exact','Data','Polynomial'));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that polynomial interpolation is rather oscillatory. Of course this will become better as we increase the number of samples. Using uniformly spaced points for polynomial interpolation would have given us even bad approximation. We will next try cubic spline interpolation. First we implement a function which constructs the matrix and right hand side." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def spline_matrix(x,y):\n", " N = len(x) - 1\n", " h = x[1:] - x[0:-1]\n", " A = zeros((N+1,N+1))\n", " A[0][0] = h[0]/3; A[0][1] = h[0]/6\n", " for i in range(1,N):\n", " A[i][i-1] = h[i-1]/6; A[i][i] = (h[i-1]+h[i])/3; A[i][i+1] = h[i]/6\n", " A[N][N-1] = h[N-1]/6; A[N][N] = h[N-1]/3\n", " \n", " r = zeros(N+1)\n", " r[0] = (y[1] - y[0])/h[0]\n", " for i in range(1,N):\n", " r[i] = (y[i+1]-y[i])/h[i] - (y[i]-y[i-1])/h[i-1]\n", " r[N] = -(y[N]-y[N-1])/h[N-1]\n", " return A,r" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will use uniformly spaced points for spline interpolation. We first construct the matrix and right hand side and solve the matrix problem. Finally, we evaluate the spline in each interval to make a plot." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x = linspace(xmin,xmax,N+1)\n", "y = f(x)\n", "\n", "A,r = spline_matrix(x,y)\n", "m = solve(A,r)\n", "\n", "plot(xe,ye,x,y,'o')\n", "for i in range(0,N):\n", " h = x[i+1] - x[i]\n", " xs = linspace(x[i],x[i+1],6)\n", " ys = ((x[i+1]-xs)**3*m[i] + (xs-x[i])**3*m[i+1])/(6*h) \\\n", " + ((x[i+1]-xs)*y[i] + (xs-x[i])*y[i+1])/h \\\n", " - ((x[i+1]-xs)*m[i] + (xs-x[i])*m[i+1])*h/6\n", " plot(xs,ys,'r-')\n", "\n", "legend(('Exact','Data','Spline'));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that spline gives a better approximation and does not have too much oscillation." ] } ], "metadata": { "kernelspec": { "display_name": "Python [default]", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.13" } }, "nbformat": 4, "nbformat_minor": 2 }