In [1]:
%pylab inline
from sympy.utilities.lambdify import lambdify
from sympy import init_session
init_session()
from numpy import arange, abs
from pandas import DataFrame
from pylab import plot
import matplotlib.pyplot as plt
In [2]:
function = x*tanh(sin(x)+pi/3)
dfunction = function.diff(x)
d2function = dfunction.diff(x)
numFunction = lambdify(x, function, "numpy")
dnumFunction = lambdify(x, dfunction, "numpy")
d2numFunction = lambdify(x, d2function, "numpy")
d4f = lambdify(x, d2function.diff(x, 2), "numpy")
order = 200
t = pylab.linspace(-10, 10, 1000)
fig, axes = plt.subplots(figsize=(15,8))
axes.plot(t, numFunction(t))
Out[2]:
In [3]:
def maximum(f, start, end):
from scipy.optimize import minimize_scalar
return f(minimize_scalar(lambda x: -f(x), bounds=(start, end), method='bounded').x)
In [7]:
def calcM(points, values, lengths, diffOnEdges=False):
order = len(points)
W = ma.zeros((order-2, order-2))
for k in range(order-3):
W[k][k+1] = lengths[k+1]/6
W[k+1][k] = lengths[k+1]/6
W[k][k] = (lengths[k+1] + lengths[k])/3
W[order-3][order-3] = (lengths[order-2] + lengths[order-3])/3
B = ma.array([(values[m] - values[m+1])/lengths[m] + (values[m+2] - values[m+1])/lengths[m+1] for m in range(order-2)])
Q = dot(linalg.inv(W), B)
M = None
if diffOnEdges:
M = ma.concatenate([[d2numFunction(points[0])], ma.copy(Q), [d2numFunction(points[order - 1])]])
else:
M = ma.concatenate([[0], ma.copy(Q), [0]])
return M
def splinePoly(x, F, h, M, k, y):
if k < 0 or k > len(M)-2:
return 0
else:
return (M[k]*(x[k+1]-y)**3/(6*h[k])
+ M[k+1]*(y-x[k])**3/(6*h[k])
+ (F[k]-M[k]*h[k]**2/6)*(x[k+1]-y)/h[k]
+ (F[k+1]-M[k+1]*h[k]**2/6)*(y-x[k])/h[k])
def splinePolyDiff(x, F, h, M, k, y):
if k < 0 or k > len(M)-2:
return 0
else:
return (-M[k]*(x[k+1]-y)**2/(2*h[k])
+ M[k+1]*(y-x[k])**2/(2*h[k])
+ -(F[k]-M[k]*h[k]**2/6)/h[k]
+ (F[k+1]-M[k+1]*h[k]**2/6)/h[k])
def defineSpline(start, end, order, diffOnEdges=False):
points = pylab.linspace(start, end, order)
values = numFunction(points)
lengths = [points[i+1] - points[i] for i in range(0, len(points) - 1)]
M = calcM(points, values, lengths, diffOnEdges)
maxK = len(M) - 2
def spline(y):
k = int(floor((y - points[0])/lengths[0]))
return splinePoly(points, values, lengths, M, k, y)
def splineDiff(y):
k = int(floor((y - points[0])/lengths[0]))
return splinePolyDiff(points, values, lengths, M, k, y)
return pylab.vectorize(spline), pylab.vectorize(splineDiff)
In [9]:
spline = None
results = []
for order in range(9, 200, 10):
spline, splineDiff = defineSpline(-10, 10, order)
spline1, spline1Diff = defineSpline(-10, 10, order, True)
err = maximum(lambda x: abs(numFunction(x) - spline(x)), -10, 10)
err1 = maximum(lambda x: abs(numFunction(x) - spline1(x)), -10, 10)
derr = maximum(lambda x: abs(dnumFunction(x) - splineDiff(x)), -10, 10)
derr1 = maximum(lambda x: abs(dnumFunction(x) - spline1Diff(x)), -10, 10)
maxd4 = maximum(lambda x: abs(d4f(x)), -10, 10)
tv1 = ((10. - (-10.))/order)**4*maxd4
tv2 = ((10. - (-10.))/order)**3*maxd4
results.append((err, err1, derr, derr1, tv1, tv2))
fig, axes = plt.subplots(figsize=(15,8))
axes.plot(t, spline(t))
DataFrame(results, index=range(9,200,10), columns=("ΔS1", "ΔS2", "dΔS1", "dΔS2", "h^(4-0)*max|f^(4)|", "h^(4-1)*max|f^(4)|"))
Out[9]: