In [1]:

    

%matplotlib inline
from pylab import *
from numpy import *
from numpy.random import normal, seed


    

In [2]:

    

nelem = 40
noise = 0.2

x = linspace(0, 10, nelem)
f = sin(x)
seed(1234)
fn = f + normal(0.0, noise, size=nelem)

subplot(121)
plot(x, fn)

subplot(122)
plot(x[:-1], diff(f))
plot(x[:-1], diff(fn))


    

    
Out[2]:





[<matplotlib.lines.Line2D at 0x7fd4140f4358>]






    

In [3]:

    

# the antiderivative matrix
def A(n):
    antideriv = ones((n, n))
    antideriv[triu_indices(n, k=1)] = 0.
    return antideriv


    

In [4]:

    

print(A(10))


    

    




[[ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 1.  1.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 1.  1.  1.  0.  0.  0.  0.  0.  0.  0.]
 [ 1.  1.  1.  1.  0.  0.  0.  0.  0.  0.]
 [ 1.  1.  1.  1.  1.  0.  0.  0.  0.  0.]
 [ 1.  1.  1.  1.  1.  1.  0.  0.  0.  0.]
 [ 1.  1.  1.  1.  1.  1.  1.  0.  0.  0.]
 [ 1.  1.  1.  1.  1.  1.  1.  1.  0.  0.]
 [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  0.]
 [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]]

In [5]:

    

# check that the inverse of antiderivative approximates derivative
inv(A(10))


    

    
Out[5]:





array([[ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [-1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.]])

In [7]:

    

AI = A(nelem)
AD = inv(AI)


    

In [8]:

    

df = inv(A(nelem)).dot(f)
dfn = inv(A(nelem)).dot(fn)


    

In [9]:

    

plot(x, df)
plot(x, dfn)


    

    
Out[9]:





[<matplotlib.lines.Line2D at 0x7fd4141b90f0>]






    

In [10]:

    

#calculate the condition number
u, s, v = svd(inv(A(nelem)))
k = max(s)/min(s)
print(k)


    

    




51.5306511313

In [15]:

    

# ruzne mu 1. 5.
mu = 5.
xmu = inv(AI.transpose().dot(AI)+mu**2*eye(nelem)).dot(AI.transpose()).dot(fn)


    

In [16]:

    

plot(x, xmu)
plot(x, df)
plot(x, dfn)


    

    
Out[16]:





[<matplotlib.lines.Line2D at 0x7fd40f471e48>]






    

In [17]:

    

mu = 4.
xmu = inv(AI.transpose().dot(AI)+mu**2*AD.transpose().dot(AD)).dot(AI.transpose()).dot(fn)


    

In [18]:

    

plot(x, xmu)
plot(x, df)
plot(x, dfn)


    

    
Out[18]:





[<matplotlib.lines.Line2D at 0x7fd40f3dcbe0>]






    

In [19]:

    

from savitzky_golay import *
plot(x, xmu)
plot(x, df)
plot(x, dfn)
plot(x, savitzky_golay(f, 15, 3, 1, rate=1.))


    

    
Out[19]:





[<matplotlib.lines.Line2D at 0x7fd40f354828>]






    

In [ ]:

    

#another example:
nelem = 40
noise = 0.2

xmin = -1
xmax = 1
x = linspace(xmin, xmax, nelem)
f = abs(x) - abs(xmin)
fn = f + normal(0.0, noise, size=nelem)

AI = A(nelem)
AD = inv(AI)

df = AD.dot(f)
dfn = AD.dot(fn)


    

In [ ]:

    

mu = 4.
xmu = inv(AI.transpose().dot(AI)+mu**2*AD.transpose().dot(AD)).dot(AI.transpose()).dot(fn)


    

In [ ]:

    

plot(x, savitzky_golay(f, 15, 4, 1, rate=1.))
plot(x, xmu)
plot(x, df)
#plot(x, dfn)


    

In [ ]: