#!/usr/bin/python # ================================================================ # John Kerl # kerl.john.r@gmail.com # 2005-01-11 # # This is an Lagrange-interpolation module. # # Example using degree 3: # Given x_0 .. x_3, y_0 .. y_3, compute the polynomial # # (x - x_1)(x - x_2)(x - x_3) # y_0 ------------------------------------ + # (x_0-x_1)(x_0-x_2)(x_0-x_3) # # (x - x_0) (x - x_2)(x - x_3) # y_1 ------------------------------------ + # (x_1-x_0) (x_1-x_2)(x_1-x_3) # # (x - x_0)(x - x_1) (x - x_3) # y_2 ------------------------------------ + # (x_2-x_0)(x_2-x_1) (x_2-x_3) # # (x - x_0)(x - x_1)(x - x_2) # y_3 ------------------------------------ # (x_3-x_0)(x_3-x_1)(x_3-x_2) # ---------------------------------------------------------------- # This software is released under the terms of the GNU GPL. # Please see LICENSE.txt in the same directory as this file. # ================================================================ from __future__ import division # 1/2 = 0.5, not 0. from math import * from kerlutil import * # ---------------------------------------------------------------- # Cubic case: # c0 = y0 / ( (x0-x1)*(x0-x2)*(x0-x3)) # c1 = y1 / ((x1-x0) *(x1-x2)*(x1-x3)) # c2 = y2 / ((x2-x0)*(x2-x1) *(x2-x3)) # c3 = y3 / ((x3-x0)*(x3-x1)*(x3-x2)) def li_make_coeffs_from_xs_and_ys(xs, ys): n = len(xs) cs = [0] * n for i in range(0, n): denom = 1.0 for j in range(0, n): if (i != j): denom *= xs[i] - xs[j] cs[i] = ys[i] / denom return cs # ---------------------------------------------------------------- def li_make_coeffs_from_xs_and_f(xs, f): return li_make_coeffs_from_xs_and_ys(xs, map(f, xs)) # ---------------------------------------------------------------- # Cubic case: # c0 * (x-x1)*(x-x2)*(x-x3) + \ # c1 * (x-x0) *(x-x2)*(x-x3) + \ # c2 * (x-x0)*(x-x1) *(x-x3) + \ # c3 * (x-x0)*(x-x1)*(x-x2) def li_eval(x, xs, cs): n = len(xs) sum = 0.0 for i in range(0, n): term = cs[i] for j in range(0, n): if (i != j): term *= x - xs[j] sum += term return sum