Symbolic Computation — Sage 9.4 Reference Manual: Symbolic Calculus (2024)

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 trac ticket #9217

EXAMPLES:

The basic units of the calculus package are symbolic expressions whichare elements of the symbolic expression ring (SR). To create asymbolic variable object in Sage, use the var() function, whoseargument is the text of that variable. Note that Sage is intelligentabout LaTeXing variable names.

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

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

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

More complicated expressions in Sage can be built up using ordinaryarithmetic. The following are valid, and follow the rules of Pythonarithmetic: (The ‘=’ operator represents assignment, and notequality)

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 thescenes through Maxima:

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

Note that these methods usually require an explicit variable name. If noneis 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., thevariable order is specified), then Sage can calculate the matrixderivative (i.e., the gradient, Jacobian matrix, etc.) if no variablesare specified. In the example below, we use the second derivativetest to determine that there is a saddle point at (0,-1/2).

sage: f(x,y)=x^2*y+y^2+ysage: 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 anexception 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 thesum 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^4sage: 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))]

And complex exponentiation works now:

sage: M = i*matrix([[pi]])sage: e^M[-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 pythondict:

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)

It was formerly the case that if there was no ambiguity of variablenames, we didn’t have to specify them; that still works for the moment,but the behavior is deprecated:

sage: f = sin(x)sage: f(y)doctest:...: DeprecationWarning: Substitution using function-callsyntax and unnamed arguments is deprecated and will be removedfrom a future release of Sage; you can use named arguments instead,like EXPR(x=..., y=...)See http://trac.sagemath.org/5930 for details.sin(y)sage: f(pi)0

However if there is ambiguity, we should explicitly state whatvariables we’re substituting for:

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

We can also make a CallableSymbolicExpression,which is a SymbolicExpression that is a function ofspecified variables in a fixed order. EachSymbolicExpression has afunction(...) method that is used to create aCallableSymbolicExpression, 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 aCallableSymbolicExpression, which relies on theSage 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 ambiguitysubstituting 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: f5*sin(x)sage: f(x=2)5*sin(2)sage: f(x=pi)0sage: float(f(x=pi))0.0

Another example:

sage: f = integrate(1/sqrt(9+x^2), x); farcsinh(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^2sage: dy = derivative(y,x)sage: z = integral(sqrt(1 + dy^2), x, 0, 2)sage: zsqrt(17) + 1/4*arcsinh(4)sage: n(z,200)4.6467837624329358733826155674904591885104869874232887508703sage: 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*Isage: CC(2*I)2.00000000000000*Isage: ComplexField(200)(2*I)2.0000000000000000000000000000000000000000000000000000000000*Isage: ComplexField(200)(sin(I))1.1752011936438014568823818505956008151557179813340958702296*Isage: f = sin(I) + cos(I/2); fcosh(1/2) + I*sinh(1)sage: CC(f)1.12762596520638 + 1.17520119364380*Isage: ComplexField(200)(f)1.1276259652063807852262251614026720125478471180986674836290 + 1.1752011936438014568823818505956008151557179813340958702296*Isage: ComplexField(100)(f)1.1276259652063807852262251614 + 1.1752011936438014568823818506*I

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

sage: f = sum(1/var('n%s'%i)^i for i in range(10))sage: f1/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 immediatelyavailable 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 ofa symbolic expression.

EXAMPLES:

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

sage: g = x^3-3sage: from sage.calculus.calculus import atsage: at(g, x=1)-2sage: 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 (trac ticket #12796):

sage: from sage.calculus.calculus import atsage: 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_diffsage: 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); gdiff(cr(a), a)

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

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

EXAMPLES:

sage: from sage.calculus.calculus import dummy_integratesage: 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 laplacetransforms that Maxima can’t compute:

EXAMPLES:

sage: from sage.calculus.calculus import dummy_inverse_laplacesage: 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 transformsthat Maxima can’t compute:

EXAMPLES:

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

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

Return the inverse Laplace transform with respect to the variable \(t\) andtransform 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(s) = \frac{1}{2\pi i} \int_{\gamma-i\infty}^{\gamma + i\infty} e^{st} F(s) dt,\]

