#!/usr/bin/python # ---------------------------------------------------------------- # John Kerl # kerl.john.r@gmail.com # 2005-01-11 # # This is an example of Lagrange interpolation, 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) # # Program output: # # 0.00000 -0.34510 1.00000 -1.3451e+00 # 0.00391 -0.32857 1.00391 -1.3325e+00 # 0.00781 -0.31211 1.00784 -1.3200e+00 # ... # 0.98828 2.68467 2.68661 -1.9441e-03 # 0.99219 2.69585 2.69713 -1.2788e-03 # 0.99609 2.70705 2.70768 -6.3083e-04 # 1.00000 2.71828 2.71828 0.0000e+00 # 1.00391 2.72953 2.72892 6.1388e-04 # 1.00781 2.74081 2.73960 1.2110e-03 # 1.01172 2.75212 2.75032 1.7915e-03 # ... # 1.23828 3.46298 3.44968 1.3299e-02 # 1.24219 3.47640 3.46318 1.3215e-02 # 1.24609 3.48986 3.47674 1.3125e-02 # 1.25000 3.50337 3.49034 1.3029e-02 # 1.25391 3.51693 3.50400 1.2927e-02 # 1.25781 3.53054 3.51772 1.2819e-02 # 1.26172 3.54419 3.53149 1.2705e-02 # ... # 1.48828 4.43016 4.42948 6.8910e-04 # 1.49219 4.44727 4.44681 4.5811e-04 # 1.49609 4.46445 4.46422 2.2838e-04 # 1.50000 4.48169 4.48169 0.0000e+00 # 1.50391 4.49900 4.49923 -2.2694e-04 # 1.50781 4.51639 4.51684 -4.5236e-04 # 1.51172 4.53384 4.53452 -6.7618e-04 # ... # 1.98828 7.30223 7.30297 -7.4339e-04 # 1.99219 7.33106 7.33155 -4.9811e-04 # 1.99609 7.36000 7.36025 -2.5028e-04 # 2.00000 7.38906 7.38906 0.0000e+00 # 2.00391 7.41823 7.41798 2.5265e-04 # 2.00781 7.44752 7.44701 5.0758e-04 # 2.01562 7.50644 7.50542 1.0239e-03 # ... # 2.48828 12.04297 12.04056 2.4043e-03 # 2.49219 12.08932 12.08769 1.6277e-03 # 2.49609 12.13583 12.13500 8.2641e-04 # 2.50000 12.18249 12.18249 0.0000e+00 # 2.50391 12.22932 12.23017 -8.5189e-04 # 2.50781 12.27631 12.27804 -1.7296e-03 # 2.51172 12.32346 12.32610 -2.6336e-03 # ... # 2.98828 19.39450 19.85153 -4.5704e-01 # 2.99219 19.46416 19.92923 -4.6507e-01 # 2.99609 19.53404 20.00723 -4.7319e-01 # # Note that the fit is good at the interpolation points, but worse in between # them and even worse at the edges. # # When the data is plotted, examine the error curve and note that it has # appoximately the classic W shape of a quartic. This makes sense, since we've # used a cubic fit; the next-order term is quartic. # ================================================================ # 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 * import lagrinterp_m # ---------------------------------------------------------------- # Input mesh: xs = [0.20, 0.40, 0.60, 0.80] # Function to be approximated: f = atan # Interpolation mesh: xlo = 0.0 xhi = 1.0 dx = 1.0 / 100 # Coefficients for interpolating polynomial: cs = lagrinterp_m.li_make_coeffs_from_xs_and_f(xs, f) # Evaluate the function and the interpolant over the mesh: for x in mfrange(xlo, dx, xhi): y = f(x) z = lagrinterp_m.li_eval(x, xs, cs) print "%10.5f %10.5f %10.5f %10.3e" % (x, y, z, y-z)