# Import symbolics, plotting, and linear algebra functions.
from pylab import plot, savefig, grid, legend, clf, pcolor, spy
from numpy import arange, zeros
from sympy import (var, pi, factorial, sin, cos, log, exp, tan, atan, asin,
acos, sinh, cosh, tanh, asinh, acosh, lambdify)
from scipy.linalg import solve
# Use the symbol 'x' for independent variable.
var("x")
# Calculate the vector of spline coefficients based on
# a list of x-coordinates (x_list) and a list of the
# corresponding y-coordinates (y_list).
def calc_spline_coeffs(x_list, y_list, print_m, print_rhs, print_vec, spy_m):
nx = len(x_list)
ny = len(y_list)
if nx != ny:
print "Both lists must have equal length!"
print "Stopping."
return
n_pts = nx
n_elem = n_pts - 1
m = zeros((4*n_elem, 4*n_elem)) # zero matrix
rhs = zeros(4*n_elem) # zero vector
# fill the right-hand side vector
for i in range(0, n_pts-1): # loop over intervals
rhs[2*i] = y_list[i]
rhs[2*i+1] = y_list[i+1]
# fill the matrix
# Step 1: matching y-values at interval endpoints.
# This will generate the first 2*n_elem rows in the matrix.
for i in range(0, n_pts-1): # loop over intervals
xx = x_list[i]
m[2*i][4*i + 0] = 1
m[2*i][4*i + 1] = xx
m[2*i][4*i + 2] = xx**2
m[2*i][4*i + 3] = xx**3
xx = x_list[i+1]
m[2*i + 1][4*i + 0] = 1
m[2*i + 1][4*i + 1] = xx
m[2*i + 1][4*i + 2] = xx**2
m[2*i + 1][4*i + 3] = xx**3
# Step 2: matching first derivatives at all interior points.
# This will generate additional n_elem-1 rows in the matrix.
offset = 2*n_elem
for i in range(1, n_pts-1): # loop over internal points
xx = x_list[i]
m[offset + i-1][4*(i-1) + 1] = 1
m[offset + i-1][4*(i-1) + 2] = 2*xx
m[offset + i-1][4*(i-1) + 3] = 3*xx**2
m[offset + i-1][4*(i-1) + 5] = -1
m[offset + i-1][4*(i-1) + 6] = -2*xx
m[offset + i-1][4*(i-1) + 7] = -3*xx**2
# Step 3: matching second derivatives at all interior points
# This will generate additional n_elem-1 rows in the matrix.
offset = 2*n_elem + n_elem-1
for i in range(1, n_pts-1): # loop over internal points
xx = x_list[i]
m[offset + i-1][4*(i-1) + 2] = 2
m[offset + i-1][4*(i-1) + 3] = 6*xx
m[offset + i-1][4*(i-1) + 6] = -2
m[offset + i-1][4*(i-1) + 7] = -6*xx
# Step 4: Additional two conditions are needed to define
# a cubic spline. This will generate the last two rows in
# the matrix. Setting the second derivative (curvature) at both
# endpoints equal to zero will result into "natural cubic spline",
# but you can also prescribe non-zero values, or decide to
# prescribe first derivative (slope).
# Choose just one of the following two variables to be True,
# and state the corresponding value for the derivative.
set_second_derivative_left = True # if False, first derivative will be set
derivative_value_left = 0
# Choose just one of the following two variables to be True:
set_second_derivative_right = True # if False, first derivative will be set
derivative_value_right = 0
# Translate these two conditions into the last two matrix/rhs rows.
offset = 2*n_elem + 2*(n_elem-1)
xx = x_list[0] # left end-point
if set_second_derivative_left:
m[offset + 0][2] = 2
m[offset + 0][3] = 6*xx
rhs[-2] = derivative_value_left # value of the second derivative
else:
m[offset + 0][1] = 1
m[offset + 0][2] = 2*xx
m[offset + 0][3] = 3*xx**2
rhs[-2] = derivative_value_left # value of the first derivative
xx = x_list[-1] # right end-point
if set_second_derivative_right:
m[offset + 1][-2] = 2
m[offset + 1][-1] = 6*xx
rhs[-1] = derivative_value_right # value of the second derivative
else:
m[offset + 1][-3] = 1
m[offset + 1][-2] = 2*xx
m[offset + 1][-1] = 3*xx**2
rhs[-1] = derivative_value_right # value of the first derivative
# showing matrix m and right-hand side rhs
if print_m:
print "m ="
print m
if print_rhs:
print "rhs ="
print rhs
# solve the matrix system
vec = solve(m, rhs)
# showing solution vector vec
if print_vec:
print "solution vector ="
print vec
print "The cubic spline:"
for i in range(0, n_elem):
poly = vec[4*i] + vec[4*i+1]*x + vec[4*i+2]*x**2 + vec[4*i+3]*x**3
print " in (%g, %g): %s" % (x_list[i], x_list[i+1], poly)
if spy_m:
clf()
spy(m)
savefig("m.png")
return vec
# Elementary plotting function in interval (x_min, x_max).
# flag == 0: Plot C(x) only.
# flag == 1: Plot C(x) and f(x).
# flag == 2: Plot error E(x) = f(x) - C(x).
def plot_element(flag, a, b, c, d, f, x, x_min, x_max, show_label, m=100):
h = (x_max-x_min)/float(m)
x_array = arange(x_min, x_max+h, h)
if flag == 0:
y_array = [float(a + xx*(b + xx*(c + d*xx))) for xx in x_array]
if show_label:
label = "Cubic spline C(x)"
plot(x_array, y_array, "b-", label=label)
else: plot(x_array, y_array, "b-")
if flag == 1:
f_array = [f(xx) for xx in x_array]
if show_label:
label = "Function f(x)"
plot(x_array, f_array, "g-", label=label)
else:
plot(x_array, f_array, "g-")
if flag == 2:
e_array = [f(xx) - float(a + xx*(b + xx*(c + d*xx))) for xx in x_array]
if show_label:
label = "Error E(x) = f(x) - C(x)"
plot(x_array, e_array, "b-", label=label)
else:
plot(x_array, e_array, "b-")
e_abs = [abs(f(xx) - float(a + xx*(b + xx*(c + d*xx)))) for xx in x_array]
err_max = max(e_abs)
return err_max
# Plotting function
# flag == 0: Plot C(x)
# flag == 1: Plot f(x)
# flag == 2: Plot E(x) = f(x) - C(x)
def plot_spline(flag, x_list, vec, f, x, m=100):
n_pts = len(x_list)
n_elem = n_pts - 1
x_left = x_list[0]
x_right = x_list[-1]
h = float(x_right-x_left)/n_elem
err_max = 0
# loop over elements
for i in range(0, n_elem):
# i-th interval
x_min = x_list[i]
x_max = x_list[i+1]
# corresponding spline coefficients
a = vec[4*i]
b = vec[4*i+1]
c = vec[4*i+2]
d = vec[4*i+3]
# add label at the end
if i == n_elem-1: show_label = True
else: show_label = False
# do the element plot
if flag == 0: # plot C(x)
plot_element(0, a, b, c, d, None, None, x_min, x_max, show_label)
if flag == 1: # plot f(x)
plot_element(1, a, b, c, d, f, x, x_min, x_max, show_label)
if flag == 2: # plot E(x)
err_val = plot_element(2, a, b, c, d, f, x, x_min, x_max, show_label)
if err_val > err_max: err_max = err_val
return err_max
# General function: Construct and plot the cubic spline C(x). Also plot the error E(x) = f(x) - C(x) if
# y-coordinates are given via a function f(x)
# case == 0: Given is list of x-coordinates and list of y-coordinates.
# case == 1: Given is list of x-coordinates and a function f(x) to provide y-coordinates.
# case == 2: Given is interval (a,b), number of equally-spaced points, and a function f(x) to provide y-coordinates.
# case == 3: Given is interval (a,b), number of Chebyshev points, and a function f(x) to provide y-coordinates.
def cubic_spline(case, x_list, y_list, a, b, n_pts, f, x, print_m=True, print_rhs=True, print_vec=True, spy_m=True, m=100):
if case != 0 and case != 1 and case != 2 and case != 3:
print "Error: the first argument must be 0, 1, 2 or 3."
print "Stopping."
return
# construct x_list and/or y_list
if case == 1:
f = lambdify(x, f, modules=["math"])
y_list = [f(xx) for xx in x_list]
if case == 2:
n_elem = n_pts - 1
h = float(b-a)/n_elem
x_list = arange(a, b+h, h)
f = lambdify(x, f, modules=["math"])
y_list = [f(xx) for xx in x_list]
if case == 3:
n_elem = n_pts - 1
x_list = [float(cos(i*pi/n_elem)*(b-a)/2) for i in range(0, n_pts)]
f = lambdify(x, f, modules=["math"])
y_list = [f(xx) for xx in x_list]
# calculate the spline coefficients
vec = calc_spline_coeffs(x_list, y_list, print_m, print_rhs, print_vec, spy_m)
# plot the spline C(x), and/or function f(x) and error E(x)
if case == 0:
clf()
plot_spline(0, x_list, vec, None, None, m) # plot C(x)
legend()
grid(True)
savefig("a.png")
else:
clf()
plot_spline(0, x_list, vec, f, x, m) # plot C(x)
plot_spline(1, x_list, vec, f, x, m) # plot f(x)
legend()
grid(True)
savefig("a.png")
clf()
err_max = plot_spline(2, x_list, vec, f, x, m) # plot E(x)
legend()
grid(True)
savefig("err.png")
print "Max error:", err_max
print "Done."
return
# Given is list of x-coordinates and list of y-coordinates.
def cubic_spline_1(x_list, y_list, print_m=True, print_rhs=True, print_vec=True, spy_m=True, m=100):
cubic_spline(0, x_list, y_list, None, None, None, None, None, print_m, print_rhs, print_vec, spy_m, m)
return
# Given is list of x-coordinates and a function f(x) to provide y-coordinates.
def cubic_spline_2(x_list, f, x, print_m=True, print_rhs=True, print_vec=True, spy_m=True, m=100):
cubic_spline(1, x_list, None, None, None, None, f, x, print_m, print_rhs, print_vec, spy_m, m)
return
# Given is interval (a,b), number of equally-spaced points, and a function f(x) to provide y-coordinates.
def cubic_spline_3(a, b, n_pts, f, x, print_m=True, print_rhs=True, print_vec=True, spy_m=True, m=100):
cubic_spline(2, None, None, a, b, n_pts, f, x, print_m, print_rhs, print_vec, spy_m, m)
return
# Given is interval (a,b), number of Chebyshev points, and a function f(x) to provide y-coordinates.
def cubic_spline_4(a, b, n_pts, f, x, print_m=True, print_rhs=True, print_vec=True, spy_m=True, m=100):
cubic_spline(3, None, None, a, b, n_pts, f, x, print_m, print_rhs, print_vec, spy_m, m)
return