where \(\gamma\) is chosen so that the contour path ofintegration 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); gdiff(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 DavidJoyner): 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 have90, but the Y-men are better at fighting, because of the higherfactor of “-16” vs “-1”, and also get an occasional reinforcement,because of the “+1” term.

sage: var('t')tsage: t = var('t')sage: x = function('x')(t)sage: y = function('y')(t)sage: de1 = x.diff(t) + 16*ysage: de2 = y.diff(t) + x - 1sage: 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) + 1sage: yt91/8*e^(4*t) + 629/8*e^(-4*t)sage: p1 = plot(xt,0,1/2,rgbcolor=(1,0,0))sage: p2 = plot(yt,0,1/2,rgbcolor=(0,1,0))sage: (p1+p2).save(os.path.join(SAGE_TMP, "de_plot.png"))

Another example:

sage: var('a,s,t')(a, s, t)sage: f = exp (2*t + a) * sin(t) * t; ft*e^(a + 2*t)*sin(t)sage: L = laplace(f, t, s); L2*(s - 2)*e^a/(s^2 - 4*s + 5)^2sage: 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 isrelated to a choice that has to be made when defining Laplacetransforms of distributions):

sage: laplace(dirac_delta(t), t, s)1sage: laplace(dirac_delta(t), t, s, algorithm='sympy')(-heaviside(0) + 1, -oo, True)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); dir 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 Taylorseries, which allows more limits to be computed (but may alsocrash 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’ (complexinfinity).

EXAMPLES:

sage: x = var('x')sage: f = (1 + 1/x)^xsage: f.limit(x=oo)esage: f.limit(x=5)7776/3125sage: f.limit(x=1.2)2.06961575467...sage: f.limit(x=I, taylor=True)(-I + 1)^Isage: f(x=1.2)2.0696157546720...sage: f(x=I)(-I + 1)^Isage: CDF(f(x=I))2.0628722350809046 + 0.7450070621797239*Isage: CDF(f.limit(x=I))2.0628722350809046 + 0.7450070621797239*I

Notice that Maxima may ask for more information:

sage: var('a')asage: limit(x^a,x=0)Traceback (most recent call last):...ValueError: Computation failed since Maxima requested additionalconstraints; 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)Traceback (most recent call last):...ValueError: Computation failed since Maxima requested additionalconstraints; 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)Traceback (most recent call last):...ValueError: Computation failed since Maxima requested additionalconstraints; 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)0sage: forget()

More examples:

sage: limit(x*log(x), x=0, dir='+')0sage: lim((x+1)^(1/x), x=0)esage: lim(e^x/x, x=oo)+Infinitysage: lim(e^x/x, x=-oo)0sage: lim(-e^x/x, x=oo)-Infinitysage: lim((cos(x))/(x^2), x=0)+Infinitysage: lim(sqrt(x^2+1) - x, x=oo)0sage: lim(x^2/(sec(x)-1), x=0)2sage: lim(cos(x)/(cos(x)-1), x=0)-Infinitysage: lim(x*sin(1/x), x=0)0sage: limit(e^(-1/x), x=0, dir='right')0sage: 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)0sage: 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 # randomsage: (exp(-x)/(2+sin(x))).limit(x=oo, algorithm='giac')0sage: limit(e^(-1/x), x=0, dir='right', algorithm='giac')0sage: limit(e^(-1/x), x=0, dir='left', algorithm='giac')+Infinitysage: (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 - fricas1sage: limit(e^(-1/x), x=0, dir='right', algorithm='fricas') # optional - fricas0sage: 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 - internetpi

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); dir 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 Taylorseries, which allows more limits to be computed (but may alsocrash 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’ (complexinfinity).

EXAMPLES:

sage: x = var('x')sage: f = (1 + 1/x)^xsage: f.limit(x=oo)esage: f.limit(x=5)7776/3125sage: f.limit(x=1.2)2.06961575467...sage: f.limit(x=I, taylor=True)(-I + 1)^Isage: f(x=1.2)2.0696157546720...sage: f(x=I)(-I + 1)^Isage: CDF(f(x=I))2.0628722350809046 + 0.7450070621797239*Isage: CDF(f.limit(x=I))2.0628722350809046 + 0.7450070621797239*I

Notice that Maxima may ask for more information:

sage: var('a')asage: limit(x^a,x=0)Traceback (most recent call last):...ValueError: Computation failed since Maxima requested additionalconstraints; 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)Traceback (most recent call last):...ValueError: Computation failed since Maxima requested additionalconstraints; 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)Traceback (most recent call last):...ValueError: Computation failed since Maxima requested additionalconstraints; 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)0sage: forget()

