def computeC2Spline(nodes): '''Compute a cubic C2-spline by the list of 3 nodes. The spline consists of 2 cubic polynomials p0(x) and p1(x). p0(x) = p00 + p01*x + p02*x^2 + p03*x^3 p1(x) = p10 + p11*x + p12*x^2 + p13*x^3 The function returns the list of quadruples: [(p00, p01, p02, p03), (p10, p11, p12, p13)] The coefficients of polynomials are computed as the solution of the linear system. Let (x0, y0), (x1, y1), (x2, y2) are the nodes of interpolation. Then the system of linear equations is p0(x0) = y0 p0''(x0) = 0 p0(x1) = y1 p1(x1) = y1 p0'(x1) = p1'(x1) p0''(x1) = p1''(x1) p1(x2) = y2 p1''(x2) = 0 ''' assert(len(nodes) == 3) n = len(nodes) a = [] # The matrix is to be composed as a list of rows matrixLine = [0.]*(4*(n-1)) # The matrix row b = [0.]*(4*(n-1)) for i in range(len(matrixLine)): # Clear the matrix row matrixLine[i] = 0 x0 = nodes[0][0]; y0 = nodes[0][1] # Equation 0: # p0(x0) = y0 matrixLine[0] = 1 matrixLine[1] = x0 matrixLine[2] = x0^2 matrixLine[3] = x0^3 a.append(matrixLine.copy()) b[0] = y0 for i in range(len(matrixLine)): # Clear the matrix row matrixLine[i] = 0 # Equation 1: # p0''(x0) = 0 # 2*p02 + 6*p03*x0 = 0 matrixLine[2] = 2 matrixLine[3] = 6*x0 a.append(matrixLine.copy()) b[1] = 0 for i in range(len(matrixLine)): # Clear the matrix row matrixLine[i] = 0 x1 = nodes[1][0]; y1 = nodes[1][1] # Equation 2: # p0(x1) = y1 matrixLine[0] = 1 matrixLine[1] = x1 matrixLine[2] = x1^2 matrixLine[3] = x1^3 a.append(matrixLine.copy()) b[2] = y1 for i in range(len(matrixLine)): # Clear the matrix row matrixLine[i] = 0 # Equation 3: # p1(x1) = y1 matrixLine[4] = 1 matrixLine[5] = x1 matrixLine[6] = x1^2 matrixLine[7] = x1^3 a.append(matrixLine.copy()) b[3] = y1 for i in range(len(matrixLine)): # Clear the matrix row matrixLine[i] = 0 # Equation 4: # p0'(x1) - p1'(x1) = 0 matrixLine[1] = 1 matrixLine[2] = 2*x1 matrixLine[3] = 3*x1^2 matrixLine[5] = -1 matrixLine[6] = -2*x1 matrixLine[7] = -3*x1^2 a.append(matrixLine.copy()) b[4] = 0 for i in range(len(matrixLine)): # Clear the matrix row matrixLine[i] = 0 # Equation 5: # p0''(x1) - p1''(x1) = 0 matrixLine[2] = 2 matrixLine[3] = 6*x1 matrixLine[6] = -2 matrixLine[7] = -6*x1 a.append(matrixLine.copy()) b[5] = 0 for i in range(len(matrixLine)): # Clear the matrix row matrixLine[i] = 0 x2 = nodes[2][0]; y2 = nodes[2][1] # Equation 6: # p1(x2) = y2 matrixLine[4] = 1 matrixLine[5] = x2 matrixLine[6] = x2^2 matrixLine[7] = x2^3 a.append(matrixLine.copy()) b[6] = y2 for i in range(len(matrixLine)): # Clear the matrix row matrixLine[i] = 0 # Equation 7: # p1''(x2) = 0 matrixLine[6] = 2 matrixLine[7] = 6*x2 a.append(matrixLine.copy()) b[7] = 0 m = matrix(a) bb = vector(b) print("Matrix of the system:") print(m) print("Column of free terms:", bb) p = m.solve_right(vector(b)) print("Solution of the system:", p) s = [p[0:4], p[4:8]] print("Resulting spline:", s) return s def splineValue(nodes, spline, x): '''The value of the spline in the point x''' i = 0 n = len(nodes) # x_i <= x < x_i+1 while i < n-1 and x >= nodes[i+1][0]: i += 1 assert( (i == 0 and x < nodes[1]) or (i == n-1 and x >= nodes[n-1]) or nodes[i][0] <= x < nodes[i+1][0] ) p = spline[i] return p[0] + p[1]*x + p[2]*x^2 + p[3]*x^3 def splineC2(nodes): '''Plot the C2-spline constructed by nodes. This is the restricted version that works only in the case len(nodes) == 3''' assert(len(nodes) == 3) s = computeC2Spline(nodes) var("x") p0 = s[0][0] + s[0][1]*x + s[0][2]*x^2 + s[0][3]*x^3 shape = plot(p0, (x, nodes[0][0], nodes[1][0])) p1 = s[1][0] + s[1][1]*x + s[1][2]*x^2 + s[1][3]*x^3 shape += plot(p1, (x, nodes[1][0], nodes[2][0])) # Draw nodes for node in nodes: shape += point(node, size=40, color="red") return shape