```python
from __future__ import print_function
from sympy import *
from eqn_manip import move_terms, move_terms_left, divide_terms, mult_eqn, get_coeff_for
from IPython.display import display
init_printing()
```

### Cubic spline equations
First step is make the a straightforward derivation of some cubic spline equations.  These could be further simplified to a tridiagonal matrix and efficiently solved.  For small numbers of knots and for testing purposes, it is sufficient to use a larger matrix and a general solver.

The first derivation will involve knots with constant separation of 1.  Later the derivation will be repeated with arbitrary knot positions.


```python
# Distance from the previous knot, for the case of uniform knot spacing
t = Symbol('t')

# Number of knots
n = Symbol('n', integer=True)
i = Symbol('i', integer=True)
# Function values to intepolated at the knots
y = IndexedBase('y')

# Coefficients of the spline function
a,b,c,d = [IndexedBase(s) for s in 'a b c d'.split()]

# Cubic spline function
s = a + b*t + c*t*t + d*t**3
display(Eq(y,s))

# With indexed variables
si = a[i] + b[i]*t + c[i]*t*t + d[i]*t**3
display(Eq(y[i],si))
```


$$y = t^{3} d + t^{2} c + t b + a$$



$${y}_{i} = t^{3} {d}_{i} + t^{2} {c}_{i} + t {b}_{i} + {a}_{i}$$


### Set up equations


```python
# Value at knots (t=0)
eq1 = Eq(si.subs(t,0), y[i])
eq1
```




$${a}_{i} = {y}_{i}$$




```python
# Value at knots (t=1)
eq2 = Eq(si.subs(t,1), y[i+1])
eq2
```




$${a}_{i} + {b}_{i} + {c}_{i} + {d}_{i} = {y}_{i + 1}$$




```python
# Continuity of first derivatives at knots 
dsi = diff(si, t)
eq3 = Eq(dsi.subs(t,1), dsi.subs(i, i+1).subs(t,0))
eq3
```




$${b}_{i} + 2 {c}_{i} + 3 {d}_{i} = {b}_{i + 1}$$




```python
# Continuity of second derivatives at knots
d2si = diff(si, t, 2)
eq4 = Eq(d2si.subs(t,1), d2si.subs(i, i+1).subs(t,0))
eq4
```




$$2 {c}_{i} + 6 {d}_{i} = 2 {c}_{i + 1}$$



### Boundary conditions


```python
# Natural BC (second derivative is zero at both ends)
eq5 = Eq(d2si.subs(i,0).subs(t,0), 0)
display(eq5)
eq6 = Eq(d2si.subs(i,n-1).subs(t,1), 0)
eq6
```


$$2 {c}_{0} = 0$$





$$2 {c}_{n - 1} + 6 {d}_{n - 1} = 0$$




```python
# Or
# Specified first derivatives at boundaries
eq5b = Eq(dsi.subs(i,0).subs(t,0),0)
display(eq5b)
eq6b = Eq(dsi.subs(i,n-1).subs(t,1),0)
eq6b
```


$${b}_{0} = 0$$





$${b}_{n - 1} + 2 {c}_{n - 1} + 3 {d}_{n - 1} = 0$$



### Set up equations

There are 4$n$ unknowns.  To get the right number of equations, we need 
 - one of eq5 and eq6 (boundary conditions) (2 equations)
 - eq1 and eq2 (values at beginning and end of each interval) for every knot ($2n$ equations)
 - eq3 and eq4 (continuity of 1st and 2nd derivatives) for every interior knot ($2(n-1)$ equations)

This gives the necessary $4n$ equations


```python
display(eq1)
display(eq2)
display(eq3)
display(eq4)
```


$${a}_{i} = {y}_{i}$$



$${a}_{i} + {b}_{i} + {c}_{i} + {d}_{i} = {y}_{i + 1}$$



$${b}_{i} + 2 {c}_{i} + 3 {d}_{i} = {b}_{i + 1}$$



$$2 {c}_{i} + 6 {d}_{i} = 2 {c}_{i + 1}$$



```python
# Equations 3 and 4 need to be rearranged so all the unknowns are the left-hand side
sym_lhs = [a[i],b[i],c[i],d[i],a[i+1],b[i+1],c[i+1],d[i+1]]
sym_rhs = [y[i]]
eq3b = divide_terms(eq3, sym_lhs, sym_rhs)
display(eq3b)
eq4b = divide_terms(eq4, sym_lhs, sym_rhs)
display(eq4b)
```