More examples:

sage: limit(x*log(x), x=0, dir='+')0sage: lim((x+1)^(1/x), x=0)esage: lim(e^x/x, x=oo)+Infinitysage: lim(e^x/x, x=-oo)0sage: lim(-e^x/x, x=oo)-Infinitysage: lim((cos(x))/(x^2), x=0)+Infinitysage: lim(sqrt(x^2+1) - x, x=oo)0sage: lim(x^2/(sec(x)-1), x=0)2sage: lim(cos(x)/(cos(x)-1), x=0)-Infinitysage: lim(x*sin(1/x), x=0)0sage: limit(e^(-1/x), x=0, dir='right')0sage: 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)0sage: 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 # randomsage: (exp(-x)/(2+sin(x))).limit(x=oo, algorithm='giac')0sage: limit(e^(-1/x), x=0, dir='right', algorithm='giac')0sage: limit(e^(-1/x), x=0, dir='left', algorithm='giac')+Infinitysage: (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 - fricas1sage: limit(e^(-1/x), x=0, dir='right', algorithm='fricas') # optional - fricas0sage: 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 - internetpi

sage.calculus.calculus.mapped_opts(v)¶

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

INPUT:

v - an object

OUTPUT: a string.

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

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’ (defaultboth, but with numerical first)
bits - the number of bits to use in numericalapprox
degree - the expected algebraic degree
epsilon - return without error as long asf(self) epsilon, in the case that the result cannot be proven.
All of the above parameters are optional, with epsilon=0, bits anddegree tested up to 1000 and 24 by default respectively. Thenumerical algorithm will be faster if bits and/or degree are givenexplicitly. The algebraic algorithm ignores the last threeparameters.

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

If the minimal polynomial could not be found, two distinct kinds oferrors are raised. If no reasonable candidate was found with thegiven bit/degree parameters, a ValueError will beraised. If a reasonable candidate was found but (perhaps due tolimits in the underlying symbolic package) was unable to be provedcorrect, a NotImplementedError will be raised.

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

Algebraic: Attempt to evaluate this expression in QQbar, usingcyclotomic fields to resolve exponential and trig functions atrational multiples of pi, field extensions to handle roots andrational exponents, and computing compositums to represent the fullexpression as an element of a number field where the minimalpolynomial can be computed exactly. The bits, degree, and epsilonparameters are ignored.

Numerical: Computes a numerical approximation ofself and use PARI’s algdep to get a candidateminpoly \(f\). If \(f(\mathtt{self})\),evaluated to a higher precision, is close enough to 0 then evaluate\(f(\mathtt{self})\) symbolically, attempting to provevanishing. 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 suitablecandidate was found) or a NotImplementedError (if alikely candidate was found but could not be proved correct).

EXAMPLES: First some simple examples:

sage: sqrt(2).minpoly()x^2 - 2sage: minpoly(2^(1/3))x^3 - 2sage: minpoly(sqrt(2) + sqrt(-1))x^4 - 2*x^2 + 9sage: 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/4sage: sin(pi/7).minpoly()x^6 - 7/4*x^4 + 7/8*x^2 - 7/64sage: 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 numberfield.

sage: (sqrt(2) + sqrt(3) + sqrt(6)).minpoly()x^4 - 22*x^2 - 48*x - 23sage: 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 creatingnumber fields:

sage: x = var('x')sage: eqn = x^3 + sqrt(2)*x + 5 == 0sage: 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 complicatedradical expansion.

