]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/interpolation.py
2 from misc
import product
4 def lagrange_coefficient(k
, x
, xs
):
6 Returns the coefficient function l_{k}(variable) of y_{k} in the
7 Lagrange polynomial of f. See,
9 http://en.wikipedia.org/wiki/Lagrange_polynomial
15 - ``k`` -- the index of the coefficient.
17 - ``x`` -- the symbolic variable to use for the first argument
20 - ``xs`` -- The list of points at which the function values are
25 A symbolic function of one variable.
29 sage: xs = [ -pi/2, -pi/6, 0, pi/6, pi/2 ]
30 sage: lagrange_coefficient(0, x, xs)
31 1/8*(pi - 6*x)*(pi - 2*x)*(pi + 6*x)*x/pi^4
34 numerator
= product([x
- xs
[j
] for j
in range(0, len(xs
)) if j
!= k
])
35 denominator
= product([xs
[k
] - xs
[j
] for j
in range(0, len(xs
)) if j
!= k
])
37 return (numerator
/ denominator
)
41 def lagrange_polynomial(x
, xs
, ys
):
43 Return the Lagrange form of the interpolation polynomial in `x` of
44 at the points (xs[k], ys[k]).
48 - ``x`` - The independent variable of the resulting polynomial.
50 - ``xs`` - The list of points at which we interpolate `f`.
52 - ``ys`` - The function values at `xs`.
56 A symbolic function (polynomial) interpolating each (xs[k], ys[k]).
60 sage: xs = [ -pi/2, -pi/6, 0, pi/6, pi/2 ]
61 sage: ys = map(sin, xs)
62 sage: L = lagrange_polynomial(x, xs, ys)
63 sage: expected = 27/16*(pi - 6*x)*(pi - 2*x)*(pi + 2*x)*x/pi^4
64 sage: expected -= 1/8*(pi - 6*x)*(pi - 2*x)*(pi + 6*x)*x/pi^4
65 sage: expected -= 1/8*(pi - 6*x)*(pi + 2*x)*(pi + 6*x)*x/pi^4
66 sage: expected += 27/16*(pi - 2*x)*(pi + 2*x)*(pi + 6*x)*x/pi^4
67 sage: bool(L == expected)
71 ls
= [ lagrange_coefficient(k
, x
, xs
) for k
in range(0, len(xs
)) ]
72 sigma
= sum([ ys
[k
] * ls
[k
] for k
in range(0, len(xs
)) ])
77 def divided_difference_coefficients(xs
):
79 Assuming some function `f`, compute the coefficients of the
80 divided difference f[xs[0], ..., xs[n]].
84 sage: divided_difference_coefficients([0])
86 sage: divided_difference_coefficients([0, pi])
88 sage: divided_difference_coefficients([0, pi, 2*pi])
89 [1/2/pi^2, -1/pi^2, 1/2/pi^2]
92 coeffs
= [ product([ (QQ(1) / (xj
- xi
)) for xi
in xs
if xi
!= xj
])
96 def divided_difference(f
, xs
):
98 Return the Newton divided difference of `f` at the points
101 http://en.wikipedia.org/wiki/Divided_differences
105 - ``f`` -- The function whose divided difference we seek.
107 - ``xs`` -- The list of points at which to compute `f`.
111 The divided difference of `f` at ``xs``.
115 sage: divided_difference(sin, [0])
117 sage: divided_difference(sin, [0, pi])
119 sage: divided_difference(sin, [0, pi, 2*pi])
122 We try something entirely symbolic::
124 sage: f = function('f', x)
125 sage: divided_difference(f, [x])
127 sage: x1,x2 = var('x1,x2')
128 sage: divided_difference(f, [x1,x2])
129 (f(x1) - f(x2))/(x1 - x2)
133 # Avoid that goddamned DeprecationWarning when we have a named
134 # argument but don't know what it is.
135 if len(f
.variables()) == 0:
139 return f({ v: xs[0] }
)
141 # Use the recursive definition.
142 numerator
= divided_difference(f
, xs
[1:])
143 numerator
-= divided_difference(f
, xs
[:-1])
144 return numerator
/ (xs
[-1] - xs
[0])
147 def newton_polynomial(f
, x
, xs
):
149 Return the Newton form of the interpolating polynomial of `f` at
150 the points `xs` in the variable `x`.
154 - ``f`` -- The function to interpolate.
156 - ``x`` -- The independent variable to use for the interpolating
159 - ``xs`` -- The list of points at which to interpolate `f`.
167 sage: xs = [ -pi/2, -pi/6, 0, pi/6, pi/2 ]
168 sage: ys = map(sin, xs)
169 sage: L = lagrange_polynomial(x, xs, ys)
170 sage: N = newton_polynomial(sin, x, xs)
179 for k
in range(0, degree
+1):
180 term
= divided_difference(f
, xs
[:k
+1])
181 term
*= product([ x
- xk
for xk
in xs
[:k
]])