$$- {b}_{i + 1} + {b}_{i} + 2 {c}_{i} + 3 {d}_{i} = 0$$



$$- 2 {c}_{i + 1} + 2 {c}_{i} + 6 {d}_{i} = 0$$



```python
# Set up matrix with concrete number of knots
nknots = 3
nval = nknots -1
neqn = 4*nval
m = Matrix.eye(neqn, neqn)
rhs = Matrix.eye(neqn,1)
unknowns = Matrix.eye(neqn,1)
```


```python
use_natural_bc = True
symlist = list()
for idx in range(nval):
    symlist.append(a[idx])
    symlist.append(b[idx])
    symlist.append(c[idx])
    symlist.append(d[idx])
    
for idx,sym in enumerate(symlist):
    unknowns[idx] = sym

# First rows
for idx,sym in enumerate(symlist):
    if use_natural_bc:    
        m[0,idx] = get_coeff_for(eq5.lhs.subs(i,0), sym , symlist)
    else:
        m[0,idx] = get_coeff_for(eq5b.lhs.subs(i,0), sym, symlist)
    m[1,idx] = get_coeff_for(eq1.lhs.subs(i,0), sym,symlist)
    m[2,idx] = get_coeff_for(eq2.lhs.subs(i,0), sym,symlist)
    
rhs[0] = eq5.rhs.subs(i,0)
rhs[1] = eq1.rhs.subs(i,0)
rhs[2] = eq2.rhs.subs(i,0)

jdx = 3
for nv in range(1,nval):
    for idx,sym in enumerate(symlist):
        m[jdx+0,idx] = get_coeff_for(eq1.lhs.subs(i,nv), sym,symlist)
        m[jdx+1,idx] = get_coeff_for(eq2.lhs.subs(i,nv), sym,symlist)
        m[jdx+2,idx] = get_coeff_for(eq3b.lhs.subs(i,nv-1), sym,symlist)
        m[jdx+3,idx] = get_coeff_for(eq4b.lhs.subs(i,nv-1), sym,symlist)
    rhs[jdx+0] = eq1.rhs.subs(i,nv)
    rhs[jdx+1] = eq2.rhs.subs(i,nv)
    rhs[jdx+2] = eq3b.rhs.subs(i,nv)
    rhs[jdx+3] = eq4b.rhs.subs(i,nv)
    jdx += 4

# Last rows
for idx,sym in enumerate(symlist):
    if use_natural_bc:
        m[jdx,idx] = get_coeff_for(eq6.lhs.subs(n,nval), sym , symlist)
    else:
        m[jdx,idx] = get_coeff_for(eq6b.lhs.subs(n,nval), sym , symlist)
rhs[jdx] = eq6.rhs.subs(n,nval)

# Display of matrix-vector multiplication is not working - pretend the first equals is a multiplication
with evaluate(False):
    display(Eq(Eq(m,unknowns),rhs))
```


$$\left[\begin{matrix}0 & 0 & 2 & 0 & 0 & 0 & 0 & 0\\1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\1 & 1 & 1 & 1 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1 & 1 & 1 & 1\\0 & 1 & 2 & 3 & 0 & -1 & 0 & 0\\0 & 0 & 2 & 6 & 0 & 0 & -2 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 2 & 6\end{matrix}\right] = \left[\begin{matrix}{a}_{0}\\{b}_{0}\\{c}_{0}\\{d}_{0}\\{a}_{1}\\{b}_{1}\\{c}_{1}\\{d}_{1}\end{matrix}\right] = \left[\begin{matrix}0\\{y}_{0}\\{y}_{1}\\{y}_{1}\\{y}_{2}\\0\\0\\0\end{matrix}\right]$$



```python
# Solve
sln = m.LUsolve(rhs)
sln
```




$$\left[\begin{matrix}{y}_{0}\\- \frac{5 {y}_{0}}{4} + \frac{3 {y}_{1}}{2} - \frac{{y}_{2}}{4}\\0\\\frac{{y}_{0}}{4} - \frac{{y}_{1}}{2} + \frac{{y}_{2}}{4}\\{y}_{1}\\- \frac{{y}_{0}}{2} + \frac{{y}_{2}}{2}\\\frac{3 {y}_{0}}{4} - \frac{3 {y}_{1}}{2} + \frac{3 {y}_{2}}{4}\\- \frac{{y}_{0}}{4} + \frac{{y}_{1}}{2} - \frac{{y}_{2}}{4}\end{matrix}\right]$$




