# Operations on holonomic functions¶

Two holonomic functions can be added or multiplied with the result also a holonomic functions.

>>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
>>> from sympy.polys.domains import QQ
>>> from sympy import symbols
>>> x = symbols('x')
>>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')


p and q here are holonomic representation of $$e^x$$ and $$\sin(x)$$ respectively.

>>> p = HolonomicFunction(Dx - 1, x, 0, [1])
>>> q = HolonomicFunction(Dx**2 + 1, x, 0, [0, 1])


Holonomic representation of $$e^x+\sin(x)$$

>>> p + q
HolonomicFunction((-1) + (1)*Dx + (-1)*Dx**2 + (1)*Dx**3, x, 0, [1, 2, 1])


Holonomic representation of $$e^x \cdot \sin(x)$$

>>> p * q
HolonomicFunction((2) + (-2)*Dx + (1)*Dx**2, x, 0, [0, 1])


## Integration and Differentiation¶

HolonomicFunction.integrate(limits, initcond=False)[source]

Integrates the given holonomic function.

Examples

>>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
>>> from sympy.polys.domains import ZZ, QQ
>>> from sympy import symbols
>>> x = symbols('x')
>>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')
>>> HolonomicFunction(Dx - 1, x, 0, [1]).integrate((x, 0, x))  # e^x - 1
HolonomicFunction((-1)*Dx + (1)*Dx**2, x, 0, [0, 1])
>>> HolonomicFunction(Dx**2 + 1, x, 0, [1, 0]).integrate((x, 0, x))
HolonomicFunction((1)*Dx + (1)*Dx**3, x, 0, [0, 1, 0])

HolonomicFunction.diff(*args)[source]

Differentiation of the given Holonomic function.

Examples

>>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
>>> from sympy.polys.domains import ZZ, QQ
>>> from sympy import symbols
>>> x = symbols('x')
>>> R, Dx = DifferentialOperators(ZZ.old_poly_ring(x),'Dx')
>>> HolonomicFunction(Dx**2 + 1, x, 0, [0, 1]).diff().to_expr()
cos(x)
>>> HolonomicFunction(Dx - 2, x, 0, [1]).diff().to_expr()
2*exp(2*x)


## Composition with polynomials¶

HolonomicFunction.composition(expr, *args, **kwargs)[source]

Returns function after composition of a holonomic function with an algebraic function. The method can’t compute initial conditions for the result by itself, so they can be also be provided.

Examples

>>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
>>> from sympy.polys.domains import ZZ, QQ
>>> from sympy import symbols
>>> x = symbols('x')
>>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')
>>> HolonomicFunction(Dx - 1, x).composition(x**2, 0, [1])  # e^(x**2)
HolonomicFunction((-2*x) + (1)*Dx, x, 0, [1])
>>> HolonomicFunction(Dx**2 + 1, x).composition(x**2 - 1, 1, [1, 0])
HolonomicFunction((4*x**3) + (-1)*Dx + (x)*Dx**2, x, 1, [1, 0])


## Convert to holonomic sequence¶

HolonomicFunction.to_sequence(lb=True)[source]

Finds recurrence relation for the coefficients in the series expansion of the function about $$x_0$$, where $$x_0$$ is the point at which the initial condition is stored.

If the point $$x_0$$ is ordinary, solution of the form $$[(R, n_0)]$$ is returned. Where $$R$$ is the recurrence relation and $$n_0$$ is the smallest n for which the recurrence holds true.

If the point $$x_0$$ is regular singular, a list of solutions in the format $$(R, p, n_0)$$ is returned, i.e. $$[(R, p, n_0), ... ]$$. Each tuple in this vector represents a recurrence relation $$R$$ associated with a root of the indicial equation p. Conditions of a different format can also be provided in this case, see the docstring of HolonomicFunction class.

If it’s not possible to numerically compute a initial condition, it is returned as a symbol $$C_j$$, denoting the coefficient of $$(x - x_0)^j$$ in the power series about $$x_0$$.

References

Examples

>>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
>>> from sympy.polys.domains import ZZ, QQ
>>> from sympy import symbols, S
>>> x = symbols('x')
>>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')
>>> HolonomicFunction(Dx - 1, x, 0, [1]).to_sequence()
[(HolonomicSequence((-1) + (n + 1)Sn, n), u(0) = 1, 0)]
>>> HolonomicFunction((1 + x)*Dx**2 + Dx, x, 0, [0, 1]).to_sequence()
[(HolonomicSequence((n**2) + (n**2 + n)Sn, n), u(0) = 0, u(1) = 1, u(2) = -1/2, 2)]
>>> HolonomicFunction(-S(1)/2 + x*Dx, x, 0, {S(1)/2: [1]}).to_sequence()
[(HolonomicSequence((n), n), u(0) = 1, 1/2, 1)]


