1 {
3   "name": ""
4  },
5  "nbformat": 3,
6  "nbformat_minor": 0,
7  "worksheets": [
8   {
9    "cells": [
10     {
12      "level": 1,
14      "source": [
15       "Splines\n"
16      ]
17     },
18     {
19      "cell_type": "markdown",
21      "source": [
22       "Piecewise continuous curves. $S:[a,b]\\to \\Re$\n",
23       "\n",
24       "$S = P_i (t) \\quad t_{i-1} \\leq t \\lt t_i$ for $i=0,...k$\n",
25       "\n",
26       "$P_i$ are polynomials. Chosen to gaurantee smoothness of $S$. $P_i^{(j)} (t_i) = P_{i+1}^{(j)} (t_i)$\n"
27      ]
28     },
29     {
30      "cell_type": "markdown",
32      "source": [
33       "The spline is defined by $\\{q_0(x) ... q_n(x)\\}$ for $n$ points $x_i$.\n",
34       "\n",
35       "For cubic splines have: \n",
36       "\\begin{align}\n",
37       "    q_i(x) &= a_i x^3 + b_i x^2 + c_i x + d_i \\\\\n",
38       "    q_i(x_i) = y_i &= q_{i+1}(x_i) \\\\\n",
39       "    &\\implies x_i^3 (a_{i} - a_{i+1}) + x_i^2 (b_{i} - b_{i+1}) + x_i (c_{i} - c_{i+1}) + (d_{i} - d_{i+1}) = 0 \\quad \\forall i\\\\\n",
40       "    q'_i(x_i) &= q'_{i+1}(x_i) \\\\\n",
41       "     &\\implies 3 x_i^2 (a_{i} - a_{i+1}) + 2 x_i (b_{i} - b_{i+1}) + (c_{i} - c_{i+1}) = 0 \\quad \\forall i\\\\\n",
42       "     q''_i(x_i) &= q''_{i+1}(x_i) \\\\\n",
43       "     &\\implies 6 x_i (a_{i} - a_{i+1}) + 2 (b_{i} - b_{i+1}) = 0 \\quad \\forall i\n",
44       "\\end{align}\n"
45      ]
46     },
47     {
48      "cell_type": "markdown",
50      "source": [
51       "**Not useful:** Summing:\n",
52       "\n",
53       "\\begin{align}\n",
54       "    0 &= a_0 (x_0^3 + 3 x_0^2 + 6 x_0) + b_0 (x_0^2 + 2x_0+2) + c_0 (x_0 + 1) + d_0 \\\\\n",
55       "      &+ \\sum_{i=1}^{n} a_i \\left[ (x_i^3 - x_{i-1}^3) + 3(x_i^2 - x_{i-1}^2) + 6(x_i - x_{i-1})\\right]\n",
56       "      + b_i\\left[(x_i^2 - x_{i-1}^2) + 2(x_i - x_{i-1})\\right] + c_i\\left[x_i - x_{i-1}\\right] \\\\\n",
57       "\\end{align}"
58      ]
59     },
60     {
61      "cell_type": "markdown",
63      "source": [
64       "$$q_i(x) = a_i x^3 + b_i x^2 + c_i x + d_i$$"
65      ]
66     },
67     {
68      "cell_type": "markdown",
70      "source": [
71       "$$q_i'(x) = 3 a_i x^2 + 2 b_i x + c_i$$"
72      ]
73     },
74     {
75      "cell_type": "markdown",
77      "source": [
78       "$$q_i''(x) = 6 a_i x + 2 b_i$$"
79      ]
80     },
81     {
82      "cell_type": "markdown",
84      "source": [
85       "$$x_i^3 a_i + x_i^2 b_i + x_i c_i + d_i = y_i = x_i^3 a_{i+1} + x_i^2 b_{i+1} + x_i c_{i+1} + d_{i+1}$$"
86      ]
87     },
88     {
89      "cell_type": "code",
90      "collapsed": false,
91      "input": [],
92      "language": "python",
94      "outputs": []
95     },
96     {
97      "cell_type": "code",
98      "collapsed": false,
99      "input": [
100       "def Spline(P):\n",
101       "    \"\"\" Make a cubic spline that goes through the knots P \"\"\"\n",
102       "    x = [1.*p[0] for p in P]\n",
103       "    y = [1.*p[1] for p in P]\n",
104       "   \n",
105       "    n = len(P)-1\n",
106       "    \n",
107       "    A = zeros((n+1,)*2)\n",
108       "    B = 3. * ones(n+1)\n",
109       "    for i in xrange(1,n):\n",
110       "        A[i, i-1:i+2] = [1./(x[i] - x[i-1]), 2.*(1./(x[i]-x[i-1]) + 1./(x[i+1]-x[i])), 1./(x[i+1]-x[i])]\n",
111       "    \n",
112       "    A[0,0:2] = [2./(x[1]-x[0]), 1./(x[1]-x[0])]\n",
113       "    B[0] = 3. * ((y[1] - y[0]) / (x[1] - x[0])**2.)\n",
114       "    A[n,n-1:] = [1./(x[n] - x[n-1]), 2./(x[n]-x[n-1])]\n",
115       "    B[n] = 3. * ((y[n] - y[n-1]) / (x[n] - x[n-1])**2.)\n",
116       "    \n",
117       "    k = solve(A, B)[-1::-1]\n",
118       "\n",
119       "    \n",
120       "    def p(xx):\n",
121       "        for i in xrange(n):\n",
122       "            if xx >= x[i] and xx <= x[i+1]:\n",
123       "                break\n",
124       "        if (i >= n):\n",
125       "            i = n-1\n",
126       "            \n",
127       "        t = (xx - x[i]) / (x[i+1] - x[i])\n",
128       "        a = k[i]*(x[i+1] - x[i]) - (y[i+1] - y[i])\n",
129       "        b = -k[i+1]*(x[i+1]-x[i]) - (y[i+1] - y[i])\n",
130       "        return (1. - t)*y[i] + t*y[i+1] + t*(1.-t)*(a*(1.-t)+b*t)\n",
131       "        \n",
132       "    return p\n",
133       "                "
134      ],
135      "language": "python",
137      "outputs": [],
138      "prompt_number": 216
139     },
140     {
141      "cell_type": "code",
142      "collapsed": false,
143      "input": [
144       "points = [(-1,0),(0,-1),(1,0)]\n"
145      ],
146      "language": "python",
148      "outputs": [],
149      "prompt_number": 245
150     },
151     {
152      "cell_type": "code",
153      "collapsed": false,
154      "input": [
155       "p1 = Spline(points)"
156      ],
157      "language": "python",
159      "outputs": [],
160      "prompt_number": 246
161     },
162     {
163      "cell_type": "code",
164      "collapsed": false,
165      "input": [
166       "x = linspace(-1,1)\n",
167       "plot(x, map(p1, x))"
168      ],
169      "language": "python",
171      "outputs": [
172       {
174        "output_type": "pyout",
175        "prompt_number": 247,
176        "text": [
177         "[<matplotlib.lines.Line2D at 0x66545d0>]"
178        ]
179       },
180       {
182        "output_type": "display_data",
184        "text": [
185         "<matplotlib.figure.Figure at 0x63b66d0>"
186        ]
187       }
188      ],
189      "prompt_number": 247
190     },
191     {
192      "cell_type": "code",
193      "collapsed": false,
194      "input": [
195       "plot(x[:-1], diff(map(p1, x)))"
196      ],
197      "language": "python",
199      "outputs": [
200       {
202        "output_type": "pyout",
203        "prompt_number": 248,
204        "text": [
205         "[<matplotlib.lines.Line2D at 0x6837810>]"
206        ]
207       },
208       {
210        "output_type": "display_data",
212        "text": [
213         "<matplotlib.figure.Figure at 0x5bfe0d0>"
214        ]
215       }
216      ],
217      "prompt_number": 248
218     },
219     {
220      "cell_type": "code",
221      "collapsed": false,
222      "input": [
223       "plot(x[:-2], diff(diff(map(p1, x))))"
224      ],
225      "language": "python",
227      "outputs": [
228       {
230        "output_type": "pyout",
231        "prompt_number": 249,
232        "text": [
233         "[<matplotlib.lines.Line2D at 0x690d350>]"
234        ]
235       },
236       {
238        "output_type": "display_data",
240        "text": [
241         "<matplotlib.figure.Figure at 0x5886750>"
242        ]
243       }
244      ],
245      "prompt_number": 249
246     },
247     {
248      "cell_type": "code",
249      "collapsed": false,
250      "input": [],
251      "language": "python",
253      "outputs": [],
254      "prompt_number": 249
255     },
256     {
257      "cell_type": "code",
258      "collapsed": false,
259      "input": [],
260      "language": "python",