X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Finterpolation.py;h=6e0ecb3d79c5129a2e6270177b648142085786d1;hb=HEAD;hp=4be95e2463b7d4b1592a5f4e6cf86ba54e36da7f;hpb=37a881b845ebe285b70969db3636a250e64138ca;p=sage.d.git diff --git a/mjo/interpolation.py b/mjo/interpolation.py index 4be95e2..6e0ecb3 100644 --- a/mjo/interpolation.py +++ b/mjo/interpolation.py @@ -18,7 +18,7 @@ def lagrange_denominator(k, xs): The product of all xs[j] with j != k. """ - return product([xs[k] - xs[j] for j in range(0, len(xs)) if j != k]) + return product( xs[k] - xs[j] for j in range(len(xs)) if j != k ) def lagrange_coefficient(k, x, xs): @@ -44,6 +44,10 @@ def lagrange_coefficient(k, x, xs): A symbolic expression of one variable. + SETUP:: + + sage: from mjo.interpolation import lagrange_coefficient + TESTS:: sage: xs = [ -pi/2, -pi/6, 0, pi/6, pi/2 ] @@ -75,10 +79,14 @@ def lagrange_polynomial(x, xs, ys): A symbolic expression (polynomial) interpolating each (xs[k], ys[k]). + SETUP:: + + sage: from mjo.interpolation import lagrange_polynomial + TESTS:: sage: xs = [ -pi/2, -pi/6, 0, pi/6, pi/2 ] - sage: ys = map(sin, xs) + sage: ys = list(map(sin, xs)) sage: L = lagrange_polynomial(x, xs, ys) sage: expected = 27/16*(pi - 6*x)*(pi - 2*x)*(pi + 2*x)*x/pi^4 sage: expected -= 1/8*(pi - 6*x)*(pi - 2*x)*(pi + 6*x)*x/pi^4 @@ -88,9 +96,8 @@ def lagrange_polynomial(x, xs, ys): True """ - ls = [ lagrange_coefficient(k, x, xs) for k in range(0, len(xs)) ] - sigma = sum([ ys[k] * ls[k] for k in range(0, len(xs)) ]) - return sigma + ls = [ lagrange_coefficient(k, x, xs) for k in range(len(xs)) ] + return sum( ys[k] * ls[k] for k in range(len(xs)) ) @@ -111,6 +118,10 @@ def lagrange_interpolate(f, x, xs): A polynomial in ``x`` which interpolates ``f`` at ``xs``. + SETUP:: + + sage: from mjo.interpolation import lagrange_interpolate + EXAMPLES: We're exact on polynomials of degree `n` if we use `n+1` points:: @@ -135,7 +146,11 @@ def divided_difference_coefficients(xs): Assuming some function `f`, compute the coefficients of the divided difference f[xs[0], ..., xs[n]]. - TESTS: + SETUP:: + + sage: from mjo.interpolation import divided_difference_coefficients + + TESTS:: sage: divided_difference_coefficients([0]) [1] @@ -145,8 +160,7 @@ def divided_difference_coefficients(xs): [1/2/pi^2, -1/pi^2, 1/2/pi^2] """ - coeffs = [ QQ(1)/lagrange_denominator(k, xs) for k in range(0, len(xs)) ] - return coeffs + return [ ~lagrange_denominator(k, xs) for k in range(len(xs)) ] def divided_difference(xs, ys): @@ -166,18 +180,22 @@ def divided_difference(xs, ys): The (possibly symbolic) divided difference function. + SETUP:: + + sage: from mjo.interpolation import divided_difference + TESTS:: sage: xs = [0] - sage: ys = map(sin, xs) + sage: ys = list(map(sin, xs)) sage: divided_difference(xs, ys) 0 sage: xs = [0, pi] - sage: ys = map(sin, xs) + sage: ys = list(map(sin, xs)) sage: divided_difference(xs, ys) 0 sage: xs = [0, pi, 2*pi] - sage: ys = map(sin, xs) + sage: ys = list(map(sin, xs)) sage: divided_difference(xs, ys) 0 @@ -215,26 +233,22 @@ def newton_polynomial(x, xs, ys): A symbolic expression. - TESTS: + SETUP:: + + sage: from mjo.interpolation import lagrange_polynomial, newton_polynomial + + TESTS:: sage: xs = [ -pi/2, -pi/6, 0, pi/6, pi/2 ] - sage: ys = map(sin, xs) + sage: ys = list(map(sin, xs)) sage: L = lagrange_polynomial(x, xs, ys) sage: N = newton_polynomial(x, xs, ys) sage: bool(N == L) True """ - degree = len(xs) - 1 - - N = SR(0) - - for k in range(0, degree+1): - term = divided_difference(xs[:k+1], ys[:k+1]) - term *= lagrange_psi(x, xs[:k]) - N += term - - return N + return sum( divided_difference(xs[:k+1], ys[:k+1])*lagrange_psi(x, xs[:k]) + for k in range(len(xs)) ) def hermite_coefficient(k, x, xs): @@ -304,11 +318,15 @@ def hermite_interpolant(x, xs, ys, y_primes): A symbolic expression. - TESTS: + SETUP:: + + sage: from mjo.interpolation import hermite_interpolant + + TESTS:: sage: xs = [ 0, pi/6, pi/2 ] - sage: ys = map(sin, xs) - sage: y_primes = map(cos, xs) + sage: ys = list(map(sin, xs)) + sage: y_primes = list(map(cos, xs)) sage: H = hermite_interpolant(x, xs, ys, y_primes) sage: expected = -27/4*(pi - 6*x)*(pi - 2*x)^2*sqrt(3)*x^2/pi^4 sage: expected += (5*(pi - 2*x)/pi + 1)*(pi - 6*x)^2*x^2/pi^4 @@ -318,11 +336,11 @@ def hermite_interpolant(x, xs, ys, y_primes): True """ - s1 = sum([ ys[k] * hermite_coefficient(k, x, xs) - for k in range(0, len(xs)) ]) + s1 = sum( ys[k] * hermite_coefficient(k, x, xs) + for k in range(len(xs)) ) - s2 = sum([ y_primes[k] * hermite_deriv_coefficient(k, x, xs) - for k in range(0, len(xs)) ]) + s2 = sum( y_primes[k] * hermite_deriv_coefficient(k, x, xs) + for k in range(len(xs)) ) return (s1 + s2) @@ -347,4 +365,4 @@ def lagrange_psi(x, xs): """ - return product([ (x - xj) for xj in xs ]) + return product( (x - xj) for xj in xs )