```python
# Can substitute concrete knot values
ysubs = {y[0]:0.5, y[1]:1.1, y[2]:1.5}
ysubs
```




$$\left \{ {y}_{0} : 0.5, \quad {y}_{1} : 1.1, \quad {y}_{2} : 1.5\right \}$$




```python
# Create substitutions for spline coefficients
subslist = dict(zip(symlist,sln))
subslist
```




$$\left \{ {a}_{0} : {y}_{0}, \quad {a}_{1} : {y}_{1}, \quad {b}_{0} : - \frac{5 {y}_{0}}{4} + \frac{3 {y}_{1}}{2} - \frac{{y}_{2}}{4}, \quad {b}_{1} : - \frac{{y}_{0}}{2} + \frac{{y}_{2}}{2}, \quad {c}_{0} : 0, \quad {c}_{1} : \frac{3 {y}_{0}}{4} - \frac{3 {y}_{1}}{2} + \frac{3 {y}_{2}}{4}, \quad {d}_{0} : \frac{{y}_{0}}{4} - \frac{{y}_{1}}{2} + \frac{{y}_{2}}{4}, \quad {d}_{1} : - \frac{{y}_{0}}{4} + \frac{{y}_{1}}{2} - \frac{{y}_{2}}{4}\right \}$$




```python
# For a particular interval
si0 = si.subs(i,0).subs(subslist)
display(si0)
si0b = si0.subs(ysubs)
display(si0b)
```


$$t^{3} \left(\frac{{y}_{0}}{4} - \frac{{y}_{1}}{2} + \frac{{y}_{2}}{4}\right) + t \left(- \frac{5 {y}_{0}}{4} + \frac{3 {y}_{1}}{2} - \frac{{y}_{2}}{4}\right) + {y}_{0}$$



$$- 0.05 t^{3} + 0.65 t + 0.5$$



```python
# Can take the seconnd derivative, to compare with the output values from some spline solvers
# At t=0, values is zero, as expected from natural boundary conditions
display(diff(si0b, t, 2).subs(t,0))
# Compute at t=1
display(diff(si0b, t, 2).subs(t,1))
```


$$0$$



$$-0.3$$


# Derive for general knot spacing


```python
# Distance from the previous knot
# t[i] ranges from 0 to x[i+1]-x[i]
t = IndexedBase('t')

# Number of knots
n = Symbol('n', integer=True)
i = Symbol('i', integer=True)
# Location of knots
x = IndexedBase('x')
# Function values to intepolated at the knots
y = IndexedBase('y')

# Coefficients of the spline function
a,b,c,d = [IndexedBase(s) for s in 'a b c d'.split()]


# With indexed variables
si = a[i] + b[i]*t[i] + c[i]*t[i]*t[i] + d[i]*t[i]**3
display(Eq(y[i],si))
```


$${y}_{i} = {a}_{i} + {b}_{i} {t}_{i} + {c}_{i} {t}_{i}^{2} + {d}_{i} {t}_{i}^{3}$$



```python
# Value at knots (t[i]=0)
eq1 = Eq(si.subs(t[i],0), y[i])
eq1
```




$${a}_{i} = {y}_{i}$$




```python
# Value at knots (t[i]=x[i+1])
eq2 = Eq(si.subs(t[i],x[i+1]-x[i]), y[i+1])
eq2
```




$$\left({x}_{i + 1} - {x}_{i}\right)^{3} {d}_{i} + \left({x}_{i + 1} - {x}_{i}\right)^{2} {c}_{i} + \left({x}_{i + 1} - {x}_{i}\right) {b}_{i} + {a}_{i} = {y}_{i + 1}$$




```python
# Continuity of first derivatives at knots 
dsi = diff(si, t[i])
eq3 = Eq(dsi.subs(t[i],x[i+1]-x[i]), dsi.subs(i, i+1).subs(t[i+1],0))
eq3
```




$$3 \left({x}_{i + 1} - {x}_{i}\right)^{2} {d}_{i} + 2 \left({x}_{i + 1} - {x}_{i}\right) {c}_{i} + {b}_{i} = {b}_{i + 1}$$




