Symbolic Computation#

AUTHORS:

Bobby Moretti and William Stein (2006-2007)
Robert Bradshaw (2007-10): minpoly(), numerical algorithm
Robert Bradshaw (2008-10): minpoly(), algebraic algorithm
Golam Mortuza Hossain (2009-06-15): _limit_latex()
Golam Mortuza Hossain (2009-06-22): _laplace_latex(), _inverse_laplace_latex()
Tom Coates (2010-06-11): fixed github issue #9217

EXAMPLES:

The basic units of the calculus package are symbolic expressions which are elements of the symbolic expression ring (SR). To create a symbolic variable object in Sage, use the var() function, whose argument is the text of that variable. Note that Sage is intelligent about LaTeXing variable names.

sage: x1 = var('x1'); x1
x1
sage: latex(x1)
x_{1}
sage: theta = var('theta'); theta
theta
sage: latex(theta)
\theta

Sage predefines x to be a global indeterminate. Thus the following works:

sage: x^2
x^2
sage: type(x)
<class 'sage.symbolic.expression.Expression'>

More complicated expressions in Sage can be built up using ordinary arithmetic. The following are valid, and follow the rules of Python arithmetic: (The ‘=’ operator represents assignment, and not equality)

sage: var('x,y,z')
(x, y, z)
sage: f = x + y + z/(2*sin(y*z/55))
sage: g = f^f; g
(x + y + 1/2*z/sin(1/55*y*z))^(x + y + 1/2*z/sin(1/55*y*z))

Differentiation and integration are available, but behind the scenes through Maxima:

sage: f = sin(x)/cos(2*y)
sage: f.derivative(y)
2*sin(x)*sin(2*y)/cos(2*y)^2
sage: g = f.integral(x); g
-cos(x)/cos(2*y)

Note that these methods usually require an explicit variable name. If none is given, Sage will try to find one for you.

sage: f = sin(x); f.derivative()
cos(x)

If the expression is a callable symbolic expression (i.e., the variable order is specified), then Sage can calculate the matrix derivative (i.e., the gradient, Jacobian matrix, etc.) if no variables are specified. In the example below, we use the second derivative test to determine that there is a saddle point at (0,-1/2).

sage: f(x,y) = x^2*y + y^2 + y
sage: f.diff()  # gradient
(x, y) |--> (2*x*y, x^2 + 2*y + 1)
sage: solve(list(f.diff()), [x,y])
[[x == -I, y == 0], [x == I, y == 0], [x == 0, y == (-1/2)]]
sage: H=f.diff(2); H  # Hessian matrix
[(x, y) |--> 2*y (x, y) |--> 2*x]
[(x, y) |--> 2*x   (x, y) |--> 2]
sage: H(x=0, y=-1/2)
[-1  0]
[ 0  2]
sage: H(x=0, y=-1/2).eigenvalues()
[-1, 2]

Here we calculate the Jacobian for the polar coordinate transformation:

sage: T(r,theta) = [r*cos(theta),r*sin(theta)]
sage: T
(r, theta) |--> (r*cos(theta), r*sin(theta))
sage: T.diff() # Jacobian matrix
[   (r, theta) |--> cos(theta) (r, theta) |--> -r*sin(theta)]
[   (r, theta) |--> sin(theta)  (r, theta) |--> r*cos(theta)]
sage: diff(T) # Jacobian matrix
[   (r, theta) |--> cos(theta) (r, theta) |--> -r*sin(theta)]
[   (r, theta) |--> sin(theta)  (r, theta) |--> r*cos(theta)]
sage: T.diff().det() # Jacobian
(r, theta) |--> r*cos(theta)^2 + r*sin(theta)^2

When the order of variables is ambiguous, Sage will raise an exception when differentiating:

sage: f = sin(x+y); f.derivative()
Traceback (most recent call last):
...
ValueError: No differentiation variable specified.

Simplifying symbolic sums is also possible, using the sum() command, which also uses Maxima in the background:

sage: k, m = var('k, m')
sage: sum(1/k^4, k, 1, oo)
1/90*pi^4
sage: sum(binomial(m,k), k, 0, m)
2^m

Symbolic matrices can be used as well in various ways, including exponentiation:

sage: M = matrix([[x,x^2],[1/x,x]])
sage: M^2
[x^2 + x   2*x^3]
[      2 x^2 + x]
sage: e^M
[          1/2*(e^(2*sqrt(x)) + 1)*e^(x - sqrt(x)) 1/2*(x*e^(2*sqrt(x)) - x)*sqrt(x)*e^(x - sqrt(x))]
[  1/2*(e^(2*sqrt(x)) - 1)*e^(x - sqrt(x))/x^(3/2)           1/2*(e^(2*sqrt(x)) + 1)*e^(x - sqrt(x))]

