# [SciPy-User] Crossing of Splines

denis denis-bz-gg@t-online...
Thu Jun 10 09:16:47 CDT 2010

```On Jun 8, 8:22 pm, Marco <gae...@gmail.com> wrote:
> Hi all!
>
> I have 2 different datasets which I fit using interpolate.splrep().
>
> I am interested in finding the point where the splines cross: as of
> now I use interpolate.splev() to evaluate each spline and then look
> for a zero in the difference of the evaluated splines.

Following Anne's idea, if you extend class UnivariateSpline to use
given knots,
you could then subtract the whole splines / all their coefficients
and run roots() on the difference, along the lines

a = myUniSpline( x, ya, s=s, knots= )
b = myUniSpline( x, yb, s=s, knots= )
a.minus( b.getcoeffs() )  # spline a - b
roots = a.roots()

Marco, are you interpolating, s=0, or smoothing, s > 0 ?
If smoothing, knots wobble around with s and y (try the little test
below) --
good, adaptive, for interpolating 1 function, not so good for your
problem..

cheers
-- denis

""" scipy UnivariateSpline sensitivity to s """

from __future__ import division
import sys
import numpy as np
from scipy.interpolate import UnivariateSpline  # \$scipy/
interpolate/fitpack2.py
import pylab as pl

np.set_printoptions( 2, threshold=10, edgeitems=3, suppress=True )

N = 100
H = 5
cycle = 4
plot = 0
exec "\n".join( sys.argv[1:] )  # N= ...

x = np.arange(N+1)
xup = np.arange( 0, N + 1e-6, 1/H )
if cycle == 0:
y = np.zeros(N+1)
y[N//2] = 1
else:
y = np.sin( 2*np.pi * np.arange(N+1) / cycle )

#...............................................................................
title = "UnivariateSpline N=%d H=%d cycle=%.2g" % (N, H, cycle)
print title
for s in (0, .1, .5, 10, None ):
uspline = UnivariateSpline( x, y, s=s )  # s=0 interpolates
yup = uspline( xup )
res = uspline.get_residual()  # == s == |y - yup[::H]|**2
label = "s: %s  res: %.2g" % (s, res)
print label
knots = uspline.get_knots()
print "%d knots: %s" % (len(knots), knots )
roots = uspline.roots()
print "%d roots: %s" % (len(roots), roots )
print ""
if plot:
pl.plot( xup, yup, label=label )
if cycle == 0:
pl.xlim( N//2 - 10, N//2 + 10 )

if plot:
pl.title(title)
pl.legend()
pl.show()
```