```python
# Continuity of second derivatives at knots
d2si = diff(si, t[i], 2)
eq4 = Eq(d2si.subs(t[i],x[i+1]-x[i]), d2si.subs(i, i+1).subs(t[i+1],0))
eq4
```




$$6 \left({x}_{i + 1} - {x}_{i}\right) {d}_{i} + 2 {c}_{i} = 2 {c}_{i + 1}$$




```python
# Natural BC (second derivative is zero at both ends)
eq5 = Eq(d2si.subs(i,0).subs(t[0],0), 0)
display(eq5)
eq6 = Eq(d2si.subs(t[i],x[i+1]-x[i]).subs(i,n-1), 0)
eq6
```


$$2 {c}_{0} = 0$$





$$6 \left(- {x}_{n - 1} + {x}_{n}\right) {d}_{n - 1} + 2 {c}_{n - 1} = 0$$




```python
# Or
# Specified first derivatives at boundaries
eq5b = Eq(dsi.subs(i,0).subs(t[0],0),0)
display(eq5b)
eq6b = Eq(dsi.subs(t[i],x[i+1]-x[i]).subs(i,n-1),0)
eq6b
```


$${b}_{0} = 0$$





$$3 \left(- {x}_{n - 1} + {x}_{n}\right)^{2} {d}_{n - 1} + 2 \left(- {x}_{n - 1} + {x}_{n}\right) {c}_{n - 1} + {b}_{n - 1} = 0$$




```python
# Equations 3 and 4 need to be rearranged so all the unknowns are the left-hand side
sym_lhs = [a[i],b[i],c[i],d[i],a[i+1],b[i+1],c[i+1],d[i+1]]
sym_rhs = [y[i]]
eq3b = divide_terms(eq3, sym_lhs, sym_rhs)
display(eq3b)
eq4b = divide_terms(eq4, sym_lhs, sym_rhs)
display(eq4b)
```


$$3 \left({x}_{i + 1} - {x}_{i}\right)^{2} {d}_{i} + 2 \left({x}_{i + 1} - {x}_{i}\right) {c}_{i} - {b}_{i + 1} + {b}_{i} = 0$$



$$6 \left({x}_{i + 1} - {x}_{i}\right) {d}_{i} - 2 {c}_{i + 1} + 2 {c}_{i} = 0$$



```python
display(eq1)
display(eq2)
display(eq3b)
display(eq4b)
```


$${a}_{i} = {y}_{i}$$



$$\left({x}_{i + 1} - {x}_{i}\right)^{3} {d}_{i} + \left({x}_{i + 1} - {x}_{i}\right)^{2} {c}_{i} + \left({x}_{i + 1} - {x}_{i}\right) {b}_{i} + {a}_{i} = {y}_{i + 1}$$



$$3 \left({x}_{i + 1} - {x}_{i}\right)^{2} {d}_{i} + 2 \left({x}_{i + 1} - {x}_{i}\right) {c}_{i} - {b}_{i + 1} + {b}_{i} = 0$$



$$6 \left({x}_{i + 1} - {x}_{i}\right) {d}_{i} - 2 {c}_{i + 1} + 2 {c}_{i} = 0$$



```python
# Set up matrix with concrete number of knots
nval = 2
nknots = nval + 1
neqn = 4*nval
m = Matrix.eye(neqn, neqn)
rhs = Matrix.eye(neqn,1)
unknowns = Matrix.eye(neqn,1)
```


