]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/interpolation.py
2 load('~/.sage/init.sage')
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(f
, x
, xs
):
43 Return the Lagrange form of the interpolation polynomial in `x` of
44 `f` at the points `xs`.
48 - ``f`` - The function to interpolate.
50 - ``x`` - The independent variable of the resulting polynomial.
52 - ``xs`` - The list of points at which we interpolate `f`.
56 A symbolic function (polynomial) interpolating `f` at `xs`.
60 sage: xs = [ -pi/2, -pi/6, 0, pi/6, pi/2 ]
61 sage: L = lagrange_polynomial(sin, x, xs)
62 sage: expected = 27/16*(pi - 6*x)*(pi - 2*x)*(pi + 2*x)*x/pi^4
63 sage: expected -= 1/8*(pi - 6*x)*(pi - 2*x)*(pi + 6*x)*x/pi^4
64 sage: expected -= 1/8*(pi - 6*x)*(pi + 2*x)*(pi + 6*x)*x/pi^4
65 sage: expected += 27/16*(pi - 2*x)*(pi + 2*x)*(pi + 6*x)*x/pi^4
66 sage: bool(L == expected)
70 ys
= [ f(xs
[k
]) for k
in range(0, len(xs
)) ]
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
)) ])
76 def divided_difference(f
, xs
):
78 Return the Newton divided difference of `f` at the points
81 http://en.wikipedia.org/wiki/Divided_differences
85 - ``f`` -- The function whose divided difference we seek.
87 - ``xs`` -- The list of points at which to compute `f`.
91 The divided difference of `f` at ``xs``.
95 sage: divided_difference(sin, [0])
97 sage: divided_difference(sin, [0, pi])
99 sage: divided_difference(sin, [0, pi, 2*pi])
102 We try something entirely symbolic::
104 sage: f = function('f', x)
105 sage: divided_difference(f, [x])
107 sage: x1,x2 = var('x1,x2')
108 sage: divided_difference(f, [x1,x2])
109 (f(x1) - f(x2))/(x1 - x2)
113 # Avoid that goddamned DeprecationWarning when we have a named
114 # argument but don't know what it is.
115 if len(f
.variables()) == 0:
119 return f({ v: xs[0] }
)
121 # Use the recursive definition.
122 numerator
= divided_difference(f
, xs
[1:])
123 numerator
-= divided_difference(f
, xs
[:-1])
124 return numerator
/ (xs
[-1] - xs
[0])
127 def newton_polynomial(f
, x
, xs
):
129 Return the Newton form of the interpolating polynomial of `f` at
130 the points `xs` in the variable `x`.
134 - ``f`` -- The function to interpolate.
136 - ``x`` -- The independent variable to use for the interpolating
139 - ``xs`` -- The list of points at which to interpolate `f`.
147 sage: xs = [ -pi/2, -pi/6, 0, pi/6, pi/2 ]
148 sage: L = lagrange_polynomial(sin, x, xs)
149 sage: N = newton_polynomial(sin, x, xs)
158 for k
in range(0, degree
+1):
159 term
= divided_difference(f
, xs
[:k
+1])
160 term
*= product([ x
- xk
for xk
in xs
[:k
]])