理论
#a,b,c,d are list of float that have same length
#b0 c0
#a1 b1 c1
# a2 b2 c2
# a3 b3
def Thomas(a, b, c, d):
length = len(a)
u = [b[0]]
l = [0.]
for i in range(1, length):
l.append(a[i]/u[i-1])
u.append(b[i]-l[i]*c[i-1])
y = [d[0]]
for i in range(1, length):
y.append(d[i] - l[i]*y[i-1])
x = [0.] * length
x[length-1] = y[length-1] / u[length-1]
for i in range(length-2, -1, -1):
x[i] = (y[i] - c[i]*x[i+1])/u[i]
return x
'''
x =
0.65909
1.63636
0.90909
'''
def test_Thomas():
a = [0., -4., -1.]
b = [4., 4., 4.]
c = [-1., -1., 0.]
d = [1., 3., 2.]
print(Thomas(a, b, c, d))
#x,y list of float, xq(x_query) float
def spline_natural(x, y, xq):
assert(x[0] <= xq and xq <= x[-1])
length = len(x)
d = [0.] * length
mu = [0.] * length
lamda = [0.] * length
h = [0.] * (length-1)
f = [0.] * (length-1)
for i in range(0, length-1):
h[i] = x[i+1] - x[i]
for i in range(0, length-1):
f[i] = (y[i+1] - y[i])/h[i]
for i in range(1, length-2):
lamda[i] = h[i]/(h[i-1]+h[i])
for i in range(2, length-1):
mu[i] = h[i-1]/(h[i-1]+h[i])
for i in range(1, length-1):
d[i] = 6 /(h[i-1]+h[i]) * (f[i] - f[i-1])
m = Thomas(mu[1:-1], [2]*(length-2), lamda[1:-1], d[1:-1])
m = [0] + m + [0]
idx = 0
for i in range(0, length-1):
if xq <= x[i+1]:
idx = i
break
i = idx
s = m[i] * (x[i+1] - xq)**3 /(6*h[i])
s += m[i+1] * (xq - x[i])**3 /(6*h[i])
s += (y[i] - m[i]*h[i]*h[i]/6) * (x[i+1] - xq)/h[i]
s += (y[i+1] - m[i+1]*h[i]*h[i]/6) * (xq -x[i])/h[i]
return s
def floatrange(start,stop,steps):
''' Computes a range of floating value.
Input:
start (float) : Start value.
end (float) : End value
steps (integer): Number of values
Output:
A list of floats
Example:
>>> print floatrange(0.25, 1.3, 5)
[0.25, 0.51249999999999996, 0.77500000000000002, 1.0375000000000001, 1.3]
'''
return [start+float(i)*(stop-start)/(float(steps)-1) for i in range(steps)]
def test_spline():
x = [1, 2, 4, 5]
y = [1, 3, 4, 2]
x_test = floatrange(1, 5, 100)
y_pre = []
for x_value in x_test:
y_pre.append(spline_natural(x,y,x_value))
import numpy as np
import matplotlib.pyplot as plt
plt.figure()
plt.plot(x,y,'r',x_test,y_pre,'b')
plt.show()
test_spline()