|
| 1 | +"""Jacobi polynomials.""" |
| 2 | + |
| 3 | +import typing |
| 4 | +from functools import cache |
| 5 | + |
| 6 | +from math import prod, factorial |
| 7 | +import sympy |
| 8 | + |
| 9 | +from symfem.functions import ScalarFunction |
| 10 | +from symfem.symbols import x |
| 11 | + |
| 12 | +__all__: typing.List[str] = [] |
| 13 | + |
| 14 | + |
| 15 | +def _jrc( |
| 16 | + a: int, b: int, n: int |
| 17 | +) -> typing.Tuple[ |
| 18 | + sympy.core.expr.Expr, |
| 19 | + sympy.core.expr.Expr, |
| 20 | + sympy.core.expr.Expr, |
| 21 | +]: |
| 22 | + """Get the Jacobi recurrence relation coefficients. |
| 23 | +
|
| 24 | + Args: |
| 25 | + a: The parameter a |
| 26 | + b: The parameter b |
| 27 | + n: The parameter n |
| 28 | +
|
| 29 | + Returns: |
| 30 | + The Jacobi coefficients |
| 31 | + """ |
| 32 | + if b != 0: |
| 33 | + raise NotImplementedError() |
| 34 | + return ( |
| 35 | + sympy.Rational((a + 2 * n + 1) * (a + 2 * n + 2), 2 * (n + 1) * (a + n + 1)), |
| 36 | + sympy.Rational(a * a * (a + 2 * n + 1), 2 * (n + 1) * (a + n + 1) * (a + 2 * n)), |
| 37 | + sympy.Rational(n * (a + n) * (a + 2 * n + 2), (n + 1) * (a + n + 1) * (a + 2 * n)), |
| 38 | + ) |
| 39 | + |
| 40 | + |
| 41 | +def _jrc_monic( |
| 42 | + a: int, b: int, n: int |
| 43 | +) -> typing.Tuple[ |
| 44 | + sympy.core.expr.Expr, |
| 45 | + sympy.core.expr.Expr, |
| 46 | +]: |
| 47 | + """Get the Jacobi recurrence relation coefficients for monic polynomials. |
| 48 | +
|
| 49 | + Args: |
| 50 | + a: The parameter a |
| 51 | + b: The parameter b |
| 52 | + n: The parameter n |
| 53 | +
|
| 54 | + Returns: |
| 55 | + The Jacobi coefficients |
| 56 | + """ |
| 57 | + return ( |
| 58 | + sympy.Rational(a**2 - b**2, (2 * n + a + b) * (2 * n + a + b + 2)), |
| 59 | + sympy.Rational( |
| 60 | + 4 * n * (n + a) * (n + b) * (n + a + b), |
| 61 | + (2 * n + a + b + 1) * (2 * n + a + b) ** 2 * (2 * n + a + b - 1), |
| 62 | + ), |
| 63 | + ) |
| 64 | + |
| 65 | + |
| 66 | +@cache |
| 67 | +def jacobi_polynomial_x(n: int, a: int, b: int) -> sympy.core.expr.Expr: |
| 68 | + """Get a Jacobi polynomial. |
| 69 | +
|
| 70 | + Args: |
| 71 | + n: Polynomial degree |
| 72 | + a: The parameter a |
| 73 | + b: The parameter b |
| 74 | + """ |
| 75 | + if a < 0: |
| 76 | + if b <= 0 and n + a >= 0: |
| 77 | + # See https://mathoverflow.net/questions/298661/jacobi-polynomials-with-negative-integer-parameters |
| 78 | + return ( |
| 79 | + ((x[0] - 1) / 2) ** -a |
| 80 | + * ((x[0] + 1) / 2) ** -b |
| 81 | + * jacobi_polynomial(n + a + b, -a, -b) |
| 82 | + ) |
| 83 | + raise NotImplementedError() |
| 84 | + if n == 0: |
| 85 | + return sympy.Integer(1) |
| 86 | + if n == 1: |
| 87 | + return (a + 1) + (a + b + 2) * (x[0] - 1) / 2 |
| 88 | + i, j, k = _jrc(a, b, n - 1) |
| 89 | + return (i * x[0] + j) * jacobi_polynomial(n - 1, a, b) - k * jacobi_polynomial(n - 2, a, b) |
| 90 | + |
| 91 | + |
| 92 | +def jacobi_polynomial(n: int, a: int, b: int, variable: sympy.Symbol = x[0]) -> ScalarFunction: |
| 93 | + """Get a Jacobi polynomial. |
| 94 | +
|
| 95 | + Args: |
| 96 | + n: Polynomial degree |
| 97 | + a: The parameter a |
| 98 | + b: The parameter b |
| 99 | + variable: Variable to use |
| 100 | + """ |
| 101 | + return ScalarFunction(jacobi_polynomial_x(n, a, b).subs(x[0], variable)) |
| 102 | + |
| 103 | + |
| 104 | +def choose(n: int, r: int) -> int: |
| 105 | + """Choose function.""" |
| 106 | + return prod(range(n - r + 1, n + 1)) // factorial(r) |
| 107 | + |
| 108 | + |
| 109 | +@cache |
| 110 | +def monic_jacobi_polynomial_x(n: int, a: int, b: int) -> sympy.core.expr.Expr: |
| 111 | + """Get a monic Jacobi polynomial. |
| 112 | +
|
| 113 | + Args: |
| 114 | + n: Polynomial degree |
| 115 | + a: The parameter a |
| 116 | + b: The parameter b |
| 117 | + """ |
| 118 | + if n == 0: |
| 119 | + return sympy.Integer(1) |
| 120 | + if a + b == -n: |
| 121 | + return sum( |
| 122 | + choose(n + a, i) * choose(n + b, n - i) * (x[0] - 1) ** (n - i) * (x[0] + 1) ** i |
| 123 | + for i in range(n + 1) |
| 124 | + ) / choose(2 * n + a + b, n) |
| 125 | + if n == 1: |
| 126 | + return x[0] + sympy.Rational(2 * (a + 1), a + b + 2) - 1 |
| 127 | + i, j = _jrc_monic(a, b, n - 1) |
| 128 | + return (x[0] + i) * monic_jacobi_polynomial(n - 1, a, b) - j * monic_jacobi_polynomial( |
| 129 | + n - 2, a, b |
| 130 | + ) |
| 131 | + |
| 132 | + |
| 133 | +def monic_jacobi_polynomial( |
| 134 | + n: int, a: int, b: int, variable: sympy.Symbol = x[0] |
| 135 | +) -> ScalarFunction: |
| 136 | + """Get a monic Jacobi polynomial. |
| 137 | +
|
| 138 | + Args: |
| 139 | + n: Polynomial degree |
| 140 | + a: The parameter a |
| 141 | + b: The parameter b |
| 142 | + variable: Variable to use |
| 143 | + """ |
| 144 | + return ScalarFunction(monic_jacobi_polynomial_x(n, a, b).subs(x[0], variable)) |
0 commit comments