Complex exponentiation works, but may require a patched version of maxima (github issue #32898) for now:

sage: M = i*matrix([[pi]])
sage: e^M  # not tested, requires patched maxima
[-1]
sage: M = i*matrix([[pi,0],[0,2*pi]])
sage: e^M
[-1  0]
[ 0  1]
sage: M = matrix([[0,pi],[-pi,0]])
sage: e^M
[-1  0]
[ 0 -1]

Substitution works similarly. We can substitute with a python dict:

sage: f = sin(x*y - z)
sage: f({x: var('t'), y: z})
sin(t*z - z)

Also we can substitute with keywords:

sage: f = sin(x*y - z)
sage: f(x=t, y=z)
sin(t*z - z)

Another example:

sage: f = sin(2*pi*x/y)
sage: f(x=4)
sin(8*pi/y)

It is no longer allowed to call expressions with positional arguments:

sage: f = sin(x)
sage: f(y)
Traceback (most recent call last):
...
TypeError: Substitution using function-call syntax and unnamed
arguments has been removed. You can use named arguments instead, like
EXPR(x=..., y=...)
sage: f(x=pi)
0

We can also make a CallableSymbolicExpression, which is a SymbolicExpression that is a function of specified variables in a fixed order. Each SymbolicExpression has a function(...) method that is used to create a CallableSymbolicExpression, as illustrated below:

sage: u = log((2-x)/(y+5))
sage: f = u.function(x, y); f
(x, y) |--> log(-(x - 2)/(y + 5))

There is an easier way of creating a CallableSymbolicExpression, which relies on the Sage preparser.

sage: f(x,y) = log(x)*cos(y); f
(x, y) |--> cos(y)*log(x)

Then we have fixed an order of variables and there is no ambiguity substituting or evaluating:

sage: f(x,y) = log((2-x)/(y+5))
sage: f(7,t)
log(-5/(t + 5))

Some further examples:

sage: f = 5*sin(x)
sage: f
5*sin(x)
sage: f(x=2)
5*sin(2)
sage: f(x=pi)
0
sage: float(f(x=pi))
0.0

Another example:

sage: f = integrate(1/sqrt(9+x^2), x); f
arcsinh(1/3*x)
sage: f(x=3)
arcsinh(1)
sage: f.derivative(x)
1/sqrt(x^2 + 9)

We compute the length of the parabola from 0 to 2:

sage: x = var('x')
sage: y = x^2
sage: dy = derivative(y,x)
sage: z = integral(sqrt(1 + dy^2), x, 0, 2)
sage: z
sqrt(17) + 1/4*arcsinh(4)
sage: n(z,200)
4.6467837624329358733826155674904591885104869874232887508703
sage: float(z)
4.646783762432936

We test pickling:

sage: x, y = var('x,y')
sage: f = -sqrt(pi)*(x^3 + sin(x/cos(y)))
sage: bool(loads(dumps(f)) == f)
True

Coercion examples:

We coerce various symbolic expressions into the complex numbers:

sage: CC(I)
1.00000000000000*I
sage: CC(2*I)
2.00000000000000*I
sage: ComplexField(200)(2*I)
2.0000000000000000000000000000000000000000000000000000000000*I
sage: ComplexField(200)(sin(I))
1.1752011936438014568823818505956008151557179813340958702296*I
sage: f = sin(I) + cos(I/2); f
cosh(1/2) + I*sinh(1)
sage: CC(f)
1.12762596520638 + 1.17520119364380*I
sage: ComplexField(200)(f)
1.1276259652063807852262251614026720125478471180986674836290
 + 1.1752011936438014568823818505956008151557179813340958702296*I
sage: ComplexField(100)(f)
1.1276259652063807852262251614 + 1.1752011936438014568823818506*I

We illustrate construction of an inverse sum where each denominator has a new variable name:

sage: f = sum(1/var('n%s'%i)^i for i in range(10))
sage: f
1/n1 + 1/n2^2 + 1/n3^3 + 1/n4^4 + 1/n5^5 + 1/n6^6 + 1/n7^7 + 1/n8^8 + 1/n9^9 + 1

Note that after calling var, the variables are immediately available for use:

sage: (n1 + n2)^5
(n1 + n2)^5

We can, of course, substitute:

sage: f(n9=9, n7=n6)
1/n1 + 1/n2^2 + 1/n3^3 + 1/n4^4 + 1/n5^5 + 1/n6^6 + 1/n6^7 + 1/n8^8
 + 387420490/387420489

sage.calculus.calculus.at(ex, *args, **kwds)#

Parses at formulations from other systems, such as Maxima. Replaces evaluation ‘at’ a point with substitution method of a symbolic expression.

EXAMPLES:

We do not import at at the top level, but we can use it as a synonym for substitution if we import it:

sage: g = x^3 - 3
sage: from sage.calculus.calculus import at
sage: at(g, x=1)
-2
sage: g.subs(x=1)
-2

We find a formal Taylor expansion:

sage: h,x = var('h,x')
sage: u = function('u')
sage: u(x + h)
u(h + x)
sage: diff(u(x+h), x)
D[0](u)(h + x)
sage: taylor(u(x+h), h, 0, 4)
1/24*h^4*diff(u(x), x, x, x, x) + 1/6*h^3*diff(u(x), x, x, x)
 + 1/2*h^2*diff(u(x), x, x) + h*diff(u(x), x) + u(x)

We compute a Laplace transform:

sage: var('s,t')
(s, t)
sage: f = function('f')(t)
sage: f.diff(t, 2)
diff(f(t), t, t)
sage: f.diff(t,2).laplace(t,s)
s^2*laplace(f(t), t, s) - s*f(0) - D[0](f)(0)

We can also accept a non-keyword list of expression substitutions, like Maxima does (github issue #12796):

sage: from sage.calculus.calculus import at
sage: f = function('f')
sage: at(f(x), [x == 1])
f(1)

sage.calculus.calculus.dummy_diff(*args)#

This function is called when ‘diff’ appears in a Maxima string.

EXAMPLES:

sage: from sage.calculus.calculus import dummy_diff
sage: x,y = var('x,y')
sage: dummy_diff(sin(x*y), x, SR(2), y, SR(1))
-x*y^2*cos(x*y) - 2*y*sin(x*y)

Here the function is used implicitly:

sage: a = var('a')
sage: f = function('cr')(a)
sage: g = f.diff(a); g
diff(cr(a), a)

sage.calculus.calculus.dummy_integrate(*args)#

This function is called to create formal wrappers of integrals that Maxima can’t compute:

EXAMPLES:

sage: from sage.calculus.calculus import dummy_integrate
sage: f = function('f')
sage: dummy_integrate(f(x), x)
integrate(f(x), x)
sage: a,b = var('a,b')
sage: dummy_integrate(f(x), x, a, b)
integrate(f(x), x, a, b)

sage.calculus.calculus.dummy_inverse_laplace(*args)#

This function is called to create formal wrappers of inverse Laplace transforms that Maxima can’t compute:

EXAMPLES:

sage: from sage.calculus.calculus import dummy_inverse_laplace
sage: s,t = var('s,t')
sage: F = function('F')
sage: dummy_inverse_laplace(F(s), s, t)
ilt(F(s), s, t)

sage.calculus.calculus.dummy_laplace(*args)#

This function is called to create formal wrappers of laplace transforms that Maxima can’t compute:

EXAMPLES:

sage: from sage.calculus.calculus import dummy_laplace
sage: s,t = var('s,t')
sage: f = function('f')
sage: dummy_laplace(f(t), t, s)
laplace(f(t), t, s)

sage.calculus.calculus.dummy_pochhammer(*args)#

This function is called to create formal wrappers of Pochhammer symbols

EXAMPLES:

sage: from sage.calculus.calculus import dummy_pochhammer
sage: s,t = var('s,t')
sage: dummy_pochhammer(s, t)
gamma(s + t)/gamma(s)

sage.calculus.calculus.inverse_laplace(ex, s, t, algorithm='maxima')#

Return the inverse Laplace transform with respect to the variable \(t\) and transform parameter \(s\), if possible.

If this function cannot find a solution, a formal function is returned. The function that is returned may be viewed as a function of \(t\).

DEFINITION:

The inverse Laplace transform of a function \(F(s)\) is the function \(f(t)\), defined by

\[f(t) = \frac{1}{2\pi i} \int_{\gamma-i\infty}^{\gamma + i\infty} e^{st} F(s) ds,\]

where \(\gamma\) is chosen so that the contour path of integration is in the region of convergence of \(F(s)\).

INPUT:

ex – a symbolic expression
s – transform parameter
t – independent variable
algorithm – (default: 'maxima') one of
- 'maxima' – use Maxima (the default)
- 'sympy' – use SymPy
- 'giac' – use Giac

See also

inverse_laplace()

EXAMPLES:

We compute a few Laplace transforms:

sage: var('x, s, z, t, t0')
(x, s, z, t, t0)
sage: sin(x).laplace(x, s)
1/(s^2 + 1)
sage: (z + exp(x)).laplace(x, s)
z/s + 1/(s - 1)
sage: log(t/t0).laplace(t, s)
-(euler_gamma + log(s) + log(t0))/s

We do a formal calculation:

sage: f = function('f')(x)
sage: g = f.diff(x); g
diff(f(x), x)
sage: g.laplace(x, s)
s*laplace(f(x), x, s) - f(0)

A BATTLE BETWEEN the X-women and the Y-men (by David Joyner): Solve

\[x' = -16y, x(0)=270, y' = -x + 1, y(0) = 90.\]

This models a fight between two sides, the “X-women” and the “Y-men”, where the X-women have 270 initially and the Y-men have 90, but the Y-men are better at fighting, because of the higher factor of “-16” vs “-1”, and also get an occasional reinforcement, because of the “+1” term.

sage: var('t')
t
sage: t = var('t')
sage: x = function('x')(t)
sage: y = function('y')(t)
sage: de1 = x.diff(t) + 16*y
sage: de2 = y.diff(t) + x - 1
sage: de1.laplace(t, s)
s*laplace(x(t), t, s) + 16*laplace(y(t), t, s) - x(0)
sage: de2.laplace(t, s)
s*laplace(y(t), t, s) - 1/s + laplace(x(t), t, s) - y(0)

Next we form the augmented matrix of the above system:

sage: A = matrix([[s, 16, 270], [1, s, 90+1/s]])
sage: E = A.echelon_form()
sage: xt = E[0,2].inverse_laplace(s,t)
sage: yt = E[1,2].inverse_laplace(s,t)
sage: xt
-91/2*e^(4*t) + 629/2*e^(-4*t) + 1
sage: yt
91/8*e^(4*t) + 629/8*e^(-4*t)
sage: p1 = plot(xt, 0, 1/2, rgbcolor=(1,0,0))                                   # needs sage.plot
sage: p2 = plot(yt, 0, 1/2, rgbcolor=(0,1,0))                                   # needs sage.plot
sage: import tempfile
sage: with tempfile.NamedTemporaryFile(suffix=".png") as f:                     # needs sage.plot
....:     (p1 + p2).save(f.name)

Another example:

sage: var('a,s,t')
(a, s, t)
sage: f = exp (2*t + a) * sin(t) * t; f
t*e^(a + 2*t)*sin(t)
sage: L = laplace(f, t, s); L
2*(s - 2)*e^a/(s^2 - 4*s + 5)^2
sage: inverse_laplace(L, s, t)
t*e^(a + 2*t)*sin(t)

The Laplace transform of the exponential function:

sage: laplace(exp(x), x, s)
1/(s - 1)

Dirac’s delta distribution is handled (the output of SymPy is related to a choice that has to be made when defining Laplace transforms of distributions):

sage: laplace(dirac_delta(t), t, s)
1
sage: F, a, cond = laplace(dirac_delta(t), t, s, algorithm='sympy')
sage: a, cond  # random - sympy <1.10 gives (-oo, True)
(0, True)
sage: F        # random - sympy <1.9 includes undefined heaviside(0) in answer
1
sage: laplace(dirac_delta(t), t, s, algorithm='giac')
1

Heaviside step function can be handled with different interfaces. Try with Maxima:

sage: laplace(heaviside(t-1), t, s)
e^(-s)/s

Try with giac:

sage: laplace(heaviside(t-1), t, s, algorithm='giac')
e^(-s)/s

Try with SymPy:

sage: laplace(heaviside(t-1), t, s, algorithm='sympy')
(e^(-s)/s, 0, True)

sage.calculus.calculus.lim(ex, dir=None, taylor=False, algorithm='maxima', **argv)#

Return the limit as the variable \(v\) approaches \(a\) from the given direction.

expr.limit(x = a)
expr.limit(x = a, dir='+')

INPUT:

dir – (default: None); may have the value 'plus' (or '+' or 'right' or 'above') for a limit from above, 'minus' (or '-' or 'left' or 'below') for a limit from below, or may be omitted (implying a two-sided limit is to be computed).
taylor – (default: False); if True, use Taylor series, which allows more limits to be computed (but may also crash in some obscure cases due to bugs in Maxima).
**argv - 1 named parameter

Note

The output may also use und (undefined), ind (indefinite but bounded), and infinity (complex infinity).

EXAMPLES:

sage: x = var('x')
sage: f = (1 + 1/x)^x
sage: f.limit(x=oo)
e
sage: f.limit(x=5)
7776/3125

Domain to real, a regression in 5.46.0, see https://sf.net/p/maxima/bugs/4138

sage: maxima_calculus.eval("domain:real")
...
sage: f.limit(x=1.2).n()
2.06961575467...
sage: maxima_calculus.eval("domain:complex");
...

Otherwise, it works

sage: f.limit(x=I, taylor=True)
(-I + 1)^I
sage: f(x=1.2)
2.0696157546720...
sage: f(x=I)
(-I + 1)^I
sage: CDF(f(x=I))
2.0628722350809046 + 0.7450070621797239*I
sage: CDF(f.limit(x=I))
2.0628722350809046 + 0.7450070621797239*I

Notice that Maxima may ask for more information:

sage: var('a')
a
sage: limit(x^a,x=0)
Traceback (most recent call last):
...
ValueError: Computation failed since Maxima requested additional
constraints; using the 'assume' command before evaluation
*may* help (example of legal syntax is 'assume(a>0)', see
`assume?` for more details)
Is a positive, negative or zero?

With this example, Maxima is looking for a LOT of information:

sage: assume(a>0)
sage: limit(x^a,x=0)  # random - maxima 5.46.0 does not need extra assumption
Traceback (most recent call last):
...
ValueError: Computation failed since Maxima requested additional
constraints; using the 'assume' command before evaluation *may* help
(example of legal syntax is 'assume(a>0)', see `assume?` for
 more details)
Is a an integer?
sage: assume(a,'integer')
sage: limit(x^a, x=0)  # random - maxima 5.46.0 does not need extra assumption
Traceback (most recent call last):
...
ValueError: Computation failed since Maxima requested additional
constraints; using the 'assume' command before evaluation *may* help
(example of legal syntax is 'assume(a>0)', see `assume?` for
 more details)
Is a an even number?
sage: assume(a, 'even')
sage: limit(x^a, x=0)
0
sage: forget()

More examples:

sage: limit(x*log(x), x=0, dir='+')
0
sage: lim((x+1)^(1/x), x=0)
e
sage: lim(e^x/x, x=oo)
+Infinity
sage: lim(e^x/x, x=-oo)
0
sage: lim(-e^x/x, x=oo)
-Infinity
sage: lim((cos(x))/(x^2), x=0)
+Infinity
sage: lim(sqrt(x^2+1) - x, x=oo)
0
sage: lim(x^2/(sec(x)-1), x=0)
2
sage: lim(cos(x)/(cos(x)-1), x=0)
-Infinity
sage: lim(x*sin(1/x), x=0)
0
sage: limit(e^(-1/x), x=0, dir='right')
0
sage: limit(e^(-1/x), x=0, dir='left')
+Infinity

sage: f = log(log(x)) / log(x)
sage: forget(); assume(x < -2); lim(f, x=0, taylor=True)
0
sage: forget()

Here ind means “indefinite but bounded”:

sage: lim(sin(1/x), x = 0)
ind

We can use other packages than maxima, namely “sympy”, “giac”, “fricas”.

With the standard package Giac:

sage: from sage.libs.giac.giac import libgiac     # random
sage: (exp(-x)/(2+sin(x))).limit(x=oo, algorithm='giac')
0
sage: limit(e^(-1/x), x=0, dir='right', algorithm='giac')
0
sage: limit(e^(-1/x), x=0, dir='left', algorithm='giac')
+Infinity
sage: (x / (x+2^x+cos(x))).limit(x=-infinity, algorithm='giac')
1

With the optional package FriCAS:

sage: (x / (x+2^x+cos(x))).limit(x=-infinity, algorithm='fricas')       # optional - fricas
1
sage: limit(e^(-1/x), x=0, dir='right', algorithm='fricas')             # optional - fricas
0
sage: limit(e^(-1/x), x=0, dir='left', algorithm='fricas')              # optional - fricas
+Infinity

One can also call Mathematica’s online interface:

sage: limit(pi+log(x)/x,x=oo, algorithm='mathematica_free') # optional - internet
pi

sage.calculus.calculus.limit(ex, dir=None, taylor=False, algorithm='maxima', **argv)#

Return the limit as the variable \(v\) approaches \(a\) from the given direction.

expr.limit(x = a)
expr.limit(x = a, dir='+')

INPUT:

dir – (default: None); may have the value 'plus' (or '+' or 'right' or 'above') for a limit from above, 'minus' (or '-' or 'left' or 'below') for a limit from below, or may be omitted (implying a two-sided limit is to be computed).
taylor – (default: False); if True, use Taylor series, which allows more limits to be computed (but may also crash in some obscure cases due to bugs in Maxima).
**argv - 1 named parameter

Note

The output may also use und (undefined), ind (indefinite but bounded), and infinity (complex infinity).

EXAMPLES:

sage: x = var('x')
sage: f = (1 + 1/x)^x
sage: f.limit(x=oo)
e
sage: f.limit(x=5)
7776/3125

Domain to real, a regression in 5.46.0, see https://sf.net/p/maxima/bugs/4138

sage: maxima_calculus.eval("domain:real")
...
sage: f.limit(x=1.2).n()
2.06961575467...
sage: maxima_calculus.eval("domain:complex");
...

Otherwise, it works

sage: f.limit(x=I, taylor=True)
(-I + 1)^I
sage: f(x=1.2)
2.0696157546720...
sage: f(x=I)
(-I + 1)^I
sage: CDF(f(x=I))
2.0628722350809046 + 0.7450070621797239*I
sage: CDF(f.limit(x=I))
2.0628722350809046 + 0.7450070621797239*I

Notice that Maxima may ask for more information:

sage: var('a')
a
sage: limit(x^a,x=0)
Traceback (most recent call last):
...
ValueError: Computation failed since Maxima requested additional
constraints; using the 'assume' command before evaluation
*may* help (example of legal syntax is 'assume(a>0)', see
`assume?` for more details)
Is a positive, negative or zero?

With this example, Maxima is looking for a LOT of information:

sage: assume(a>0)
sage: limit(x^a,x=0)  # random - maxima 5.46.0 does not need extra assumption
Traceback (most recent call last):
...
ValueError: Computation failed since Maxima requested additional
constraints; using the 'assume' command before evaluation *may* help
(example of legal syntax is 'assume(a>0)', see `assume?` for
 more details)
Is a an integer?
sage: assume(a,'integer')
sage: limit(x^a, x=0)  # random - maxima 5.46.0 does not need extra assumption
Traceback (most recent call last):
...
ValueError: Computation failed since Maxima requested additional
constraints; using the 'assume' command before evaluation *may* help
(example of legal syntax is 'assume(a>0)', see `assume?` for
 more details)
Is a an even number?
sage: assume(a, 'even')
sage: limit(x^a, x=0)
0
sage: forget()

More examples:

sage: limit(x*log(x), x=0, dir='+')
0
sage: lim((x+1)^(1/x), x=0)
e
sage: lim(e^x/x, x=oo)
+Infinity
sage: lim(e^x/x, x=-oo)
0
sage: lim(-e^x/x, x=oo)
-Infinity
sage: lim((cos(x))/(x^2), x=0)
+Infinity
sage: lim(sqrt(x^2+1) - x, x=oo)
0
sage: lim(x^2/(sec(x)-1), x=0)
2
sage: lim(cos(x)/(cos(x)-1), x=0)
-Infinity
sage: lim(x*sin(1/x), x=0)
0
sage: limit(e^(-1/x), x=0, dir='right')
0
sage: limit(e^(-1/x), x=0, dir='left')
+Infinity

sage: f = log(log(x)) / log(x)
sage: forget(); assume(x < -2); lim(f, x=0, taylor=True)
0
sage: forget()

Here ind means “indefinite but bounded”:

sage: lim(sin(1/x), x = 0)
ind

We can use other packages than maxima, namely “sympy”, “giac”, “fricas”.

With the standard package Giac:

sage: from sage.libs.giac.giac import libgiac     # random
sage: (exp(-x)/(2+sin(x))).limit(x=oo, algorithm='giac')
0
sage: limit(e^(-1/x), x=0, dir='right', algorithm='giac')
0
sage: limit(e^(-1/x), x=0, dir='left', algorithm='giac')
+Infinity
sage: (x / (x+2^x+cos(x))).limit(x=-infinity, algorithm='giac')
1

With the optional package FriCAS:

sage: (x / (x+2^x+cos(x))).limit(x=-infinity, algorithm='fricas')       # optional - fricas
1
sage: limit(e^(-1/x), x=0, dir='right', algorithm='fricas')             # optional - fricas
0
sage: limit(e^(-1/x), x=0, dir='left', algorithm='fricas')              # optional - fricas
+Infinity

One can also call Mathematica’s online interface:

sage: limit(pi+log(x)/x,x=oo, algorithm='mathematica_free') # optional - internet
pi

sage.calculus.calculus.mapped_opts(v)#

Used internally when creating a string of options to pass to Maxima.

INPUT:

v – an object

OUTPUT: a string.

The main use of this is to turn Python bools into lower case strings.

EXAMPLES:

sage: sage.calculus.calculus.mapped_opts(True)
'true'
sage: sage.calculus.calculus.mapped_opts(False)
'false'
sage: sage.calculus.calculus.mapped_opts('bar')
'bar'

sage.calculus.calculus.maxima_options(**kwds)#

Used internally to create a string of options to pass to Maxima.

EXAMPLES:

sage: sage.calculus.calculus.maxima_options(an_option=True, another=False, foo='bar')
'an_option=true,another=false,foo=bar'

sage.calculus.calculus.minpoly(ex, var='x', algorithm=None, bits=None, degree=None, epsilon=0)#

Return the minimal polynomial of self, if possible.

INPUT:

var – polynomial variable name (default ‘x’)
algorithm – 'algebraic' or 'numerical' (default both, but with numerical first)
bits – the number of bits to use in numerical approx
degree – the expected algebraic degree
epsilon – return without error as long as f(self) epsilon, in the case that the result cannot be proven.

All of the above parameters are optional, with epsilon=0, bits and degree tested up to 1000 and 24 by default respectively. The numerical algorithm will be faster if bits and/or degree are given explicitly. The algebraic algorithm ignores the last three parameters.

OUTPUT: The minimal polynomial of self. If the numerical algorithm is used, then it is proved symbolically when epsilon=0 (default).

If the minimal polynomial could not be found, two distinct kinds of errors are raised. If no reasonable candidate was found with the given bits/degree parameters, a ValueError will be raised. If a reasonable candidate was found but (perhaps due to limits in the underlying symbolic package) was unable to be proved correct, a NotImplementedError will be raised.

ALGORITHM: Two distinct algorithms are used, depending on the algorithm parameter. By default, the numerical algorithm is attempted first, then the algebraic one.

Algebraic: Attempt to evaluate this expression in QQbar, using cyclotomic fields to resolve exponential and trig functions at rational multiples of \(\pi\), field extensions to handle roots and rational exponents, and computing compositums to represent the full expression as an element of a number field where the minimal polynomial can be computed exactly. The bits, degree, and epsilon parameters are ignored.

Numerical: Computes a numerical approximation of self and use PARI’s pari:algdep to get a candidate minpoly \(f\). If \(f(\mathtt{self})\), evaluated to a higher precision, is close enough to 0 then evaluate \(f(\mathtt{self})\) symbolically, attempting to prove vanishing. If this fails, and epsilon is non-zero, return \(f\) if and only if \(f(\mathtt{self}) < \mathtt{epsilon}\). Otherwise raise a ValueError (if no suitable candidate was found) or a NotImplementedError (if a likely candidate was found but could not be proved correct).

EXAMPLES: First some simple examples:

sage: sqrt(2).minpoly()
x^2 - 2
sage: minpoly(2^(1/3))
x^3 - 2
sage: minpoly(sqrt(2) + sqrt(-1))
x^4 - 2*x^2 + 9
sage: minpoly(sqrt(2)-3^(1/3))
x^6 - 6*x^4 + 6*x^3 + 12*x^2 + 36*x + 1

Works with trig and exponential functions too.

sage: sin(pi/3).minpoly()
x^2 - 3/4
sage: sin(pi/7).minpoly()
x^6 - 7/4*x^4 + 7/8*x^2 - 7/64
sage: minpoly(exp(I*pi/17))
x^16 - x^15 + x^14 - x^13 + x^12 - x^11 + x^10 - x^9 + x^8
 - x^7 + x^6 - x^5 + x^4 - x^3 + x^2 - x + 1

Here we verify it gives the same result as the abstract number field.

sage: (sqrt(2) + sqrt(3) + sqrt(6)).minpoly()
x^4 - 22*x^2 - 48*x - 23
sage: K.<a,b> = NumberField([x^2-2, x^2-3])
sage: (a+b+a*b).absolute_minpoly()
x^4 - 22*x^2 - 48*x - 23

The minpoly() function is used implicitly when creating number fields:

sage: x = var('x')
sage: eqn =  x^3 + sqrt(2)*x + 5 == 0
sage: a = solve(eqn, x)[0].rhs()
sage: QQ[a]
Number Field in a with defining polynomial x^6 + 10*x^3 - 2*x^2 + 25
 with a = 0.7185272465828846? - 1.721353471724806?*I

Here we solve a cubic and then recover it from its complicated radical expansion.

sage: f = x^3 - x + 1
sage: a = f.solve(x)[0].rhs(); a
-1/2*(1/18*sqrt(23)*sqrt(3) - 1/2)^(1/3)*(I*sqrt(3) + 1)
 - 1/6*(-I*sqrt(3) + 1)/(1/18*sqrt(23)*sqrt(3) - 1/2)^(1/3)
sage: a.minpoly()
x^3 - x + 1

Note that simplification may be necessary to see that the minimal polynomial is correct.

sage: a = sqrt(2)+sqrt(3)+sqrt(5)
sage: f = a.minpoly(); f
x^8 - 40*x^6 + 352*x^4 - 960*x^2 + 576
sage: f(a)
(sqrt(5) + sqrt(3) + sqrt(2))^8 - 40*(sqrt(5) + sqrt(3) + sqrt(2))^6
 + 352*(sqrt(5) + sqrt(3) + sqrt(2))^4 - 960*(sqrt(5) + sqrt(3) + sqrt(2))^2
 + 576
sage: f(a).expand()
0

sage: a = sin(pi/7)
sage: f = a.minpoly(algorithm='numerical'); f
x^6 - 7/4*x^4 + 7/8*x^2 - 7/64
sage: f(a).horner(a).numerical_approx(100)
0.00000000000000000000000000000

The degree must be high enough (default tops out at 24).

sage: a = sqrt(3) + sqrt(2)
sage: a.minpoly(algorithm='numerical', bits=100, degree=3)
Traceback (most recent call last):
...
ValueError: Could not find minimal polynomial (100 bits, degree 3).
sage: a.minpoly(algorithm='numerical', bits=100, degree=10)
x^4 - 10*x^2 + 1

sage: cos(pi/33).minpoly(algorithm='algebraic')
x^10 + 1/2*x^9 - 5/2*x^8 - 5/4*x^7 + 17/8*x^6 + 17/16*x^5
 - 43/64*x^4 - 43/128*x^3 + 3/64*x^2 + 3/128*x + 1/1024
sage: cos(pi/33).minpoly(algorithm='numerical')
x^10 + 1/2*x^9 - 5/2*x^8 - 5/4*x^7 + 17/8*x^6 + 17/16*x^5
 - 43/64*x^4 - 43/128*x^3 + 3/64*x^2 + 3/128*x + 1/1024

Sometimes it fails, as it must given that some numbers aren’t algebraic:

sage: sin(1).minpoly(algorithm='numerical')
Traceback (most recent call last):
...
ValueError: Could not find minimal polynomial (1000 bits, degree 24).

Note

Of course, failure to produce a minimal polynomial does not necessarily indicate that this number is transcendental.

sage.calculus.calculus.mma_free_limit(expression, v, a, dir=None)#

Limit using Mathematica’s online interface.

INPUT:

expression – symbolic expression
v – variable
a – value where the variable goes to
dir – '+', '-' or None (optional, default: None)

EXAMPLES:

sage: from sage.calculus.calculus import mma_free_limit
sage: mma_free_limit(sin(x)/x, x, a=0) # optional - internet
1

Another simple limit:

sage: mma_free_limit(e^(-x), x, a=oo) # optional - internet
0

sage.calculus.calculus.nintegral(ex, x, a, b, desired_relative_error='1e-8', maximum_num_subintervals=200)#

Return a floating point machine precision numerical approximation to the integral of self from \(a\) to \(b\), computed using floating point arithmetic via maxima.

INPUT:

x – variable to integrate with respect to
a – lower endpoint of integration
b – upper endpoint of integration
desired_relative_error – (default: 1e-8) the desired relative error
maximum_num_subintervals – (default: 200) maximal number of subintervals

OUTPUT:

float: approximation to the integral
float: estimated absolute error of the approximation
the number of integrand evaluations
an error code:
- 0 – no problems were encountered
- 1 – too many subintervals were done
- 2 – excessive roundoff error
- 3 – extremely bad integrand behavior
- 4 – failed to converge
- 5 – integral is probably divergent or slowly convergent
- 6 – the input is invalid; this includes the case of desired_relative_error being too small to be achieved

ALIAS: nintegrate() is the same as nintegral()

REMARK: There is also a function numerical_integral() that implements numerical integration using the GSL C library. It is potentially much faster and applies to arbitrary user defined functions.

Also, there are limits to the precision to which Maxima can compute the integral due to limitations in quadpack. In the following example, remark that the last value of the returned tuple is 6, indicating that the input was invalid, in this case because of a too high desired precision.

sage: f = x
sage: f.nintegral(x, 0, 1, 1e-14)
(0.0, 0.0, 0, 6)

EXAMPLES:

sage: f(x) = exp(-sqrt(x))
sage: f.nintegral(x, 0, 1)
(0.5284822353142306, 4.163...e-11, 231, 0)

We can also use the numerical_integral() function, which calls the GSL C library.

sage: numerical_integral(f, 0, 1)
(0.528482232253147, 6.83928460...e-07)

Note that in exotic cases where floating point evaluation of the expression leads to the wrong value, then the output can be completely wrong:

sage: f = exp(pi*sqrt(163)) - 262537412640768744

Despite appearance, \(f\) is really very close to 0, but one gets a nonzero value since the definition of float(f) is that it makes all constants inside the expression floats, then evaluates each function and each arithmetic operation using float arithmetic:

sage: float(f)
-480.0

Computing to higher precision we see the truth:

sage: f.n(200)
-7.4992740280181431112064614366622348652078895136533593355718e-13
sage: f.n(300)
-7.49927402801814311120646143662663009137292462589621789352095066181709095575681963967103004e-13

Now numerically integrating, we see why the answer is wrong:

sage: f.nintegrate(x,0,1)
(-480.000000000000..., 5.32907051820075...e-12, 21, 0)

It is just because every floating point evaluation of \(f\) returns \(-480.0\) in floating point.

Important note: using PARI/GP one can compute numerical integrals to high precision:

sage: gp.eval('intnum(x=17,42,exp(-x^2)*log(x))')
'2.565728500561051474934096410 E-127'            # 32-bit
'2.5657285005610514829176211363206621657 E-127'  # 64-bit
sage: old_prec = gp.set_real_precision(50)
sage: gp.eval('intnum(x=17,42,exp(-x^2)*log(x))')
'2.5657285005610514829173563961304957417746108003917 E-127'
sage: gp.set_real_precision(old_prec)
57

Note that the input function above is a string in PARI syntax.

sage.calculus.calculus.nintegrate(ex, x, a, b, desired_relative_error='1e-8', maximum_num_subintervals=200)#

Return a floating point machine precision numerical approximation to the integral of self from \(a\) to \(b\), computed using floating point arithmetic via maxima.

INPUT:

x – variable to integrate with respect to
a – lower endpoint of integration
b – upper endpoint of integration
desired_relative_error – (default: 1e-8) the desired relative error
maximum_num_subintervals – (default: 200) maximal number of subintervals

OUTPUT:

float: approximation to the integral
float: estimated absolute error of the approximation
the number of integrand evaluations
an error code:
- 0 – no problems were encountered
- 1 – too many subintervals were done
- 2 – excessive roundoff error
- 3 – extremely bad integrand behavior
- 4 – failed to converge
- 5 – integral is probably divergent or slowly convergent
- 6 – the input is invalid; this includes the case of desired_relative_error being too small to be achieved

ALIAS: nintegrate() is the same as nintegral()

REMARK: There is also a function numerical_integral() that implements numerical integration using the GSL C library. It is potentially much faster and applies to arbitrary user defined functions.

sage: f = x
sage: f.nintegral(x, 0, 1, 1e-14)
(0.0, 0.0, 0, 6)

EXAMPLES:

sage: f(x) = exp(-sqrt(x))
sage: f.nintegral(x, 0, 1)
(0.5284822353142306, 4.163...e-11, 231, 0)

We can also use the numerical_integral() function, which calls the GSL C library.

sage: numerical_integral(f, 0, 1)
(0.528482232253147, 6.83928460...e-07)

Note that in exotic cases where floating point evaluation of the expression leads to the wrong value, then the output can be completely wrong:

sage: f = exp(pi*sqrt(163)) - 262537412640768744

sage: float(f)
-480.0

Computing to higher precision we see the truth:

sage: f.n(200)
-7.4992740280181431112064614366622348652078895136533593355718e-13
sage: f.n(300)
-7.49927402801814311120646143662663009137292462589621789352095066181709095575681963967103004e-13

Now numerically integrating, we see why the answer is wrong:

sage: f.nintegrate(x,0,1)
(-480.000000000000..., 5.32907051820075...e-12, 21, 0)

It is just because every floating point evaluation of \(f\) returns \(-480.0\) in floating point.

Important note: using PARI/GP one can compute numerical integrals to high precision:

sage: gp.eval('intnum(x=17,42,exp(-x^2)*log(x))')
'2.565728500561051474934096410 E-127'            # 32-bit
'2.5657285005610514829176211363206621657 E-127'  # 64-bit
sage: old_prec = gp.set_real_precision(50)
sage: gp.eval('intnum(x=17,42,exp(-x^2)*log(x))')
'2.5657285005610514829173563961304957417746108003917 E-127'
sage: gp.set_real_precision(old_prec)
57

Note that the input function above is a string in PARI syntax.

sage.calculus.calculus.symbolic_expression_from_maxima_string(x, equals_sub=False, maxima=Maxima_lib)#

Given a string representation of a Maxima expression, parse it and return the corresponding Sage symbolic expression.

INPUT:

x – a string
equals_sub – (default: False) if True, replace ‘=’ by ‘==’ in self
maxima – (default: the calculus package’s copy of Maxima) the Maxima interpreter to use.

EXAMPLES:

sage: from sage.calculus.calculus import symbolic_expression_from_maxima_string as sefms
sage: sefms('x^%e + %e^%pi + %i + sin(0)')
x^e + e^pi + I
sage: f = function('f')(x)
sage: sefms('?%at(f(x),x=2)#1')
f(2) != 1
sage: a = sage.calculus.calculus.maxima("x#0"); a
x # 0
sage: a.sage()
x != 0

sage.calculus.calculus.symbolic_expression_from_string(s, syms, accept_sequence=None, parser=False)#

Given a string, (attempt to) parse it and return the corresponding Sage symbolic expression. Normally used to return Maxima output to the user.

INPUT:

s – a string
syms – (default: {}) dictionary of strings to be regarded as symbols or functions; keys are pairs (string, number of arguments)
accept_sequence – (default: False) controls whether to allow a (possibly nested) set of lists and tuples as input
parser – (default: SR_parser) parser for internal use

EXAMPLES:

sage: from sage.calculus.calculus import symbolic_expression_from_string
sage: y = var('y')
sage: symbolic_expression_from_string('[sin(0)*x^2,3*spam+e^pi]',
....:                                 syms={('spam',0): y}, accept_sequence=True)
[0, 3*y + e^pi]

sage.calculus.calculus.symbolic_product(expression, v, a, b, algorithm='maxima', hold=False)#

Return the symbolic product \(\prod_{v = a}^b expression\) with respect to the variable \(v\) with endpoints \(a\) and \(b\).

INPUT:

expression – a symbolic expression
v – a variable or variable name
a – lower endpoint of the product
b – upper endpoint of the prduct
algorithm – (default: 'maxima') one of
- 'maxima' – use Maxima (the default)
- 'giac' – use Giac
- 'sympy' – use SymPy
- 'mathematica' – (optional) use Mathematica
hold - (default: False) if True, don’t evaluate

EXAMPLES:

sage: i, k, n = var('i,k,n')
sage: from sage.calculus.calculus import symbolic_product
sage: symbolic_product(k, k, 1, n)
factorial(n)
sage: symbolic_product(x + i*(i+1)/2, i, 1, 4)
x^4 + 20*x^3 + 127*x^2 + 288*x + 180
sage: symbolic_product(i^2, i, 1, 7)
25401600
sage: f = function('f')
sage: symbolic_product(f(i), i, 1, 7)
f(7)*f(6)*f(5)*f(4)*f(3)*f(2)*f(1)
sage: symbolic_product(f(i), i, 1, n)
product(f(i), i, 1, n)
sage: assume(k>0)
sage: symbolic_product(integrate (x^k, x, 0, 1), k, 1, n)
1/factorial(n + 1)
sage: symbolic_product(f(i), i, 1, n).log().log_expand()
sum(log(f(i)), i, 1, n)

sage.calculus.calculus.symbolic_sum(expression, v, a, b, algorithm='maxima', hold=False)#

Return the symbolic sum \(\sum_{v = a}^b expression\) with respect to the variable \(v\) with endpoints \(a\) and \(b\).

INPUT:

expression – a symbolic expression
v – a variable or variable name
a – lower endpoint of the sum
b – upper endpoint of the sum
algorithm – (default: 'maxima') one of
- 'maxima' – use Maxima (the default)
- 'maple' – (optional) use Maple
- 'mathematica' – (optional) use Mathematica
- 'giac' – (optional) use Giac
- 'sympy' – use SymPy
hold – (default: False) if True, don’t evaluate

EXAMPLES:

sage: k, n = var('k,n')
sage: from sage.calculus.calculus import symbolic_sum
sage: symbolic_sum(k, k, 1, n).factor()
1/2*(n + 1)*n

sage: symbolic_sum(1/k^4, k, 1, oo)
1/90*pi^4

sage: symbolic_sum(1/k^5, k, 1, oo)
zeta(5)

A well known binomial identity:

sage: symbolic_sum(binomial(n,k), k, 0, n)
2^n

And some truncations thereof:

sage: assume(n>1)
sage: symbolic_sum(binomial(n,k), k, 1, n)
2^n - 1
sage: symbolic_sum(binomial(n,k), k, 2, n)
2^n - n - 1
sage: symbolic_sum(binomial(n,k), k, 0, n-1)
2^n - 1
sage: symbolic_sum(binomial(n,k), k, 1, n-1)
2^n - 2

The binomial theorem:

sage: x, y = var('x, y')
sage: symbolic_sum(binomial(n,k) * x^k * y^(n-k), k, 0, n)
(x + y)^n

sage: symbolic_sum(k * binomial(n, k), k, 1, n)
2^(n - 1)*n

sage: symbolic_sum((-1)^k*binomial(n,k), k, 0, n)
0

sage: symbolic_sum(2^(-k)/(k*(k+1)), k, 1, oo)
-log(2) + 1

Summing a hypergeometric term:

sage: symbolic_sum(binomial(n, k) * factorial(k) / factorial(n+1+k), k, 0, n)
1/2*sqrt(pi)/factorial(n + 1/2)

We check a well known identity:

sage: bool(symbolic_sum(k^3, k, 1, n) == symbolic_sum(k, k, 1, n)^2)
True

A geometric sum:

sage: a, q = var('a, q')
sage: symbolic_sum(a*q^k, k, 0, n)
(a*q^(n + 1) - a)/(q - 1)

For the geometric series, we will have to assume the right values for the sum to converge:

sage: assume(abs(q) < 1)
sage: symbolic_sum(a*q^k, k, 0, oo)
-a/(q - 1)

A divergent geometric series. Don’t forget to forget your assumptions:

sage: forget()
sage: assume(q > 1)
sage: symbolic_sum(a*q^k, k, 0, oo)
Traceback (most recent call last):
...
ValueError: Sum is divergent.
sage: forget()
sage: assumptions()  # check the assumptions were really forgotten
[]

A summation performed by Mathematica:

sage: symbolic_sum(1/(1+k^2), k, -oo, oo, algorithm='mathematica')     # optional - mathematica
pi*coth(pi)

An example of this summation with Giac:

sage: symbolic_sum(1/(1+k^2), k, -oo, oo, algorithm='giac')
(pi*e^(2*pi) - pi*e^(-2*pi))/(e^(2*pi) + e^(-2*pi) - 2)

The same summation is solved by SymPy:

sage: symbolic_sum(1/(1+k^2), k, -oo, oo, algorithm='sympy')
pi/tanh(pi)

SymPy and Maxima 5.39.0 can do the following (see github issue #22005):

sage: sum(1/((2*n+1)^2-4)^2, n, 0, Infinity, algorithm='sympy')
1/64*pi^2
sage: sum(1/((2*n+1)^2-4)^2, n, 0, Infinity)
1/64*pi^2

Use Maple as a backend for summation:

sage: symbolic_sum(binomial(n,k)*x^k, k, 0, n, algorithm='maple')      # optional - maple
(x + 1)^n

If you don’t want to evaluate immediately give the hold keyword:

sage: s = sum(n, n, 1, k, hold=True); s
sum(n, n, 1, k)
sage: s.unhold()
1/2*k^2 + 1/2*k
sage: s.subs(k == 10)
sum(n, n, 1, 10)
sage: s.subs(k == 10).unhold()
55
sage: s.subs(k == 10).n()
55.0000000000000

Note

Sage can currently only understand a subset of the output of Maxima, Maple and Mathematica, so even if the chosen backend can perform the summation the result might not be convertible into a Sage expression.