{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Start notebook with --profile=sympy flag.\n",
      "\n",
      "Following the 1968 paper by Spath. Consider special case with equally spaced abscissae"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%qtconsole"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "(A, B, C, D)=symbols('A:D') # coefficients\n",
      "p = symbols('p') # tension parameter\n",
      "y1p = symbols('y1p') # derivative a left endpoint\n",
      "ynp = symbols('ynp') # derivative a right endpoint\n",
      "xcp = symbols('xcp') # x control point\n",
      "ycp = symbols('ycp') # y control point"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is the $k^{th}$ piece of the $n$-piece interpolant, \n",
      "\n",
      "$$f_k(x)  = A_k+B_k (x-x_k) + C_k \\exp(p_k (x-x_k)) +D_k \\exp(-p_k (x-x_k))$$\n",
      "\n",
      "given the derivative of the target function at $x(0)$ and at the other end $x(n-1)$.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f = A(k)+B(k)*(x-x(k)) + C(k)*exp(p(k)*(x-x(k))) +D(k)*exp(-p(k)*(x-x(k)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Sample data to interpolate"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X =[0,xcp,1]\n",
      "Y =[0,ycp,1]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Set up each piece of interpolant"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "c=[f.subs(x(k),X[i]).subs(k,i) for i in range(2)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Left-side Interpolation conditions"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cond_i=[(Y[i]-f.subs(x,x(k)).subs(k,i)) for i in range(2)] # conditions for interpolation\n",
      "cond_i+= [ Y[2]- c[1].subs(x,X[2])]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Match end-point 1st derivatives from givens"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cond_end=[ diff(f,x).subs(k,0).subs(x(0),X[0]).subs(x,X[0]) - y1p,\n",
      "           diff(f,x).subs(k,1).subs(x(1),X[1]).subs(x,X[2]) - ynp,\n",
      " ]\n",
      "cond_end"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "[-y1p + B(0) + C(0)*p(0) - D(0)*p(0),\n",
        " -ynp + B(1) + C(1)*p(1)*exp((-xcp + 1)*p(1)) - D(1)*p(1)*exp(-(-xcp + 1)*p(1))]"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Inner continuity conditions"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cond_cont=[] # continuity conditions\n",
      "for i,j,v in zip(c[:-2],c[2:],range(1,3)):\n",
      "    cond_cont.append( (i-j).subs(x,v) )\n",
      "cond_cont\n",
      "\n",
      "cond_cont=[(c[0]-c[1]).subs(x,X[1])]\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Inner second derivatives must match"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "d2=[i.diff(x,2) for i in c]\n",
      "cond2_cont=[] # 2nd derivative continuity conditions\n",
      "for i,j,v in zip(d2[:-2],d2[2:],range(3)):\n",
      "    cond2_cont.append( (i-j).subs(x,v) )\n",
      "    \n",
      "cond2_cont= [diff(c[0]-c[1],x,2).subs(x,X[1])]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Inner first derivatives must match"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "d=[i.diff(x) for i in c]\n",
      "cond1_cont=[] # 2nd derivative continuity conditions\n",
      "for i,j,v in zip(d[:-2],d[2:],range(3)):\n",
      "    cond1_cont.append( (i-j).subs(x,v) )\n",
      "    \n",
      "cond1_cont=[diff(c[0]-c[1],x).subs(x,X[1])]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "len(cond_i)+len(cond_cont)+len(cond_end)+len(cond2_cont)+len(cond1_cont)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "8"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i in (cond_i+cond_cont+cond_end+cond2_cont):\n",
      "    print i.subs(p(k),0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-A(0) - C(0) - D(0)\n",
        "ycp - A(1) - C(1) - D(1)\n",
        "-(-xcp + 1)*B(1) - A(1) - C(1)*exp((-xcp + 1)*p(1)) - D(1)*exp(-(-xcp + 1)*p(1)) + 1\n",
        "xcp*B(0) + A(0) - A(1) + C(0)*exp(xcp*p(0)) - C(1) + D(0)*exp(-xcp*p(0)) - D(1)\n",
        "-y1p + B(0) + C(0)*p(0) - D(0)*p(0)\n",
        "-ynp + B(1) + C(1)*p(1)*exp((-xcp + 1)*p(1)) - D(1)*p(1)*exp(-(-xcp + 1)*p(1))\n",
        "C(0)*p(0)**2*exp(xcp*p(0)) - C(1)*p(1)**2 + D(0)*p(0)**2*exp(-xcp*p(0)) - D(1)*p(1)**2\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "linsys=(cond_i+cond_cont+cond_end+cond2_cont+cond1_cont)\n",
      "M=Matrix([[ l.coeff(i) for i in flatten([[A(j),B(j),C(j),D(j)] for j in range(2)]) ] for l in linsys])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sum([abs(diff(i,x,2)) for i in c]) # curvature metric"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "Abs(C(0)*p(0)**2*exp(x*p(0)) + D(0)*p(0)**2*exp(-x*p(0))) + Abs(C(1)*p(1)**2*exp((x - xcp)*p(1)) + D(1)*p(1)**2*exp(-(x - xcp)*p(1)))"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print c[0]\n",
      "print c[1]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "x*B(0) + A(0) + C(0)*exp(x*p(0)) + D(0)*exp(-x*p(0))\n",
        "(x - xcp)*B(1) + A(1) + C(1)*exp((x - xcp)*p(1)) + D(1)*exp(-(x - xcp)*p(1))\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    }
   ],
   "metadata": {}
  }
 ]
}