]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/interpolation.py
1d594c38913595c213c7a0ffe2c37d42a658eae7
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(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`.
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
]])