```python
use_natural_bc = True
symlist = list()
for idx in range(nval):
    symlist.append(a[idx])
    symlist.append(b[idx])
    symlist.append(c[idx])
    symlist.append(d[idx])
    
for idx,sym in enumerate(symlist):
    unknowns[idx] = sym

# First rows
for idx,sym in enumerate(symlist):
    if use_natural_bc:    
        m[0,idx] = get_coeff_for(eq5.lhs.subs(i,0), sym , symlist)
    else:
        m[0,idx] = get_coeff_for(eq5b.lhs.subs(i,0), sym, symlist)
    m[1,idx] = get_coeff_for(eq1.lhs.subs(i,0), sym,symlist)
    m[2,idx] = get_coeff_for(eq2.lhs.subs(i,0), sym,symlist)
    
rhs[0] = eq5.rhs.subs(i,0)
rhs[1] = eq1.rhs.subs(i,0)
rhs[2] = eq2.rhs.subs(i,0)

jdx = 3
for nv in range(1,nval):
    for idx,sym in enumerate(symlist):
        m[jdx+0,idx] = get_coeff_for(eq1.lhs.subs(i,nv), sym,symlist)
        m[jdx+1,idx] = get_coeff_for(eq2.lhs.subs(i,nv), sym,symlist)
        m[jdx+2,idx] = get_coeff_for(eq3b.lhs.subs(i,nv-1), sym,symlist)
        m[jdx+3,idx] = get_coeff_for(eq4b.lhs.subs(i,nv-1), sym,symlist)
    rhs[jdx+0] = eq1.rhs.subs(i,nv)
    rhs[jdx+1] = eq2.rhs.subs(i,nv)
    rhs[jdx+2] = eq3b.rhs.subs(i,nv)
    rhs[jdx+3] = eq4b.rhs.subs(i,nv)
    jdx += 4

# Last rows
for idx,sym in enumerate(symlist):
    if use_natural_bc:
        m[jdx,idx] = get_coeff_for(eq6.lhs.subs(n,nval), sym , symlist)
    else:
        m[jdx,idx] = get_coeff_for(eq6b.lhs.subs(n,nval), sym , symlist)
rhs[jdx] = eq6.rhs.subs(n,nval)

# Display just the matrix, displaying other pieces doesn't work, and they are the same as previously
m
```




$$\left[\begin{matrix}0 & 0 & 2 & 0 & 0 & 0 & 0 & 0\\1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\1 & - {x}_{0} + {x}_{1} & \left(- {x}_{0} + {x}_{1}\right)^{2} & \left(- {x}_{0} + {x}_{1}\right)^{3} & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1 & - {x}_{1} + {x}_{2} & \left(- {x}_{1} + {x}_{2}\right)^{2} & \left(- {x}_{1} + {x}_{2}\right)^{3}\\0 & 1 & - 2 {x}_{0} + 2 {x}_{1} & 3 \left(- {x}_{0} + {x}_{1}\right)^{2} & 0 & -1 & 0 & 0\\0 & 0 & 2 & - 6 {x}_{0} + 6 {x}_{1} & 0 & 0 & -2 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 2 & - 6 {x}_{1} + 6 {x}_{2}\end{matrix}\right]$$




```python
# Solve
sln = m.LUsolve(rhs)
# Display is large - skip for now, uncomment if you want to see it
#sln
```


```python
# Create substitutions for spline coefficients
subslist = dict(zip(symlist,sln))
#subslist
```


```python
# Can substitute concrete knot values
ysubs = {y[0]:0.5, y[1]:1.1, y[2]:1.5}
display(ysubs)
# Use uniform knots to verify answer is the same as before
xsubs = {x[0]:0.0, x[1]:1.0, x[2]:2.0}
display(xsubs)
```


$$\left \{ {y}_{0} : 0.5, \quad {y}_{1} : 1.1, \quad {y}_{2} : 1.5\right \}$$



$$\left \{ {x}_{0} : 0.0, \quad {x}_{1} : 1.0, \quad {x}_{2} : 2.0\right \}$$



```python
# For a particular interval
si0 = si.subs(i,0).subs(subslist)
#display(si0)
si0b = si0.subs(ysubs).subs(xsubs).subs(t[0],t)
display(si0b)
```


$$- 0.05 t^{3} + 0.65 t + 0.5$$



```python

```


```python

```

a7ffb1a90d07ca75389aaa76a5669f453ad6ffc5
ipynb
Jupyter Notebook
Wavefunctions/Cubic Splines Basic.ipynb
QMCPACK/qmc_algorithms
015fd1973e94f98662149418adc6b06dcd78946d
MIT
2018-02-06T06:15:19.000Z
2019-11-26T23:54:53.000Z
Wavefunctions/Cubic Splines Basic.ipynb
chrinide/qmc_algorithms
015fd1973e94f98662149418adc6b06dcd78946d
MIT
Wavefunctions/Cubic Splines Basic.ipynb
chrinide/qmc_algorithms
015fd1973e94f98662149418adc6b06dcd78946d
MIT
2017-11-14T20:25:00.000Z
2022-02-28T06:02:01.000Z
Qwen/Qwen-72B
1. YES
2. YES