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

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