sage: f = x^3 - x + 1sage: 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 minimalpolynomial is correct.

sage: a = sqrt(2)+sqrt(3)+sqrt(5)sage: f = a.minpoly(); fx^8 - 40*x^6 + 352*x^4 - 960*x^2 + 576sage: 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 + 576sage: f(a).expand()0

sage: a = sin(pi/7)sage: f = a.minpoly(algorithm='numerical'); fx^6 - 7/4*x^4 + 7/8*x^2 - 7/64sage: 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/1024sage: 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 notnecessarily indicate that this number is transcendental.

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

Limit using Mathematica’s online interface.

EXAMPLES:

sage: from sage.calculus.calculus import mma_free_limitsage: mma_free_limit(sin(x)/x, x=0) # optional - internet1

Another simple limit:

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

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

Return a floating point machine precision numerical approximationto 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’) thedesired relative error
maximum_num_subintervals - (default: 200)maxima number of subintervals

OUTPUT:

float: approximation to the integral
float: estimated absolute error of theapproximation
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 slowlyconvergent
- 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 functionnumerical_integral that implements numericalintegration using the GSL C library. It is potentially much fasterand applies to arbitrary user defined functions.

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

sage: f = xsage: 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 theexpression leads to the wrong value, then the output can becompletely wrong:

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

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

sage: float(f)-480.0

Computing to higher precision we see the truth:

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

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

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

It is just because every floating point evaluation of return -480.0in floating point.

Important note: using PARI/GP one can compute numerical integralsto 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-bitsage: 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 approximationto 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’) thedesired relative error
maximum_num_subintervals - (default: 200)maxima number of subintervals

OUTPUT:

float: approximation to the integral
float: estimated absolute error of theapproximation
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 slowlyconvergent
- 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 functionnumerical_integral that implements numericalintegration using the GSL C library. It is potentially much fasterand applies to arbitrary user defined functions.

sage: f = xsage: 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 theexpression leads to the wrong value, then the output can becompletely 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-13sage: f.n(300)-7.49927402801814311120646143662663009137292462589621789352095066181709095575681963967103004e-13

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

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

It is just because every floating point evaluation of return -480.0in floating point.

Important note: using PARI/GP one can compute numerical integralsto 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-bitsage: 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 andreturn the corresponding Sage symbolic expression.

INPUT:

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

EXAMPLES:

sage: from sage.calculus.calculus import symbolic_expression_from_maxima_string as sefmssage: sefms('x^%e + %e^%pi + %i + sin(0)')x^e + e^pi + Isage: f = function('f')(x)sage: sefms('?%at(f(x),x=2)#1')f(2) != 1sage: a = sage.calculus.calculus.maxima("x#0"); ax#0sage: 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 thecorresponding Sage symbolic expression. Normally usedto return Maxima output to the user.

INPUT:

s - a string
syms - (default: {}) dictionary ofstrings to be regarded as symbols or functions
accept_sequence - (default: False) controls whetherto allow a (possibly nested) set of lists and tuplesas input
parser – (default: SR_parser) parser for internal use

EXAMPLES:

sage: y = var('y')sage: sage.calculus.calculus.symbolic_expression_from_string('[sin(0)*x^2,3*spam+e^pi]',syms={'spam':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 respectto 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_productsage: 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 + 180sage: symbolic_product(i^2, i, 1, 7)25401600sage: 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 respectto 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_sumsage: 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 - 1sage: symbolic_sum(binomial(n,k),k,2,n)2^n - n - 1sage: symbolic_sum(binomial(n,k),k,0,n-1)2^n - 1sage: 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 assumethe 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 forgetto 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 - mathematicapi*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 (seetrac ticket #22005):

sage: sum(1/((2*n+1)^2-4)^2, n, 0, Infinity, algorithm='sympy')1/64*pi^2sage: 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); ssum(n, n, 1, k)sage: s.unhold()1/2*k^2 + 1/2*ksage: s.subs(k == 10)sum(n, n, 1, 10)sage: s.subs(k == 10).unhold()55sage: 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 performthe summation the result might not be convertible into a Sageexpression.