]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/interpolation.py
70c35c1a3ac9e812865185c71b1f4d2cc88be1ed
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 expression 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 expression (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(xs
, ys
):
98 Return the Newton divided difference of the points (xs[k],
101 http://en.wikipedia.org/wiki/Divided_differences
105 - ``xs`` -- The list of x-values.
107 - ``ys`` -- The function values at `xs`.
111 The (possibly symbolic) divided difference function.
116 sage: ys = map(sin, xs)
117 sage: divided_difference(xs, ys)
120 sage: ys = map(sin, xs)
121 sage: divided_difference(xs, ys)
123 sage: xs = [0, pi, 2*pi]
124 sage: ys = map(sin, xs)
125 sage: divided_difference(xs, ys)
128 We try something entirely symbolic::
130 sage: f = function('f', x)
131 sage: divided_difference([x], [f(x=x)])
133 sage: x1,x2 = var('x1,x2')
134 sage: divided_difference([x1,x2], [f(x=x1),f(x=x2)])
135 f(x1)/(x1 - x2) - f(x2)/(x1 - x2)
138 coeffs
= divided_difference_coefficients(xs
)
139 v_cs
= vector(coeffs
)
141 return v_cs
.dot_product(v_ys
)
144 def newton_polynomial(x
, xs
, ys
):
146 Return the Newton form of the interpolating polynomial of the
147 points (xs[k], ys[k]) in the variable `x`.
151 - ``x`` -- The independent variable to use for the interpolating
154 - ``xs`` -- The list of x-values.
156 - ``ys`` -- The function values at `xs`.
160 A symbolic expression.
164 sage: xs = [ -pi/2, -pi/6, 0, pi/6, pi/2 ]
165 sage: ys = map(sin, xs)
166 sage: L = lagrange_polynomial(x, xs, ys)
167 sage: N = newton_polynomial(x, xs, ys)
176 for k
in range(0, degree
+1):
177 term
= divided_difference(xs
[:k
+1], ys
[:k
+1])
178 term
*= product([ x
- xk
for xk
in xs
[:k
]])
184 def hermite_coefficient(k
, x
, xs
):
186 Return the Hermite coefficient h_{k}(x). See Atkinson, p. 160.
190 - ``k`` -- The index of the coefficient.
192 - ``x`` -- The symbolic variable to use as the argument of h_{k}.
194 - ``xs`` -- The list of points at which the function values are
199 A symbolic expression.
202 lk
= lagrange_coefficient(k
, x
, xs
)
203 return (1 - 2*lk
.diff(x
)(x
=xs
[k
])*(x
- xs
[k
]))*(lk
**2)
206 def hermite_deriv_coefficient(k
, x
, xs
):
208 Return the Hermite derivative coefficient, \tilde{h}_{k}(x). See
213 - ``k`` -- The index of the coefficient.
215 - ``x`` -- The symbolic variable to use as the argument of h_{k}.
217 - ``xs`` -- The list of points at which the function values are
222 A symbolic expression.
225 lk
= lagrange_coefficient(k
, x
, xs
)
226 return (x
- xs
[k
])*(lk
**2)
229 def hermite_interpolant(x
, xs
, ys
, y_primes
):
231 Return the Hermite interpolant `H(x)` such that H(xs[k]) = ys[k]
232 and H'(xs[k]) = y_primes[k] for each k.
234 Reference: Atkinson, p. 160.
238 - ``x`` -- The symbolic variable to use as the argument of H(x).
240 - ``xs`` -- The list of points at which the function values are
243 - ``ys`` -- The function values at the `xs`.
245 - ``y_primes`` -- The derivatives at the `xs`.
249 A symbolic expression.
253 sage: xs = [ 0, pi/6, pi/2 ]
254 sage: ys = map(sin, xs)
255 sage: y_primes = map(cos, xs)
256 sage: H = hermite_interpolant(x, xs, ys, y_primes)
257 sage: expected = -27/4*(pi - 6*x)*(pi - 2*x)^2*sqrt(3)*x^2/pi^4
258 sage: expected += (5*(pi - 2*x)/pi + 1)*(pi - 6*x)^2*x^2/pi^4
259 sage: expected += 81/2*((pi - 6*x)/pi + 1)*(pi - 2*x)^2*x^2/pi^4
260 sage: expected += (pi - 6*x)^2*(pi - 2*x)^2*x/pi^4
261 sage: bool(H == expected)
265 s1
= sum([ ys
[k
] * hermite_coefficient(k
, x
, xs
)
266 for k
in range(0, len(xs
)) ])
268 s2
= sum([ y_primes
[k
] * hermite_deriv_coefficient(k
, x
, xs
)
269 for k
in range(0, len(xs
)) ])