Representation of holonomic functions in SymPy¶
Class DifferentialOperator
is used to represent the annihilator
but we create differential operators easily using the function
DifferentialOperators()
. Class HolonomicFunction
represents a holonomic function.
Let’s explain this with an example:
Take \(\sin(x)\) for instance, the differential equation satisfied by it is \(y^{(2)}(x) + y(x) = 0\). By definition we conclude it is a holonomic function. The general solution of this ODE is \(C_{1} \cdot \sin(x) + C_{2} \cdot \cos(x)\) but to get \(\sin(x)\) we need to provide initial conditions i.e. \(y(0) = 0, y^{(1)}(0) = 1\).
To represent the same in this module one needs to provide the differential
equation in the form of annihilator. Basically a differential operator is an
operator on functions that differentiates them. So \(D^{n} \cdot y(x) = y^{(n)}(x)\)
where \(y^{(n)}(x)\) denotes n
times differentiation of \(y(x)\) with
respect to x
.
So the differential equation can also be written as \(D^{2} \cdot y(x) + y(x) = 0\) or \((D^{2} + 1) \cdot y(x) = 0\). The part left of \(y(x)\) is the annihilator i.e. \(D^{2}+1\).
So this is how one will represent \(\sin(x)\) as a Holonomic Function:
>>> from sympy.holonomic import DifferentialOperators, HolonomicFunction
>>> from sympy.abc import x
>>> from sympy import ZZ
>>> R, D = DifferentialOperators(ZZ.old_poly_ring(x), 'D')
>>> HolonomicFunction(D**2 + 1, x, 0, [0, 1])
HolonomicFunction((1) + (1)*D**2, x, 0, [0, 1])
The polynomial coefficients will be members of the ring ZZ[x]
in the example.
The D
operator returned by the function DifferentialOperators()
can
be used to create annihilators just like SymPy expressions.
We currently use the older implementations of rings in SymPy for priority
mechanism.
HolonomicFunction¶
- class sympy.holonomic.holonomic.HolonomicFunction(annihilator, x, x0=0, y0=None)[source]¶
A Holonomic Function is a solution to a linear homogeneous ordinary differential equation with polynomial coefficients. This differential equation can also be represented by an annihilator i.e. a Differential Operator
L
such that \(L.f = 0\). For uniqueness of these functions, initial conditions can also be provided along with the annihilator.Explanation
Holonomic functions have closure properties and thus forms a ring. Given two Holonomic Functions f and g, their sum, product, integral and derivative is also a Holonomic Function.
For ordinary points initial condition should be a vector of values of the derivatives i.e. \([y(x_0), y'(x_0), y''(x_0) ... ]\).
For regular singular points initial conditions can also be provided in this format: \({s0: [C_0, C_1, ...], s1: [C^1_0, C^1_1, ...], ...}\) where s0, s1, … are the roots of indicial equation and vectors \([C_0, C_1, ...], [C^0_0, C^0_1, ...], ...\) are the corresponding initial terms of the associated power series. See Examples below.
Examples
>>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators >>> from sympy.polys.domains import QQ >>> from sympy import symbols, S >>> x = symbols('x') >>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')
>>> p = HolonomicFunction(Dx - 1, x, 0, [1]) # e^x >>> q = HolonomicFunction(Dx**2 + 1, x, 0, [0, 1]) # sin(x)
>>> p + q # annihilator of e^x + sin(x) HolonomicFunction((-1) + (1)*Dx + (-1)*Dx**2 + (1)*Dx**3, x, 0, [1, 2, 1])
>>> p * q # annihilator of e^x * sin(x) HolonomicFunction((2) + (-2)*Dx + (1)*Dx**2, x, 0, [0, 1])
An example of initial conditions for regular singular points, the indicial equation has only one root \(1/2\).
>>> HolonomicFunction(-S(1)/2 + x*Dx, x, 0, {S(1)/2: [1]}) HolonomicFunction((-1/2) + (x)*Dx, x, 0, {1/2: [1]})
>>> HolonomicFunction(-S(1)/2 + x*Dx, x, 0, {S(1)/2: [1]}).to_expr() sqrt(x)
To plot a Holonomic Function, one can use \(.evalf()\) for numerical computation. Here’s an example on \(sin(x)**2/x\) using numpy and matplotlib.
>>> import sympy.holonomic >>> from sympy import var, sin >>> import matplotlib.pyplot as plt >>> import numpy as np >>> var("x") >>> r = np.linspace(1, 5, 100) >>> y = sympy.holonomic.expr_to_holonomic(sin(x)**2/x, x0=1).evalf(r) >>> plt.plot(r, y, label="holonomic function") >>> plt.show()
DifferentialOperator¶
- class sympy.holonomic.holonomic.DifferentialOperator(list_of_poly, parent)[source]¶
Differential Operators are elements of Weyl Algebra. The Operators are defined by a list of polynomials in the base ring and the parent ring of the Operator i.e. the algebra it belongs to.
Explanation
Takes a list of polynomials for each power of
Dx
and the parent ring which must be an instance of DifferentialOperatorAlgebra.A Differential Operator can be created easily using the operator
Dx
. See examples below.Examples
>>> from sympy.holonomic.holonomic import DifferentialOperator, DifferentialOperators >>> from sympy.polys.domains import ZZ >>> from sympy import symbols >>> x = symbols('x') >>> R, Dx = DifferentialOperators(ZZ.old_poly_ring(x),'Dx')
>>> DifferentialOperator([0, 1, x**2], R) (1)*Dx + (x**2)*Dx**2
>>> (x*Dx*x + 1 - Dx**2)**2 (2*x**2 + 2*x + 1) + (4*x**3 + 2*x**2 - 4)*Dx + (x**4 - 6*x - 2)*Dx**2 + (-2*x**2)*Dx**3 + (1)*Dx**4
See also
DifferentialOperators¶
- sympy.holonomic.holonomic.DifferentialOperators(base, generator)[source]¶
This function is used to create annihilators using
Dx
.- Parameters
base:
Base polynomial ring for the algebra. The base polynomial ring is the ring of polynomials in \(x\) that will appear as coefficients in the operators.
generator:
Generator of the algebra which can be either a noncommutative
Symbol
or a string. e.g. “Dx” or “D”.
Explanation
Returns an Algebra of Differential Operators also called Weyl Algebra and the operator for differentiation i.e. the
Dx
operator.Examples
>>> from sympy.polys.domains import ZZ >>> from sympy.abc import x >>> from sympy.holonomic.holonomic import DifferentialOperators >>> R, Dx = DifferentialOperators(ZZ.old_poly_ring(x), 'Dx') >>> R Univariate Differential Operator Algebra in intermediate Dx over the base ring ZZ[x] >>> Dx*x (1) + (x)*Dx
DifferentialOperatorAlgebra¶
- class sympy.holonomic.holonomic.DifferentialOperatorAlgebra(base, generator)[source]¶
An Ore Algebra is a set of noncommutative polynomials in the intermediate
Dx
and coefficients in a base polynomial ring \(A\). It follows the commutation rule:\[Dxa = \sigma(a)Dx + \delta(a)\]for \(a \subset A\).
Where \(\sigma: A \Rightarrow A\) is an endomorphism and \(\delta: A \rightarrow A\) is a skew-derivation i.e. \(\delta(ab) = \delta(a) b + \sigma(a) \delta(b)\).
If one takes the sigma as identity map and delta as the standard derivation then it becomes the algebra of Differential Operators also called a Weyl Algebra i.e. an algebra whose elements are Differential Operators.
This class represents a Weyl Algebra and serves as the parent ring for Differential Operators.
Examples
>>> from sympy.polys.domains import ZZ >>> from sympy import symbols >>> from sympy.holonomic.holonomic import DifferentialOperators >>> x = symbols('x') >>> R, Dx = DifferentialOperators(ZZ.old_poly_ring(x), 'Dx') >>> R Univariate Differential Operator Algebra in intermediate Dx over the base ring ZZ[x]
See also