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.

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 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() 

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 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

DifferentialOperatorAlgebra

is_singular(x0)[source]¶: Checks if the differential equation is singular at x0.

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 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

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 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]