## Series expansion¶

HolonomicFunction.series(n=6, coefficient=False, order=True, _recur=None)[source]

Finds the power series expansion of given holonomic function about $$x_0$$.

A list of series might be returned if $$x_0$$ is a regular point with multiple roots of the indicial equation.

Examples

>>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
>>> from sympy.polys.domains import ZZ, QQ
>>> from sympy import symbols
>>> x = symbols('x')
>>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')
>>> HolonomicFunction(Dx - 1, x, 0, [1]).series()  # e^x
1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)
>>> HolonomicFunction(Dx**2 + 1, x, 0, [0, 1]).series(n=8)  # sin(x)
x - x**3/6 + x**5/120 - x**7/5040 + O(x**8)


## Numerical evaluation¶

HolonomicFunction.evalf(points, method='RK4', h=0.05, derivatives=False)[source]

Finds numerical value of a holonomic function using numerical methods. (RK4 by default). A set of points (real or complex) must be provided which will be the path for the numerical integration.

The path should be given as a list $$[x_1, x_2, ... x_n]$$. The numerical values will be computed at each point in this order $$x_1 --> x_2 --> x_3 ... --> x_n$$.

Returns values of the function at $$x_1, x_2, ... x_n$$ in a list.

Examples

>>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
>>> from sympy.polys.domains import ZZ, QQ
>>> from sympy import symbols
>>> x = symbols('x')
>>> R, Dx = DifferentialOperators(QQ.old_poly_ring(x),'Dx')


A straight line on the real axis from (0 to 1)

>>> r = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]


Runge-Kutta 4th order on e^x from 0.1 to 1. Exact solution at 1 is 2.71828182845905

>>> HolonomicFunction(Dx - 1, x, 0, [1]).evalf(r)
[1.10517083333333, 1.22140257085069, 1.34985849706254, 1.49182424008069,
1.64872063859684, 1.82211796209193, 2.01375162659678, 2.22553956329232,
2.45960141378007, 2.71827974413517]


Euler’s method for the same

>>> HolonomicFunction(Dx - 1, x, 0, [1]).evalf(r, method='Euler')
[1.1, 1.21, 1.331, 1.4641, 1.61051, 1.771561, 1.9487171, 2.14358881,
2.357947691, 2.5937424601]


One can also observe that the value obtained using Runge-Kutta 4th order is much more accurate than Euler’s method.

## Convert to a linear combination of hypergeometric functions¶

HolonomicFunction.to_hyper(as_list=False, _recur=None)[source]

Returns a hypergeometric function (or linear combination of them) representing the given holonomic function.

Returns an answer of the form: $$a_1 \cdot x^{b_1} \cdot{hyper()} + a_2 \cdot x^{b_2} \cdot{hyper()} ...$$

This is very useful as one can now use hyperexpand to find the symbolic expressions/functions.

Examples

>>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
>>> from sympy.polys.domains import ZZ, QQ
>>> from sympy import symbols
>>> x = symbols('x')
>>> R, Dx = DifferentialOperators(ZZ.old_poly_ring(x),'Dx')
>>> # sin(x)
>>> HolonomicFunction(Dx**2 + 1, x, 0, [0, 1]).to_hyper()
x*hyper((), (3/2,), -x**2/4)
>>> # exp(x)
>>> HolonomicFunction(Dx - 1, x, 0, [1]).to_hyper()
hyper((), (), x)


## Convert to a linear combination of Meijer G-functions¶

HolonomicFunction.to_meijerg()[source]

Returns a linear combination of Meijer G-functions.

Examples

>>> from sympy.holonomic import expr_to_holonomic
>>> from sympy import sin, cos, hyperexpand, log, symbols
>>> x = symbols('x')
>>> hyperexpand(expr_to_holonomic(cos(x) + sin(x)).to_meijerg())
sin(x) + cos(x)
>>> hyperexpand(expr_to_holonomic(log(x)).to_meijerg()).simplify()
log(x)


## Convert to expressions¶

HolonomicFunction.to_expr()[source]

Converts a Holonomic Function back to elementary functions.

Examples

>>> from sympy.holonomic.holonomic import HolonomicFunction, DifferentialOperators
>>> from sympy.polys.domains import ZZ, QQ
>>> from sympy import symbols, S
>>> x = symbols('x')
>>> R, Dx = DifferentialOperators(ZZ.old_poly_ring(x),'Dx')
>>> HolonomicFunction(x**2*Dx**2 + x*Dx + (x**2 - 1), x, 0, [0, S(1)/2]).to_expr()
besselj(1, x)
>>> HolonomicFunction((1 + x)*Dx**3 + Dx**2, x, 0, [1, 1, 1]).to_expr()
x*log(x + 1) + log(x + 1) + 1