%matplotlib inline
from pylab import *
from numpy import *
from numpy.random import normal, seed
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))
[<matplotlib.lines.Line2D at 0x7fd4140f4358>]
# the antiderivative matrix
def A(n):
antideriv = ones((n, n))
antideriv[triu_indices(n, k=1)] = 0.
return antideriv
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.]]
# check that the inverse of antiderivative approximates derivative
inv(A(10))
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.]])
AI = A(nelem)
AD = inv(AI)
df = inv(A(nelem)).dot(f)
dfn = inv(A(nelem)).dot(fn)
plot(x, df)
plot(x, dfn)
[<matplotlib.lines.Line2D at 0x7fd4141b90f0>]
#calculate the condition number
u, s, v = svd(inv(A(nelem)))
k = max(s)/min(s)
print(k)
51.5306511313
# ruzne mu 1. 5.
mu = 5.
xmu = inv(AI.transpose().dot(AI)+mu**2*eye(nelem)).dot(AI.transpose()).dot(fn)
plot(x, xmu)
plot(x, df)
plot(x, dfn)
[<matplotlib.lines.Line2D at 0x7fd40f471e48>]
mu = 4.
xmu = inv(AI.transpose().dot(AI)+mu**2*AD.transpose().dot(AD)).dot(AI.transpose()).dot(fn)
plot(x, xmu)
plot(x, df)
plot(x, dfn)
[<matplotlib.lines.Line2D at 0x7fd40f3dcbe0>]
from savitzky_golay import *
plot(x, xmu)
plot(x, df)
plot(x, dfn)
plot(x, savitzky_golay(f, 15, 3, 1, rate=1.))
[<matplotlib.lines.Line2D at 0x7fd40f354828>]
#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)
mu = 4.
xmu = inv(AI.transpose().dot(AI)+mu**2*AD.transpose().dot(AD)).dot(AI.transpose()).dot(fn)
plot(x, savitzky_golay(f, 15, 4, 1, rate=1.))
plot(x, xmu)
plot(x, df)
#plot(x, dfn)