diff --git a/src/sage/calculus/calculus.py b/src/sage/calculus/calculus.py index 6ffe90f7dda..34abf305f70 100644 --- a/src/sage/calculus/calculus.py +++ b/src/sage/calculus/calculus.py @@ -80,23 +80,23 @@ :: - sage: f(x,y)=x^2*y+y^2+y - sage: f.diff() # gradient + 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]) + 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) + sage: H(x=0, y=-1/2) [-1 0] [ 0 2] - sage: H(x=0,y=-1/2).eigenvalues() + 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 (r, theta) |--> (r*cos(theta), r*sin(theta)) sage: T.diff() # Jacobian matrix @@ -117,7 +117,7 @@ ValueError: No differentiation variable specified. Simplifying symbolic sums is also possible, using the -sum command, which also uses Maxima in the background:: +:func:`sum` command, which also uses Maxima in the background:: sage: k, m = var('k, m') sage: sum(1/k^4, k, 1, oo) @@ -182,19 +182,19 @@ sage: f(x=pi) 0 -We can also make a ``CallableSymbolicExpression``, -which is a ``SymbolicExpression`` that is a function of +We can also make a :class:`CallableSymbolicExpression`, +which is a :class:`SymbolicExpression` that is a function of specified variables in a fixed order. Each -``SymbolicExpression`` has a +:class:`SymbolicExpression` has a ``function(...)`` method that is used to create a -``CallableSymbolicExpression``, as illustrated below:: +:class:`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 +:class:`CallableSymbolicExpression`, which relies on the Sage preparser. :: @@ -267,7 +267,8 @@ sage: CC(f) 1.12762596520638 + 1.17520119364380*I sage: ComplexField(200)(f) - 1.1276259652063807852262251614026720125478471180986674836290 + 1.1752011936438014568823818505956008151557179813340958702296*I + 1.1276259652063807852262251614026720125478471180986674836290 + + 1.1752011936438014568823818505956008151557179813340958702296*I sage: ComplexField(100)(f) 1.1276259652063807852262251614 + 1.1752011936438014568823818506*I @@ -286,8 +287,9 @@ 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: 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 TESTS: @@ -317,7 +319,7 @@ sage: maxima.eval('expand((x+y)^3)') '27' -If the copy of maxima used by the symbolic calculus package were +If the copy of Maxima used by the symbolic calculus package were the same as the default one, then the following would return 27, which would be very confusing indeed! @@ -446,27 +448,27 @@ def symbolic_sum(expression, v, a, b, algorithm='maxima', hold=False): INPUT: - - ``expression`` - a symbolic expression + - ``expression`` -- a symbolic expression - - ``v`` - a variable or variable name + - ``v`` -- a variable or variable name - - ``a`` - lower endpoint of the sum + - ``a`` -- lower endpoint of the sum - - ``b`` - upper endpoint of the sum + - ``b`` -- upper endpoint of the sum - - ``algorithm`` - (default: ``'maxima'``) one of + - ``algorithm`` -- (default: ``'maxima'``) one of - - ``'maxima'`` - use Maxima (the default) + - ``'maxima'`` -- use Maxima (the default) - - ``'maple'`` - (optional) use Maple + - ``'maple'`` -- (optional) use Maple - - ``'mathematica'`` - (optional) use Mathematica + - ``'mathematica'`` -- (optional) use Mathematica - - ``'giac'`` - (optional) use Giac + - ``'giac'`` -- (optional) use Giac - - ``'sympy'`` - use SymPy + - ``'sympy'`` -- use SymPy - - ``hold`` - (default: ``False``) if ``True`` don't evaluate + - ``hold`` -- (default: ``False``) if ``True``, don't evaluate EXAMPLES:: @@ -493,13 +495,13 @@ def symbolic_sum(expression, v, a, b, algorithm='maxima', hold=False): And some truncations thereof:: sage: assume(n>1) - sage: symbolic_sum(binomial(n,k),k,1,n) + sage: symbolic_sum(binomial(n,k), k, 1, n) 2^n - 1 - sage: symbolic_sum(binomial(n,k),k,2,n) + sage: symbolic_sum(binomial(n,k), k, 2, n) 2^n - n - 1 - sage: symbolic_sum(binomial(n,k),k,0,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) + sage: symbolic_sum(binomial(n,k), k, 1, n-1) 2^n - 2 The binomial theorem:: @@ -556,22 +558,22 @@ def symbolic_sum(expression, v, a, b, algorithm='maxima', hold=False): ... ValueError: Sum is divergent. sage: forget() - sage: assumptions() # check the assumptions were really forgotten + 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 + 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') + 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') + 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 @@ -584,7 +586,7 @@ def symbolic_sum(expression, v, a, b, algorithm='maxima', hold=False): Use Maple as a backend for summation:: - sage: symbolic_sum(binomial(n,k)*x^k, k, 0, n, algorithm = 'maple') # optional - maple + 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:: @@ -686,17 +688,17 @@ def nintegral(ex, x, a, b, INPUT: - - ``x`` - variable to integrate with respect to + - ``x`` -- variable to integrate with respect to - - ``a`` - lower endpoint of integration + - ``a`` -- lower endpoint of integration - - ``b`` - upper endpoint of integration + - ``b`` -- upper endpoint of integration - - ``desired_relative_error`` - (default: '1e-8') the + - ``desired_relative_error`` -- (default: ``1e-8``) the desired relative error - - ``maximum_num_subintervals`` - (default: 200) - maxima number of subintervals + - ``maximum_num_subintervals`` -- (default: 200) + maximal number of subintervals OUTPUT: @@ -709,26 +711,26 @@ def nintegral(ex, x, a, b, - an error code: - - ``0`` - no problems were encountered + - ``0`` -- no problems were encountered - - ``1`` - too many subintervals were done + - ``1`` -- too many subintervals were done - - ``2`` - excessive roundoff error + - ``2`` -- excessive roundoff error - - ``3`` - extremely bad integrand behavior + - ``3`` -- extremely bad integrand behavior - - ``4`` - failed to converge + - ``4`` -- failed to converge - - ``5`` - integral is probably divergent or slowly + - ``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 + - ``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 + ALIAS: :func:`nintegrate` is the same as :func:`nintegral` REMARK: There is also a function - ``numerical_integral`` that implements numerical + :func:`numerical_integral` that implements numerical integration using the GSL C library. It is potentially much faster and applies to arbitrary user defined functions. @@ -741,7 +743,7 @@ def nintegral(ex, x, a, b, :: sage: f = x - sage: f.nintegral(x,0,1,1e-14) + sage: f.nintegral(x, 0, 1, 1e-14) (0.0, 0.0, 0, 6) EXAMPLES:: @@ -750,7 +752,7 @@ def nintegral(ex, x, a, b, sage: f.nintegral(x, 0, 1) (0.5284822353142306, 4.163...e-11, 231, 0) - We can also use the ``numerical_integral`` function, + We can also use the :func:`numerical_integral` function, which calls the GSL C library. :: @@ -785,7 +787,7 @@ def nintegral(ex, x, a, b, sage: f.nintegrate(x,0,1) (-480.000000000000..., 5.32907051820075...e-12, 21, 0) - It is just because every floating point evaluation of return -480.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 @@ -831,25 +833,25 @@ def symbolic_product(expression, v, a, b, algorithm='maxima', hold=False): INPUT: - - ``expression`` - a symbolic expression + - ``expression`` -- a symbolic expression - - ``v`` - a variable or variable name + - ``v`` -- a variable or variable name - - ``a`` - lower endpoint of the product + - ``a`` -- lower endpoint of the product - - ``b`` - upper endpoint of the prduct + - ``b`` -- upper endpoint of the prduct - - ``algorithm`` - (default: ``'maxima'``) one of + - ``algorithm`` -- (default: ``'maxima'``) one of - - ``'maxima'`` - use Maxima (the default) + - ``'maxima'`` -- use Maxima (the default) - - ``'giac'`` - use Giac + - ``'giac'`` -- use Giac - - ``'sympy'`` - use SymPy + - ``'sympy'`` -- use SymPy - - ``'mathematica'`` - (optional) use Mathematica + - ``'mathematica'`` -- (optional) use Mathematica - - ``hold`` - (default: ``False``) if ``True`` don't evaluate + - ``hold`` - (default: ``False``) if ``True``, don't evaluate EXAMPLES:: @@ -935,50 +937,50 @@ def symbolic_product(expression, v, a, b, algorithm='maxima', hold=False): def minpoly(ex, var='x', algorithm=None, bits=None, degree=None, epsilon=0): r""" - Return the minimal polynomial of self, if possible. + Return the minimal polynomial of ``self``, if possible. INPUT: - - ``var`` - polynomial variable name (default 'x') + - ``var`` -- polynomial variable name (default 'x') - - ``algorithm`` - 'algebraic' or 'numerical' (default + - ``algorithm`` -- ``'algebraic'`` or ``'numerical'`` (default both, but with numerical first) - - ``bits`` - the number of bits to use in numerical + - ``bits`` -- the number of bits to use in numerical approx - - ``degree`` - the expected algebraic degree + - ``degree`` -- the expected algebraic degree - - ``epsilon`` - return without error as long as + - ``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 + 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). + 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 bit/degree parameters, a ``ValueError`` will be + given ``bits``/``degree`` parameters, a :class:`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. + correct, a :class:`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 + 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 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 + polynomial can be computed exactly. The ``bits``, ``degree``, and ``epsilon`` parameters are ignored. Numerical: Computes a numerical approximation of @@ -989,8 +991,8 @@ def minpoly(ex, var='x', algorithm=None, bits=None, degree=None, epsilon=0): 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 + Otherwise raise a :class:`ValueError` (if no suitable + candidate was found) or a :class:`NotImplementedError` (if a likely candidate was found but could not be proved correct). EXAMPLES: First some simple examples:: @@ -1014,7 +1016,8 @@ def minpoly(ex, var='x', algorithm=None, bits=None, degree=None, epsilon=0): 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 + 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. @@ -1027,14 +1030,15 @@ def minpoly(ex, var='x', algorithm=None, bits=None, degree=None, epsilon=0): sage: (a+b+a*b).absolute_minpoly() x^4 - 22*x^2 - 48*x - 23 - The minpoly function is used implicitly when creating + The :func:`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 + 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. @@ -1043,7 +1047,8 @@ def minpoly(ex, var='x', algorithm=None, bits=None, degree=None, epsilon=0): 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) + -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 @@ -1056,7 +1061,9 @@ def minpoly(ex, var='x', algorithm=None, bits=None, degree=None, epsilon=0): 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 + (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 @@ -1083,9 +1090,11 @@ def minpoly(ex, var='x', algorithm=None, bits=None, degree=None, epsilon=0): :: 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 + 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 + 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:: @@ -1162,12 +1171,12 @@ def limit(ex, dir=None, taylor=False, algorithm='maxima', **argv): 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 + - ``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 + - ``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). @@ -1175,8 +1184,8 @@ def limit(ex, dir=None, taylor=False, algorithm='maxima', **argv): .. note:: - The output may also use 'und' (undefined), 'ind' - (indefinite but bounded), and 'infinity' (complex + The output may also use ``und`` (undefined), ``ind`` + (indefinite but bounded), and ``infinity`` (complex infinity). EXAMPLES:: @@ -1488,7 +1497,7 @@ def mma_free_limit(expression, v, a, dir=None): - ``expression`` -- symbolic expression - ``v`` -- variable - ``a`` -- value where the variable goes to - - ``dir`` -- ``'+'``, ``'-'`` or ``None`` (optional, default:``None``) + - ``dir`` -- ``'+'``, ``'-'`` or ``None`` (optional, default: ``None``) EXAMPLES:: @@ -1547,19 +1556,19 @@ def laplace(ex, t, s, algorithm='maxima'): INPUT: - - ``ex`` - a symbolic expression + - ``ex`` -- a symbolic expression - - ``t`` - independent variable + - ``t`` -- independent variable - - ``s`` - transform parameter + - ``s`` -- transform parameter - - ``algorithm`` - (default: ``'maxima'``) one of + - ``algorithm`` -- (default: ``'maxima'``) one of - - ``'maxima'`` - use Maxima (the default) + - ``'maxima'`` -- use Maxima (the default) - - ``'sympy'`` - use SymPy + - ``'sympy'`` -- use SymPy - - ``'giac'`` - use Giac + - ``'giac'`` -- use Giac .. NOTE:: @@ -1622,7 +1631,7 @@ def laplace(ex, t, s, algorithm='maxima'): Next we form the augmented matrix of the above system:: - sage: A = matrix([[s, 16, 270],[1, s, 90+1/s]]) + 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) @@ -1630,11 +1639,11 @@ def laplace(ex, t, s, algorithm='maxima'): -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)) - sage: p2 = plot(yt,0,1/2,rgbcolor=(0,1,0)) + 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: - ....: (p1+p2).save(f.name) + sage: with tempfile.NamedTemporaryFile(suffix=".png") as f: # needs sage.plot + ....: (p1 + p2).save(f.name) Another example:: @@ -1789,19 +1798,19 @@ def inverse_laplace(ex, s, t, algorithm='maxima'): INPUT: - - ``ex`` - a symbolic expression + - ``ex`` -- a symbolic expression - - ``s`` - transform parameter + - ``s`` -- transform parameter - - ``t`` - independent variable + - ``t`` -- independent variable - - ``algorithm`` - (default: ``'maxima'``) one of + - ``algorithm`` -- (default: ``'maxima'``) one of - - ``'maxima'`` - use Maxima (the default) + - ``'maxima'`` -- use Maxima (the default) - - ``'sympy'`` - use SymPy + - ``'sympy'`` -- use SymPy - - ``'giac'`` - use Giac + - ``'giac'`` -- use Giac .. SEEALSO:: @@ -1843,9 +1852,11 @@ def inverse_laplace(ex, s, t, algorithm='maxima'): Transform a rational expression:: - sage: inverse_laplace((2*s^2*exp(-2*s) - exp(-s))/(s^3+1), s, t, algorithm='giac') - -1/3*(sqrt(3)*e^(1/2*t - 1/2)*sin(1/2*sqrt(3)*(t - 1)) - cos(1/2*sqrt(3)*(t - 1))*e^(1/2*t - 1/2) + - e^(-t + 1))*heaviside(t - 1) + 2/3*(2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*heaviside(t - 2) + sage: inverse_laplace((2*s^2*exp(-2*s) - exp(-s))/(s^3+1), s, t, + ....: algorithm='giac') + -1/3*(sqrt(3)*e^(1/2*t - 1/2)*sin(1/2*sqrt(3)*(t - 1)) + - cos(1/2*sqrt(3)*(t - 1))*e^(1/2*t - 1/2) + e^(-t + 1))*heaviside(t - 1) + + 2/3*(2*cos(1/2*sqrt(3)*(t - 2))*e^(1/2*t - 1) + e^(-t + 2))*heaviside(t - 2) sage: inverse_laplace(1/(s - 1), s, x) e^x @@ -1955,7 +1966,7 @@ def at(ex, *args, **kwds): 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: g = x^3 - 3 sage: from sage.calculus.calculus import at sage: at(g, x=1) -2 @@ -1970,15 +1981,16 @@ def at(ex, *args, **kwds): 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) + 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) + 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) @@ -2078,7 +2090,7 @@ def dummy_laplace(*args): 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) + sage: dummy_laplace(f(t), t, s) laplace(f(t), t, s) """ return _laplace(args[0], var(repr(args[1])), var(repr(args[2]))) @@ -2086,7 +2098,7 @@ def dummy_laplace(*args): def dummy_inverse_laplace(*args): """ - This function is called to create formal wrappers of inverse laplace + This function is called to create formal wrappers of inverse Laplace transforms that Maxima can't compute: EXAMPLES:: @@ -2094,7 +2106,7 @@ def dummy_inverse_laplace(*args): 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) + sage: dummy_inverse_laplace(F(s), s, t) ilt(F(s), s, t) """ return _inverse_laplace(args[0], var(repr(args[1])), var(repr(args[2]))) @@ -2108,7 +2120,7 @@ def dummy_pochhammer(*args): sage: from sage.calculus.calculus import dummy_pochhammer sage: s,t = var('s,t') - sage: dummy_pochhammer(s,t) + sage: dummy_pochhammer(s, t) gamma(s + t)/gamma(s) """ x, y = args @@ -2224,12 +2236,12 @@ def symbolic_expression_from_maxima_string(x, equals_sub=False, maxima=maxima): INPUT: - - ``x`` - a string + - ``x`` -- a string - - ``equals_sub`` - (default: False) if True, replace + - ``equals_sub`` -- (default: ``False``) if ``True``, replace '=' by '==' in self - - ``maxima`` - (default: the calculus package's + - ``maxima`` -- (default: the calculus package's copy of Maxima) the Maxima interpreter to use. EXAMPLES:: @@ -2415,7 +2427,7 @@ def mapped_opts(v): INPUT: - - ``v`` - an object + - ``v`` -- an object OUTPUT: a string. @@ -2548,13 +2560,13 @@ def symbolic_expression_from_string(s, syms=None, accept_sequence=False, *, pars INPUT: - - ``s`` - a string + - ``s`` -- a string - - ``syms`` - (default: {}) dictionary of - strings to be regarded as symbols or functions ; + - ``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 + - ``accept_sequence`` -- (default: ``False``) controls whether to allow a (possibly nested) set of lists and tuples as input @@ -2562,23 +2574,25 @@ def symbolic_expression_from_string(s, syms=None, accept_sequence=False, *, pars EXAMPLES:: + sage: from sage.calculus.calculus import symbolic_expression_from_string sage: y = var('y') - sage: sage.calculus.calculus.symbolic_expression_from_string('[sin(0)*x^2,3*spam+e^pi]',syms={('spam',0):y},accept_sequence=True) + 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] TESTS: Check that the precision is preserved (:trac:`28814`):: - sage: sage.calculus.calculus.symbolic_expression_from_string(str(RealField(100)(1/3))) + sage: symbolic_expression_from_string(str(RealField(100)(1/3))) 0.3333333333333333333333333333 - sage: sage.calculus.calculus.symbolic_expression_from_string(str(RealField(100)(10^-500/3))) + sage: symbolic_expression_from_string(str(RealField(100)(10^-500/3))) 3.333333333333333333333333333e-501 The Giac interface uses a different parser (:trac:`30133`):: sage: from sage.calculus.calculus import SR_parser_giac - sage: sage.calculus.calculus.symbolic_expression_from_string(repr(giac(SR.var('e'))), parser=SR_parser_giac) + sage: symbolic_expression_from_string(repr(giac(SR.var('e'))), parser=SR_parser_giac) e """ if syms is None: diff --git a/src/sage/calculus/desolvers.py b/src/sage/calculus/desolvers.py index 6dabb5b1121..afe63ccb954 100644 --- a/src/sage/calculus/desolvers.py +++ b/src/sage/calculus/desolvers.py @@ -874,8 +874,8 @@ def desolve_system(des, vars, ics=None, ivar=None, algorithm="maxima"): :: - sage: P1 = plot([solx,soly], (0,1)) - sage: P2 = parametric_plot((solx,soly), (0,1)) + sage: P1 = plot([solx,soly], (0,1)) # needs sage.plot + sage: P2 = parametric_plot((solx,soly), (0,1)) # needs sage.plot Now type ``show(P1)``, ``show(P2)`` to view these plots. @@ -892,7 +892,8 @@ def desolve_system(des, vars, ics=None, ivar=None, algorithm="maxima"): sage: desolve_system([de1, de2], [x1, x2], ics=[1,1], ivar=t) Traceback (most recent call last): ... - ValueError: Initial conditions aren't complete: number of vars is different from number of dependent variables. Got ics = [1, 1], vars = [x1(t), x2(t)] + ValueError: Initial conditions aren't complete: number of vars is different + from number of dependent variables. Got ics = [1, 1], vars = [x1(t), x2(t)] Check that :trac:`9825` is fixed:: @@ -1013,9 +1014,9 @@ def eulers_method(f, x0, y0, h, x1, algorithm="table"): :: sage: pts = eulers_method(5*x+y-5,0,1,1/2,1,algorithm="none") - sage: P1 = list_plot(pts) - sage: P2 = line(pts) - sage: (P1+P2).show() + sage: P1 = list_plot(pts) # needs sage.plot + sage: P2 = line(pts) # needs sage.plot + sage: (P1 + P2).show() # needs sage.plot AUTHORS: @@ -1163,7 +1164,7 @@ def eulers_method_2x2_plot(f, g, t0, x0, y0, h, t1): sage: from sage.calculus.desolvers import eulers_method_2x2_plot sage: f = lambda z : z[2]; g = lambda z : -sin(z[1]) - sage: P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0) + sage: P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0) # needs sage.plot """ from sage.plot.line import line @@ -1430,13 +1431,14 @@ def desolve_system_rk4(des, vars, ics=None, ivar=None, end_points=None, step=0.1 Lotka Volterra system:: sage: from sage.calculus.desolvers import desolve_system_rk4 - sage: x,y,t=var('x y t') - sage: P=desolve_system_rk4([x*(1-y),-y*(1-x)],[x,y],ics=[0,0.5,2],ivar=t,end_points=20) - sage: Q=[ [i,j] for i,j,k in P] - sage: LP=list_plot(Q) + sage: x,y,t = var('x y t') + sage: P = desolve_system_rk4([x*(1-y),-y*(1-x)], [x,y], ics=[0,0.5,2], + ....: ivar=t, end_points=20) + sage: Q = [[i,j] for i,j,k in P] + sage: LP = list_plot(Q) # needs sage.plot - sage: Q=[ [j,k] for i,j,k in P] - sage: LP=list_plot(Q) + sage: Q = [[j,k] for i,j,k in P] + sage: LP = list_plot(Q) # needs sage.plot ALGORITHM: @@ -1570,49 +1572,52 @@ def desolve_odeint(des, ics, times, dvars, ivar=None, compute_jac=False, args=() Lotka Volterra Equations:: sage: from sage.calculus.desolvers import desolve_odeint - sage: x,y=var('x,y') - sage: f=[x*(1-y),-y*(1-x)] - sage: sol=desolve_odeint(f,[0.5,2],srange(0,10,0.1),[x,y]) - sage: p=line(zip(sol[:,0],sol[:,1])) - sage: p.show() + sage: x,y = var('x,y') + sage: f = [x*(1-y), -y*(1-x)] + sage: sol = desolve_odeint(f, [0.5,2], srange(0,10,0.1), [x,y]) # needs scipy + sage: p = line(zip(sol[:,0],sol[:,1])) # needs scipy sage.plot + sage: p.show() # needs scipy sage.plot Lorenz Equations:: - sage: x,y,z=var('x,y,z') + sage: x,y,z = var('x,y,z') sage: # Next we define the parameters - sage: sigma=10 - sage: rho=28 - sage: beta=8/3 + sage: sigma = 10 + sage: rho = 28 + sage: beta = 8/3 sage: # The Lorenz equations - sage: lorenz=[sigma*(y-x),x*(rho-z)-y,x*y-beta*z] + sage: lorenz = [sigma*(y-x),x*(rho-z)-y,x*y-beta*z] sage: # Time and initial conditions - sage: times=srange(0,50.05,0.05) - sage: ics=[0,1,1] - sage: sol=desolve_odeint(lorenz,ics,times,[x,y,z],rtol=1e-13,atol=1e-14) + sage: times = srange(0,50.05,0.05) + sage: ics = [0,1,1] + sage: sol = desolve_odeint(lorenz, ics, times, [x,y,z], # needs scipy + ....: rtol=1e-13, atol=1e-14) One-dimensional stiff system:: - sage: y= var('y') - sage: epsilon=0.01 - sage: f=y^2*(1-y) - sage: ic=epsilon - sage: t=srange(0,2/epsilon,1) - sage: sol=desolve_odeint(f,ic,t,y,rtol=1e-9,atol=1e-10,compute_jac=True) - sage: p=points(zip(t,sol[:,0])) - sage: p.show() + sage: y = var('y') + sage: epsilon = 0.01 + sage: f = y^2*(1-y) + sage: ic = epsilon + sage: t = srange(0,2/epsilon,1) + sage: sol = desolve_odeint(f, ic, t, y, # needs scipy + ....: rtol=1e-9, atol=1e-10, compute_jac=True) + sage: p = points(zip(t,sol[:,0])) # needs scipy sage.plot + sage: p.show() # needs scipy sage.plot Another stiff system with some optional parameters with no default value:: - sage: y1,y2,y3=var('y1,y2,y3') - sage: f1=77.27*(y2+y1*(1-8.375*1e-6*y1-y2)) - sage: f2=1/77.27*(y3-(1+y1)*y2) - sage: f3=0.16*(y1-y3) - sage: f=[f1,f2,f3] - sage: ci=[0.2,0.4,0.7] - sage: t=srange(0,10,0.01) - sage: v=[y1,y2,y3] - sage: sol=desolve_odeint(f,ci,t,v,rtol=1e-3,atol=1e-4,h0=0.1,hmax=1,hmin=1e-4,mxstep=1000,mxords=17) + sage: y1,y2,y3 = var('y1,y2,y3') + sage: f1 = 77.27*(y2+y1*(1-8.375*1e-6*y1-y2)) + sage: f2 = 1/77.27*(y3-(1+y1)*y2) + sage: f3 = 0.16*(y1-y3) + sage: f = [f1,f2,f3] + sage: ci = [0.2,0.4,0.7] + sage: t = srange(0,10,0.01) + sage: v = [y1,y2,y3] + sage: sol = desolve_odeint(f, ci, t, v, rtol=1e-3, atol=1e-4, # needs scipy + ....: h0=0.1, hmax=1, hmin=1e-4, mxstep=1000, mxords=17) AUTHOR: diff --git a/src/sage/calculus/functional.py b/src/sage/calculus/functional.py index 2f769a8f327..03ea8bbd294 100644 --- a/src/sage/calculus/functional.py +++ b/src/sage/calculus/functional.py @@ -1,3 +1,4 @@ +# sage.doctest: needs sage.symbolic """ Functional notation support for common calculus methods diff --git a/src/sage/calculus/functions.py b/src/sage/calculus/functions.py index 6c47b0d2fc9..11a23aac323 100644 --- a/src/sage/calculus/functions.py +++ b/src/sage/calculus/functions.py @@ -1,9 +1,12 @@ +# sage.doctest: needs sage.symbolic r""" Calculus functions """ -from sage.matrix.constructor import matrix +from sage.misc.lazy_import import lazy_import from sage.structure.element import is_Matrix, is_Vector, Expression +lazy_import('sage.matrix.constructor', 'matrix') + from .functional import diff diff --git a/src/sage/calculus/integration.pyx b/src/sage/calculus/integration.pyx index 04ef352a906..19b2dc2e4c8 100644 --- a/src/sage/calculus/integration.pyx +++ b/src/sage/calculus/integration.pyx @@ -1,3 +1,4 @@ +# sage.doctest: needs sage.symbolic """ Numerical Integration diff --git a/src/sage/calculus/interpolation.pyx b/src/sage/calculus/interpolation.pyx index e31d3606cc4..21b82388461 100644 --- a/src/sage/calculus/interpolation.pyx +++ b/src/sage/calculus/interpolation.pyx @@ -55,9 +55,9 @@ cdef class Spline: This example is in the GSL documentation:: - sage: v = [(i + sin(i)/2, i+cos(i^2)) for i in range(10)] + sage: v = [(i + RDF(i).sin()/2, i + RDF(i^2).cos()) for i in range(10)] sage: s = spline(v) - sage: show(point(v) + plot(s,0,9, hue=.8)) + sage: show(point(v) + plot(s,0,9, hue=.8)) # needs sage.plot We compute the area underneath the spline:: @@ -77,13 +77,13 @@ cdef class Spline: We compute the first and second-order derivatives at a few points:: sage: s.derivative(5) - -0.16230085261803... + -0.1623008526180... sage: s.derivative(6) - 0.20997986285714... + 0.2099798628571... sage: s.derivative(5, order=2) - -3.08747074561380... + -3.0874707456138... sage: s.derivative(6, order=2) - 2.61876848274853... + 2.6187684827485... Only the first two derivatives are supported:: @@ -118,7 +118,7 @@ cdef class Spline: Replace `0`-th point, which changes the spline:: - sage: S[0]=(0,1); S + sage: S[0] = (0,1); S [(0, 1), (2, 3), (4, 5)] sage: S(1.5) 2.5 diff --git a/src/sage/calculus/interpolators.pyx b/src/sage/calculus/interpolators.pyx index ded37516d62..fb057f0fc2e 100644 --- a/src/sage/calculus/interpolators.pyx +++ b/src/sage/calculus/interpolators.pyx @@ -1,3 +1,4 @@ +# sage.doctest: needs numpy """ Complex Interpolation @@ -49,9 +50,10 @@ def polygon_spline(pts): sage: ps = polygon_spline(pts) sage: fx = lambda x: ps.value(x).real sage: fy = lambda x: ps.value(x).imag - sage: show(parametric_plot((fx, fy), (0, 2*pi))) - sage: m = Riemann_Map([lambda x: ps.value(real(x))], [lambda x: ps.derivative(real(x))],0) - sage: show(m.plot_colored() + m.plot_spiderweb()) + sage: show(parametric_plot((fx, fy), (0, 2*pi))) # needs sage.plot + sage: m = Riemann_Map([lambda x: ps.value(real(x))], + ....: [lambda x: ps.derivative(real(x))], 0) + sage: show(m.plot_colored() + m.plot_spiderweb()) # needs sage.plot Polygon approximation of an circle:: @@ -119,7 +121,7 @@ cdef class PSpline: sage: ps = polygon_spline(pts) sage: ps.value(.5) (-0.363380227632...-1j) - sage: ps.value(0) - ps.value(2*pi) + sage: ps.value(0) - ps.value(2*RDF.pi()) 0j sage: ps.value(10) (0.26760455264...+1j) @@ -152,6 +154,7 @@ cdef class PSpline: sage: ps = polygon_spline(pts) sage: ps.derivative(1 / 3) (1.27323954473...+0j) + sage: from math import pi sage: ps.derivative(0) - ps.derivative(2*pi) 0j sage: ps.derivative(10) @@ -171,7 +174,7 @@ def complex_cubic_spline(pts): INPUT: - - ``pts`` A list or array of complex numbers, or tuples of the form + - ``pts`` -- A list or array of complex numbers, or tuples of the form `(x,y)`. EXAMPLES: @@ -182,13 +185,16 @@ def complex_cubic_spline(pts): sage: cs = complex_cubic_spline(pts) sage: fx = lambda x: cs.value(x).real sage: fy = lambda x: cs.value(x).imag - sage: show(parametric_plot((fx, fy), (0, 2*pi))) - sage: m = Riemann_Map([lambda x: cs.value(real(x))], [lambda x: cs.derivative(real(x))], 0) - sage: show(m.plot_colored() + m.plot_spiderweb()) + sage: from math import pi + sage: show(parametric_plot((fx, fy), (0, 2*pi))) # needs sage.plot + sage: m = Riemann_Map([lambda x: cs.value(real(x))], + ....: [lambda x: cs.derivative(real(x))], 0) + sage: show(m.plot_colored() + m.plot_spiderweb()) # needs sage.plot Polygon approximation of a circle:: - sage: pts = [e^(I*t / 25) for t in range(25)] + sage: from cmath import exp + sage: pts = [exp(1j * t / 25) for t in range(25)] sage: cs = complex_cubic_spline(pts) sage: cs.derivative(2) (-0.0497765406583...+0.151095006434...j) @@ -273,6 +279,7 @@ cdef class CCSpline: sage: cs = complex_cubic_spline(pts) sage: cs.value(4 / 7) (-0.303961332787...-1.34716728183...j) + sage: from math import pi sage: cs.value(0) - cs.value(2*pi) 0j sage: cs.value(-2.73452) @@ -304,6 +311,7 @@ cdef class CCSpline: sage: cs = complex_cubic_spline(pts) sage: cs.derivative(3 / 5) (1.40578892327...-0.225417136326...j) + sage: from math import pi sage: cs.derivative(0) - cs.derivative(2 * pi) 0j sage: cs.derivative(-6) diff --git a/src/sage/calculus/ode.pyx b/src/sage/calculus/ode.pyx index 74bac5b0811..a7fbb87dffb 100644 --- a/src/sage/calculus/ode.pyx +++ b/src/sage/calculus/ode.pyx @@ -109,30 +109,31 @@ cdef int c_f(double t,double* y, double* dydt,void *params): class ode_solver(): r""" - :meth:`ode_solver` is a class that wraps the GSL libraries ode - solver routines To use it instantiate a class,:: + :meth:`ode_solver` is a class that wraps the GSL library's ode solver routines. - sage: T=ode_solver() + To use it, instantiate the class:: - To solve a system of the form ``dy_i/dt=f_i(t,y)``, you must + sage: T = ode_solver() + + To solve a system of the form `dy_i/dt=f_i(t,y)`, you must supply a vector or tuple/list valued function ``f`` representing - ``f_i``. The functions ``f`` and the jacobian should have the + `f_i`. The functions `f` and the jacobian should have the form ``foo(t,y)`` or ``foo(t,y,params)``. ``params`` which is optional allows for your function to depend on one or a tuple of parameters. Note if you use it, ``params`` must be a tuple even if it only has one component. For example if you wanted to solve - `y''+y=0`. You need to write it as a first order system:: + `y''+y=0`, you would need to write it as a first order system:: y_0' = y_1 y_1' = -y_0 In code:: - sage: f = lambda t,y:[y[1],-y[0]] - sage: T.function=f + sage: f = lambda t, y: [y[1], -y[0]] + sage: T.function = f - For some algorithms the jacobian must be supplied as well, the - form of this should be a function return a list of lists of the + For some algorithms, the jacobian must be supplied as well, the + form of this should be a function returning a list of lists of the form ``[ [df_1/dy_1,...,df_1/dy_n], ..., [df_n/dy_1,...,df_n,dy_n], [df_1/dt,...,df_n/dt] ]``. @@ -143,45 +144,45 @@ class ode_solver(): There are a variety of algorithms available for different types of systems. Possible algorithms are - - ``rkf45`` - runga-kutta-felhberg (4,5) + - ``'rkf45'`` -- Runge-Kutta-Fehlberg (4,5) - - ``rk2`` - embedded runga-kutta (2,3) + - ``'rk2'`` -- embedded Runge-Kutta (2,3) - - ``rk4`` - 4th order classical runga-kutta + - ``'rk4'`` -- 4th order classical Runge-Kutta - - ``rk8pd`` - runga-kutta prince-dormand (8,9) + - ``'rk8pd'`` -- Runge-Kutta Prince-Dormand (8,9) - - ``rk2imp`` - implicit 2nd order runga-kutta at gaussian points + - ``'rk2imp'`` -- implicit 2nd order Runge-Kutta at gaussian points - - ``rk4imp`` - implicit 4th order runga-kutta at gaussian points + - ``'rk4imp'`` -- implicit 4th order Runge-Kutta at gaussian points - - ``bsimp`` - implicit burlisch-stoer (requires jacobian) + - ``'bsimp'`` -- implicit Burlisch-Stoer (requires jacobian) - - ``gear1`` - M=1 implicit gear + - ``'gear1'`` -- M=1 implicit gear - - ``gear2`` - M=2 implicit gear + - ``'gear2'`` -- M=2 implicit gear - The default algorithm is ``rkf45``. If you instead wanted to use - ``bsimp`` you would do:: + The default algorithm is ``'rkf45'``. If you instead wanted to use + ``'bsimp'`` you would do:: - sage: T.algorithm="bsimp" + sage: T.algorithm = "bsimp" - The user should supply initial conditions in y_0. For example if - your initial conditions are y_0=1,y_1=1, do:: + The user should supply initial conditions in ``y_0``. For example if + your initial conditions are `y_0=1, y_1=1`, do:: - sage: T.y_0=[1,1] + sage: T.y_0 = [1,1] The actual solver is invoked by the method :meth:`ode_solve`. It has arguments ``t_span``, ``y_0``, ``num_points``, ``params``. ``y_0`` must be supplied either as an argument or above by assignment. Params which are optional and only necessary if your - system uses params can be supplied to ``ode_solve`` or by + system uses ``params`` can be supplied to ``ode_solve`` or by assignment. ``t_span`` is the time interval on which to solve the ode. There are two ways to specify ``t_span``: - * If ``num_points`` is not specified then the sequence ``t_span`` + * If ``num_points`` is not specified, then the sequence ``t_span`` is used as the time points for the solution. Note that the first element ``t_span[0]`` is the initial time, where the initial condition ``y_0`` is the specified solution, and @@ -192,10 +193,10 @@ class ode_solver(): and the solution will be computed at ``num_points`` equally spaced points between ``t_span[0]`` and ``t_span[1]``. The initial condition is also included in the output so that - ``num_points``\ +1 total points are returned. E.g. if ``t_span + ``num_points + 1`` total points are returned. E.g. if ``t_span = [0.0, 1.0]`` and ``num_points = 10``, then solution is returned at the 11 time points ``[0.0, 0.1, 0.2, 0.3, 0.4, 0.5, - 0.6, 0.7, 0.8, 0.9, 1.0]``\ . + 0.6, 0.7, 0.8, 0.9, 1.0]``. (Note that if ``num_points`` is specified and ``t_span`` is not length 2 then ``t_span`` are used as the time points and @@ -203,13 +204,19 @@ class ode_solver(): Error is estimated via the expression ``D_i = error_abs*s_i+error_rel*(a|y_i|+a_dydt*h*|y_i'|)``. The user can - specify ``error_abs`` (1e-10 by default), ``error_rel`` (1e-10 by - default) ``a`` (1 by default), ``a_(dydt)`` (0 by default) and - ``s_i`` (as scaling_abs which should be a tuple and is 1 in all - components by default). If you specify one of ``a`` or ``a_dydt`` + specify + + - ``error_abs`` (1e-10 by default), + - ``error_rel`` (1e-10 by default), + - ``a`` (1 by default), + - ``a_dydt`` (0 by default) and + - ``s_i`` (as ``scaling_abs`` which should be a tuple and is 1 in all + components by default). + + If you specify one of ``a`` or ``a_dydt`` you must specify the other. You may specify ``a`` and ``a_dydt`` without ``scaling_abs`` (which will be taken =1 be default). - ``h`` is the initial step size which is (1e-2) by default. + ``h`` is the initial step size, which is 1e-2 by default. ``ode_solve`` solves the solution as a list of tuples of the form, ``[ (t_0,[y_1,...,y_n]),(t_1,[y_1,...,y_n]),...,(t_n,[y_1,...,y_n])]``. @@ -223,27 +230,29 @@ class ode_solver(): Consider solving the Van der Pol oscillator `x''(t) + ux'(t)(x(t)^2-1)+x(t)=0` between `t=0` and `t= 100`. As a first order system it is `x'=y`, `y'=-x+uy(1-x^2)`. Let us take `u=10` - and use initial conditions `(x,y)=(1,0)` and use the runga-kutta - prince-dormand algorithm. :: + and use initial conditions `(x,y)=(1,0)` and use the Runge-Kutta + Prince-Dormand algorithm. :: - sage: def f_1(t,y,params): - ....: return[y[1],-y[0]-params[0]*y[1]*(y[0]**2-1.0)] + sage: def f_1(t, y, params): + ....: return [y[1], -y[0] - params[0]*y[1]*(y[0]**2-1.0)] - sage: def j_1(t,y,params): - ....: return [ [0.0, 1.0],[-2.0*params[0]*y[0]*y[1]-1.0,-params[0]*(y[0]*y[0]-1.0)], [0.0, 0.0] ] + sage: def j_1(t, y, params): + ....: return [[0.0, 1.0], + ....: [-2.0*params[0]*y[0]*y[1] - 1.0, -params[0]*(y[0]*y[0]-1.0)], + ....: [0.0, 0.0]] - sage: T=ode_solver() - sage: T.algorithm="rk8pd" - sage: T.function=f_1 - sage: T.jacobian=j_1 - sage: T.ode_solve(y_0=[1,0],t_span=[0,100],params=[10.0],num_points=1000) + sage: T = ode_solver() + sage: T.algorithm = "rk8pd" + sage: T.function = f_1 + sage: T.jacobian = j_1 + sage: T.ode_solve(y_0=[1,0], t_span=[0,100], params=[10.0], num_points=1000) sage: import tempfile - sage: with tempfile.NamedTemporaryFile(suffix=".png") as f: + sage: with tempfile.NamedTemporaryFile(suffix=".png") as f: # needs sage.plot ....: T.plot_solution(filename=f.name) The solver line is equivalent to:: - sage: T.ode_solve(y_0=[1,0],t_span=[x/10.0 for x in range(1000)],params = [10.0]) + sage: T.ode_solve(y_0=[1,0], t_span=[x/10.0 for x in range(1000)], params=[10.0]) Let's try a system:: @@ -254,40 +263,42 @@ class ode_solver(): We will not use the jacobian this time and will change the error tolerances. :: - sage: g_1= lambda t,y: [y[1]*y[2],-y[0]*y[2],-0.51*y[0]*y[1]] - sage: T.function=g_1 - sage: T.y_0=[0,1,1] - sage: T.scale_abs=[1e-4,1e-4,1e-5] - sage: T.error_rel=1e-4 - sage: T.ode_solve(t_span=[0,12],num_points=100) + sage: g_1 = lambda t,y: [y[1]*y[2], -y[0]*y[2], -0.51*y[0]*y[1]] + sage: T.function = g_1 + sage: T.y_0 = [0,1,1] + sage: T.scale_abs = [1e-4, 1e-4, 1e-5] + sage: T.error_rel = 1e-4 + sage: T.ode_solve(t_span=[0,12], num_points=100) - By default T.plot_solution() plots the y_0, to plot general y_i use:: + By default ``T.plot_solution()`` plots the `y_0`; to plot general `y_i`, use:: - sage: with tempfile.NamedTemporaryFile(suffix=".png") as f: + sage: with tempfile.NamedTemporaryFile(suffix=".png") as f: # needs sage.plot ....: T.plot_solution(i=0, filename=f.name) ....: T.plot_solution(i=1, filename=f.name) ....: T.plot_solution(i=2, filename=f.name) The method interpolate_solution will return a spline interpolation - through the points found by the solver. By default y_0 is - interpolated. You can interpolate y_i through the keyword - argument i. :: + through the points found by the solver. By default, `y_0` is + interpolated. You can interpolate `y_i` through the keyword + argument ``i``. :: sage: f = T.interpolate_solution() - sage: plot(f,0,12).show() + sage: plot(f,0,12).show() # needs sage.plot sage: f = T.interpolate_solution(i=1) - sage: plot(f,0,12).show() + sage: plot(f,0,12).show() # needs sage.plot sage: f = T.interpolate_solution(i=2) - sage: plot(f,0,12).show() + sage: plot(f,0,12).show() # needs sage.plot sage: f = T.interpolate_solution() + sage: from math import pi sage: f(pi) 0.5379... The solver attributes may also be set up using arguments to ode_solver. The previous example can be rewritten as:: - sage: T = ode_solver(g_1,y_0=[0,1,1],scale_abs=[1e-4,1e-4,1e-5],error_rel=1e-4, algorithm="rk8pd") - sage: T.ode_solve(t_span=[0,12],num_points=100) + sage: T = ode_solver(g_1, y_0=[0,1,1], scale_abs=[1e-4,1e-4,1e-5], + ....: error_rel=1e-4, algorithm="rk8pd") + sage: T.ode_solve(t_span=[0,12], num_points=100) sage: f = T.interpolate_solution() sage: f(pi) 0.5379... @@ -295,7 +306,7 @@ class ode_solver(): Unfortunately because Python functions are used, this solver is slow on systems that require many function evaluations. It is possible to pass a compiled function by deriving from the - class ``ode_sysem`` and overloading ``c_f`` and ``c_j`` with C + class :class:`ode_system` and overloading ``c_f`` and ``c_j`` with C functions that specify the system. The following will work in the notebook: @@ -382,11 +393,11 @@ class ode_solver(): sage: T.function = lambda t,y: [cos(y[0]) * sin(t)] sage: T.jacobian = lambda t,y: [[-sin(y[0]) * sin(t)]] sage: T.ode_solve(y_0=[1],t_span=[0,20],num_points=1000) - sage: T.plot_solution() + sage: T.plot_solution() # needs sage.plot And with some options:: - sage: T.plot_solution(color='red', axes_labels=["t", "x(t)"]) + sage: T.plot_solution(color='red', axes_labels=["t", "x(t)"]) # needs sage.plot """ if interpolate: from sage.plot.line import line2d diff --git a/src/sage/calculus/riemann.pyx b/src/sage/calculus/riemann.pyx index 46e2964eade..05d5b1ec60a 100644 --- a/src/sage/calculus/riemann.pyx +++ b/src/sage/calculus/riemann.pyx @@ -1,3 +1,4 @@ +# sage.doctest: needs numpy sage.symbolic """ Riemann Mapping @@ -1190,9 +1191,10 @@ cpdef complex_to_spiderweb(np.ndarray[COMPLEX_T, ndim = 2] z_values, sage: from sage.calculus.riemann import complex_to_spiderweb sage: import numpy - sage: zval = numpy.array([[0, 1, 1000],[.2+.3j,1,-.3j],[0,0,0]],dtype = numpy.complex128) + sage: zval = numpy.array([[0,1,1000], [.2+.3j,1,-.3j], [0,0,0]], + ....: dtype=numpy.complex128) sage: deriv = numpy.array([[.1]],dtype = numpy.float64) - sage: complex_to_spiderweb(zval, deriv,deriv, 4,4,[0,0,0],1,False,0.001) + sage: complex_to_spiderweb(zval, deriv, deriv, 4, 4, [0,0,0], 1, False, 0.001) array([[[1., 1., 1.], [1., 1., 1.], [1., 1., 1.]], @@ -1205,7 +1207,7 @@ cpdef complex_to_spiderweb(np.ndarray[COMPLEX_T, ndim = 2] z_values, [1., 1., 1.], [1., 1., 1.]]]) - sage: complex_to_spiderweb(zval, deriv,deriv, 4,4,[0,0,0],1,True,0.001) + sage: complex_to_spiderweb(zval, deriv, deriv, 4, 4, [0,0,0], 1, True, 0.001) array([[[1. , 1. , 1. ], [1. , 0.05558355, 0.05558355], [0.17301243, 0. , 0. ]], @@ -1280,12 +1282,12 @@ cpdef complex_to_rgb(np.ndarray[COMPLEX_T, ndim = 2] z_values): sage: from sage.calculus.riemann import complex_to_rgb sage: import numpy - sage: complex_to_rgb(numpy.array([[0, 1, 1000]], dtype = numpy.complex128)) + sage: complex_to_rgb(numpy.array([[0, 1, 1000]], dtype=numpy.complex128)) array([[[1. , 1. , 1. ], [1. , 0.05558355, 0.05558355], [0.17301243, 0. , 0. ]]]) - sage: complex_to_rgb(numpy.array([[0, 1j, 1000j]], dtype = numpy.complex128)) + sage: complex_to_rgb(numpy.array([[0, 1j, 1000j]], dtype=numpy.complex128)) array([[[1. , 1. , 1. ], [0.52779177, 1. , 0.05558355], [0.08650622, 0.17301243, 0. ]]]) diff --git a/src/sage/calculus/test_sympy.py b/src/sage/calculus/test_sympy.py index 927e6ee4fb6..b99d5d7857e 100644 --- a/src/sage/calculus/test_sympy.py +++ b/src/sage/calculus/test_sympy.py @@ -1,4 +1,4 @@ -# -*- coding: utf-8 -*- +# sage.doctest: needs sage.symbolic r""" A Sample Session using SymPy @@ -102,46 +102,47 @@ And here are some actual tests of sympy:: - sage: from sympy import Symbol, cos, sympify, pprint - sage: from sympy.abc import x + sage: from sympy import Symbol, cos, sympify, pprint # needs sympy + sage: from sympy.abc import x # needs sympy :: - sage: e = (1/cos(x)^3)._sympy_(); e + sage: e = (1/cos(x)^3)._sympy_(); e # needs sympy cos(x)**(-3) - sage: f = e.series(x, 0, int(10)); f + sage: f = e.series(x, 0, int(10)); f # needs sympy 1 + 3*x**2/2 + 11*x**4/8 + 241*x**6/240 + 8651*x**8/13440 + O(x**10) And the pretty-printer. Since unicode characters are not working on some architectures, we disable it:: - sage: from sympy.printing import pprint_use_unicode - sage: prev_use = pprint_use_unicode(False) - sage: pprint(e) + sage: from sympy.printing import pprint_use_unicode # needs sympy + sage: prev_use = pprint_use_unicode(False) # needs sympy + sage: pprint(e) # needs sympy 1 ------- 3 cos (x) - sage: pprint(f) + sage: pprint(f) # needs sympy 2 4 6 8 3*x 11*x 241*x 8651*x / 10\ 1 + ---- + ----- + ------ + ------- + O\x / 2 8 240 13440 - sage: pprint_use_unicode(prev_use) + sage: pprint_use_unicode(prev_use) # needs sympy False And the functionality to convert from sympy format to Sage format:: - sage: e._sage_() + sage: e._sage_() # needs sympy cos(x)^(-3) - sage: e._sage_().taylor(x._sage_(), 0, 8) + sage: e._sage_().taylor(x._sage_(), 0, 8) # needs sympy 8651/13440*x^8 + 241/240*x^6 + 11/8*x^4 + 3/2*x^2 + 1 - sage: f._sage_() + sage: f._sage_() # needs sympy 8651/13440*x^8 + 241/240*x^6 + 11/8*x^4 + 3/2*x^2 + Order(x^10) + 1 Mixing SymPy with Sage:: + sage: # needs sympy sage: import sympy sage: var("x")._sympy_() + var("y")._sympy_() x + y @@ -155,7 +156,7 @@ sage: t1, t2 (omega + x, omega + x) - sage: e=sympy.sin(var("y"))+sage.all.cos(sympy.Symbol("x")) + sage: e = sympy.sin(var("y"))+sage.all.cos(sympy.Symbol("x")) sage: type(e) sage: e @@ -172,23 +173,24 @@ :: - sage: a = sympy.Matrix([1, 2, 3]) - sage: a[1] + sage: a = sympy.Matrix([1, 2, 3]) # needs sympy + sage: a[1] # needs sympy 2 :: - sage: sympify(1.5) + sage: sympify(1.5) # needs sympy 1.50000000000000 - sage: sympify(2) + sage: sympify(2) # needs sympy 2 - sage: sympify(-2) + sage: sympify(-2) # needs sympy -2 TESTS: This was fixed in Sympy, see :trac:`14437`:: + sage: # needs sympy sage: from sympy import Function, Symbol, rsolve sage: u = Function('u') sage: n = Symbol('n', integer=True) diff --git a/src/sage/calculus/transforms/all.py b/src/sage/calculus/transforms/all.py index 06ed525b868..379b3b69c37 100644 --- a/src/sage/calculus/transforms/all.py +++ b/src/sage/calculus/transforms/all.py @@ -1,3 +1,5 @@ -from .fft import FastFourierTransform, FFT -from .dwt import WaveletTransform, DWT +from sage.misc.lazy_import import lazy_import + +lazy_import("sage.calculus.transforms.fft", ["FastFourierTransform", "FFT"]) +lazy_import("sage.calculus.transforms.dwt", ["WaveletTransform", "DWT"]) from .dft import IndexedSequence diff --git a/src/sage/calculus/transforms/dft.py b/src/sage/calculus/transforms/dft.py index db18a20e129..b413bc0ea81 100644 --- a/src/sage/calculus/transforms/dft.py +++ b/src/sage/calculus/transforms/dft.py @@ -16,7 +16,7 @@ - plotting, printing -- :meth:`IndexedSequence.plot`, :meth:`IndexedSequence.plot_histogram`, :meth:`_repr_`, :meth:`__str__` -- dft -- computes the discrete Fourier transform for the following cases: +- :meth:`dft` -- computes the discrete Fourier transform for the following cases: * a sequence (over `\QQ` or :class:`CyclotomicField`) indexed by ``range(N)`` or `\ZZ / N \ZZ` @@ -26,21 +26,21 @@ * a sequence (as above) indexed by a complete set of representatives of the conjugacy classes of a finite matrix group -- idft -- computes the discrete Fourier transform for the following cases: +- :meth:`idft` -- computes the discrete Fourier transform for the following cases: * a sequence (over `\QQ` or CyclotomicField) indexed by ``range(N)`` or `\ZZ / N \ZZ` -- dct, dst (for discrete Fourier/Cosine/Sine transform) +- :meth:`dct`, :meth:`dst` (for discrete Fourier/Cosine/Sine transform) - convolution (in :meth:`IndexedSequence.convolution` and :meth:`IndexedSequence.convolution_periodic`) -- fft, ifft -- (fast Fourier transforms) wrapping GSL's +- :meth:`fft`, :meth:`ifft` -- (fast Fourier transforms) wrapping GSL's ``gsl_fft_complex_forward()``, ``gsl_fft_complex_inverse()``, using William Stein's :func:`FastFourierTransform` -- dwt, idwt -- (fast wavelet transforms) wrapping GSL's ``gsl_dwt_forward()``, +- :meth:`dwt`, :meth:`idwt` -- (fast wavelet transforms) wrapping GSL's ``gsl_dwt_forward()``, ``gsl_dwt_backward()`` using Joshua Kantor's :func:`WaveletTransform` class. Allows for wavelets of type: @@ -71,20 +71,20 @@ # # https://www.gnu.org/licenses/ ########################################################################## -from sage.rings.number_field.number_field import CyclotomicField +from sage.functions.trig import sin, cos from sage.misc.lazy_import import lazy_import -from sage.groups.abelian_gps.abelian_group import AbelianGroup -from sage.groups.perm_gps.permgroup_element import is_PermutationGroupElement -from sage.rings.integer_ring import ZZ from sage.rings.integer import Integer +from sage.rings.integer_ring import ZZ from sage.rings.rational_field import QQ -from sage.rings.real_mpfr import RR -from sage.functions.all import sin, cos -from sage.calculus.transforms.fft import FastFourierTransform -from sage.calculus.transforms.dwt import WaveletTransform from sage.structure.sage_object import SageObject from sage.structure.sequence import Sequence + +lazy_import("sage.calculus.transforms.dwt", "WaveletTransform") +lazy_import("sage.calculus.transforms.fft", "FastFourierTransform") +lazy_import("sage.groups.abelian_gps.abelian_group", "AbelianGroup") +lazy_import("sage.groups.perm_gps.permgroup_element", "PermutationGroupElement") lazy_import("sage.plot.all", ["polygon", "line", "text"]) +lazy_import("sage.rings.number_field.number_field", "CyclotomicField") class IndexedSequence(SageObject): @@ -218,8 +218,7 @@ def _repr_(self): indexed by [0, 1, 2] sage: I = GF(3) sage: A = [i^2 for i in I] - sage: s = IndexedSequence(A,I) - sage: s + sage: s = IndexedSequence(A,I); s Indexed sequence: [0, 1, 1] indexed by Finite Field of size 3 """ @@ -240,9 +239,10 @@ def plot_histogram(self, clr=(0, 0, 1), eps=0.4): sage: J = list(range(3)) sage: A = [ZZ(i^2)+1 for i in J] sage: s = IndexedSequence(A,J) - sage: P = s.plot_histogram() - sage: show(P) # Not tested + sage: P = s.plot_histogram() # needs sage.plot + sage: show(P) # not tested # needs sage.plot """ + from sage.rings.real_mpfr import RR # elements must be coercible into RR I = self.index_object() N = len(I) @@ -268,9 +268,10 @@ def plot(self): sage: I = list(range(3)) sage: A = [ZZ(i^2)+1 for i in I] sage: s = IndexedSequence(A,I) - sage: P = s.plot() - sage: show(P) # Not tested + sage: P = s.plot() # needs sage.plot + sage: show(P) # not tested # needs sage.plot """ + from sage.rings.real_mpfr import RR # elements must be coercible into RR I = self.index_object() S = self.list() @@ -286,30 +287,44 @@ def dft(self, chi=lambda x: x): sage: J = list(range(6)) sage: A = [ZZ(1) for i in J] sage: s = IndexedSequence(A,J) - sage: s.dft(lambda x:x^2) + sage: s.dft(lambda x: x^2) # needs sage.rings.number_field Indexed sequence: [6, 0, 0, 6, 0, 0] indexed by [0, 1, 2, 3, 4, 5] - sage: s.dft() + sage: s.dft() # needs sage.rings.number_field Indexed sequence: [6, 0, 0, 0, 0, 0] indexed by [0, 1, 2, 3, 4, 5] + + sage: # needs sage.groups sage: G = SymmetricGroup(3) sage: J = G.conjugacy_classes_representatives() - sage: s = IndexedSequence([1,2,3],J) # 1,2,3 are the values of a class fcn on G + sage: s = IndexedSequence([1,2,3], J) # 1,2,3 are the values of a class fcn on G sage: s.dft() # the "scalar-valued Fourier transform" of this class fcn Indexed sequence: [8, 2, 2] indexed by [(), (1,2), (1,2,3)] - sage: J = AbelianGroup(2,[2,3],names='ab') - sage: s = IndexedSequence([1,2,3,4,5,6],J) + sage: J = AbelianGroup(2, [2,3], names='ab') + sage: s = IndexedSequence([1,2,3,4,5,6], J) sage: s.dft() # the precision of output is somewhat random and architecture dependent. - Indexed sequence: [21.0000000000000, -2.99999999999997 - 1.73205080756885*I, -2.99999999999999 + 1.73205080756888*I, -9.00000000000000 + 0.0000000000000485744257349999*I, -0.00000000000000976996261670137 - 0.0000000000000159872115546022*I, -0.00000000000000621724893790087 - 0.0000000000000106581410364015*I] - indexed by Multiplicative Abelian group isomorphic to C2 x C3 + Indexed sequence: [21.0000000000000, + -2.99999999999997 - 1.73205080756885*I, + -2.99999999999999 + 1.73205080756888*I, + -9.00000000000000 + 0.0000000000000485744257349999*I, + -0.00000000000000976996261670137 - 0.0000000000000159872115546022*I, + -0.00000000000000621724893790087 - 0.0000000000000106581410364015*I] + indexed by Multiplicative Abelian group isomorphic to C2 x C3 sage: J = CyclicPermutationGroup(6) - sage: s = IndexedSequence([1,2,3,4,5,6],J) + sage: s = IndexedSequence([1,2,3,4,5,6], J) sage: s.dft() # the precision of output is somewhat random and architecture dependent. - Indexed sequence: [21.0000000000000, -2.99999999999997 - 1.73205080756885*I, -2.99999999999999 + 1.73205080756888*I, -9.00000000000000 + 0.0000000000000485744257349999*I, -0.00000000000000976996261670137 - 0.0000000000000159872115546022*I, -0.00000000000000621724893790087 - 0.0000000000000106581410364015*I] - indexed by Cyclic group of order 6 as a permutation group + Indexed sequence: [21.0000000000000, + -2.99999999999997 - 1.73205080756885*I, + -2.99999999999999 + 1.73205080756888*I, + -9.00000000000000 + 0.0000000000000485744257349999*I, + -0.00000000000000976996261670137 - 0.0000000000000159872115546022*I, + -0.00000000000000621724893790087 - 0.0000000000000106581410364015*I] + indexed by Cyclic group of order 6 as a permutation group + + sage: # needs sage.rings.number_field sage: p = 7; J = list(range(p)); A = [kronecker_symbol(j,p) for j in J] - sage: s = IndexedSequence(A,J) + sage: s = IndexedSequence(A, J) sage: Fs = s.dft() sage: c = Fs.list()[1]; [x/c for x in Fs.list()]; s.list() [0, 1, 1, -1, 1, -1, -1] @@ -336,7 +351,7 @@ def dft(self, chi=lambda x: x): zeta = CyclotomicField(N).gen() FT = [sum([S[i] * chi(zeta**(i * j)) for i in J]) for j in J] elif (J[0] not in ZZ) and G.is_abelian() and F == ZZ or (F.is_field() and F.base_ring() == QQ): - if is_PermutationGroupElement(J[0]): + if isinstance(J[0], PermutationGroupElement): # J is a CyclicPermGp n = G.order() a = list(n.factor()) @@ -364,13 +379,13 @@ def idft(self): sage: J = list(range(5)) sage: A = [ZZ(1) for i in J] sage: s = IndexedSequence(A,J) - sage: fs = s.dft(); fs + sage: fs = s.dft(); fs # needs sage.rings.number_field Indexed sequence: [5, 0, 0, 0, 0] indexed by [0, 1, 2, 3, 4] - sage: it = fs.idft(); it + sage: it = fs.idft(); it # needs sage.rings.number_field Indexed sequence: [1, 1, 1, 1, 1] indexed by [0, 1, 2, 3, 4] - sage: it == s + sage: it == s # needs sage.rings.number_field True """ F = self.base_ring() # elements must be coercible into QQ(zeta_N) @@ -390,18 +405,23 @@ def dct(self): EXAMPLES:: sage: J = list(range(5)) - sage: A = [exp(-2*pi*i*I/5) for i in J] - sage: s = IndexedSequence(A,J) - sage: s.dct() + sage: A = [exp(-2*pi*i*I/5) for i in J] # needs sage.symbolic + sage: s = IndexedSequence(A, J) # needs sage.symbolic + sage: s.dct() # needs sage.symbolic Indexed sequence: [0, 1/16*(sqrt(5) + I*sqrt(-2*sqrt(5) + 10) + ... indexed by [0, 1, 2, 3, 4] """ - from sage.symbolic.constants import pi F = self.base_ring() # elements must be coercible into RR + try: + pi = F.pi() + except AttributeError: + from sage.symbolic.constants import pi + pi = F(pi) + J = self.index_object() # must be = range(N) N = len(J) S = self.list() - PI = 2 * F(pi) / N + PI = 2 * pi / N FT = [sum([S[i] * cos(PI * i * j) for i in J]) for j in J] return IndexedSequence(FT, J) @@ -412,16 +432,21 @@ def dst(self): EXAMPLES:: sage: J = list(range(5)) - sage: I = CC.0; pi = CC(pi) + sage: I = CC.0; pi = CC.pi() sage: A = [exp(-2*pi*i*I/5) for i in J] - sage: s = IndexedSequence(A,J) + sage: s = IndexedSequence(A, J) sage: s.dst() # discrete sine Indexed sequence: [0.000000000000000, 1.11022302462516e-16 - 2.50000000000000*I, ...] indexed by [0, 1, 2, 3, 4] """ - from sage.symbolic.constants import pi F = self.base_ring() # elements must be coercible into RR + try: + pi = F.pi() + except AttributeError: + from sage.symbolic.constants import pi + pi = F(pi) + J = self.index_object() # must be = range(N) N = len(J) S = self.list() @@ -711,6 +736,7 @@ def dwt(self, other="haar", wavelet_k=2): Indexed sequence: [2.82842712474999, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000] indexed by [0, 1, 2, 3, 4, 5, 6, 7] """ + from sage.rings.real_mpfr import RR # elements must be coercible into RR J = self.index_object() # must be = range(N) N = len(J) # must be 1 minus a power of 2 @@ -788,6 +814,7 @@ def idwt(self, other="haar", wavelet_k=2): sage: t.idwt("bspline", 103) == s True """ + from sage.rings.real_mpfr import RR # elements must be coercible into RR J = self.index_object() # must be = range(N) N = len(J) # must be 1 minus a power of 2 diff --git a/src/sage/calculus/transforms/dwt.pyx b/src/sage/calculus/transforms/dwt.pyx index 9296afb7824..807169f2886 100644 --- a/src/sage/calculus/transforms/dwt.pyx +++ b/src/sage/calculus/transforms/dwt.pyx @@ -64,27 +64,28 @@ def WaveletTransform(n, wavelet_type, wavelet_k): sage: for i in range(1, 11): ....: a[i] = 1 ....: a[128-i] = 1 - sage: a.plot().show(ymin=0) + sage: a.plot().show(ymin=0) # needs sage.plot sage: a.forward_transform() - sage: a.plot().show() + sage: a.plot().show() # needs sage.plot sage: a = WaveletTransform(128,'haar',2) sage: for i in range(1, 11): a[i] = 1; a[128-i] = 1 sage: a.forward_transform() - sage: a.plot().show(ymin=0) + sage: a.plot().show(ymin=0) # needs sage.plot sage: a = WaveletTransform(128,'bspline_centered',103) sage: for i in range(1, 11): a[i] = 1; a[100+i] = 1 sage: a.forward_transform() - sage: a.plot().show(ymin=0) + sage: a.plot().show(ymin=0) # needs sage.plot This example gives a simple example of wavelet compression:: + sage: # needs sage.symbolic sage: a = DWT(2048,'daubechies',6) sage: for i in range(2048): a[i]=float(sin((i*5/2048)**2)) - sage: a.plot().show() # long time (7s on sage.math, 2011) + sage: a.plot().show() # long time (7s on sage.math, 2011), needs sage.plot sage: a.forward_transform() sage: for i in range(1800): a[2048-i-1] = 0 sage: a.backward_transform() - sage: a.plot().show() # long time (7s on sage.math, 2011) + sage: a.plot().show() # long time (7s on sage.math, 2011), needs sage.plot """ cdef size_t _n, _k _n = int(n) diff --git a/src/sage/calculus/transforms/fft.pyx b/src/sage/calculus/transforms/fft.pyx index 05db8f5e3e7..168698957c4 100644 --- a/src/sage/calculus/transforms/fft.pyx +++ b/src/sage/calculus/transforms/fft.pyx @@ -22,7 +22,6 @@ AUTHORS: from cysignals.memory cimport sig_malloc, sig_free -import sage.libs.pari.all from sage.rings.integer import Integer from sage.rings.complex_mpfr import ComplexNumber @@ -70,9 +69,9 @@ def FastFourierTransform(size, base_ring=None): ....: a[128-i] = 1 sage: a[:6:2] [(0.0, 0.0), (1.0, 0.0), (1.0, 0.0)] - sage: a.plot().show(ymin=0) + sage: a.plot().show(ymin=0) # needs sage.plot sage: a.forward_transform() - sage: a.plot().show() + sage: a.plot().show() # needs sage.plot """ return FastFourierTransform_complex(int(size)) @@ -152,6 +151,7 @@ cdef class FastFourierTransform_complex(FastFourierTransform_base): EXAMPLES:: + sage: # needs sage.rings.mpfr sage.symbolic sage: I = CC(I) sage: a = FastFourierTransform(4) sage: a[0] = 1 @@ -243,18 +243,19 @@ cdef class FastFourierTransform_complex(FastFourierTransform_base): EXAMPLES:: sage: a = FastFourierTransform(4) - sage: a._plot_polar(0,2) + sage: a._plot_polar(0,2) # needs sage.plot Graphics object consisting of 2 graphics primitives """ from sage.plot.point import point + from sage.symbolic.constants import pi, I cdef int i v = [] - pi = sage.symbolic.constants.pi.n() - I = sage.symbolic.constants.I.n() - s = 1/(3*pi) # so arg gets scaled between -1/3 and 1/3 + pi = pi.n() + I = I.n() + s = 1/(3*pi) # so arg gets scaled between -1/3 and 1/3. for i from xmin <= i < xmax: z = self.data[2*i] + I*self.data[2*i+1] @@ -282,15 +283,15 @@ cdef class FastFourierTransform_complex(FastFourierTransform_base): EXAMPLES:: sage: a = FastFourierTransform(4) - sage: a._plot_rect(0,3) + sage: a._plot_rect(0,3) # needs sage.plot Graphics object consisting of 3 graphics primitives """ + from sage.plot.point import point + cdef int i cdef double x, h v = [] - point = sage.plot.all.point - for i in range(xmin, xmax): x = self.data[2*i] h = self.data[2*i+1] @@ -302,10 +303,12 @@ cdef class FastFourierTransform_complex(FastFourierTransform_base): Plot a slice of the array. - ``style`` -- Style of the plot, options are ``"rect"`` or ``"polar"`` - - ``rect`` -- height represents real part, color represents - imaginary part. - - ``polar`` -- height represents absolute value, color - represents argument. + + - ``rect`` -- height represents real part, color represents + imaginary part. + - ``polar`` -- height represents absolute value, color + represents argument. + - ``xmin`` -- The lower bound of the slice to plot. 0 by default. - ``xmax`` -- The upper bound of the slice to plot. ``len(self)`` by default. - ``**args`` -- passed on to the line plotting function. @@ -318,11 +321,11 @@ cdef class FastFourierTransform_complex(FastFourierTransform_base): sage: a = FastFourierTransform(16) sage: for i in range(16): a[i] = (random(),random()) - sage: A = plot(a) - sage: B = plot(a, style='polar') - sage: type(A) + sage: A = plot(a) # needs sage.plot + sage: B = plot(a, style='polar') # needs sage.plot + sage: type(A) # needs sage.plot - sage: type(B) + sage: type(B) # needs sage.plot sage: a = FastFourierTransform(125) sage: b = FastFourierTransform(125) @@ -330,7 +333,7 @@ cdef class FastFourierTransform_complex(FastFourierTransform_base): sage: for i in range(1, 60): b[i]=1 sage: a.forward_transform() sage: a.inverse_transform() - sage: (a.plot()+b.plot()) + sage: a.plot() + b.plot() # needs sage.plot Graphics object consisting of 250 graphics primitives """ @@ -408,7 +411,7 @@ cdef class FastFourierTransform_complex(FastFourierTransform_base): sage: for i in range(1, 60): b[i]=1 sage: a.forward_transform() sage: a.inverse_transform() - sage: (a.plot()+b.plot()) + sage: a.plot() + b.plot() # needs sage.plot Graphics object consisting of 250 graphics primitives sage: abs(sum([CDF(a[i])-CDF(b[i]) for i in range(125)])) < 2**-16 True @@ -421,7 +424,7 @@ cdef class FastFourierTransform_complex(FastFourierTransform_base): sage: for i in range(1, 60): b[i]=1 sage: a.forward_transform() sage: a.inverse_transform() - sage: (a.plot()+b.plot()) + sage: a.plot() + b.plot() # needs sage.plot Graphics object consisting of 256 graphics primitives """ @@ -459,7 +462,7 @@ cdef class FastFourierTransform_complex(FastFourierTransform_base): sage: for i in range(1, 60): b[i]=1 sage: a.forward_transform() sage: a.backward_transform() - sage: (a.plot() + b.plot()).show(ymin=0) # long time (2s on sage.math, 2011) + sage: (a.plot() + b.plot()).show(ymin=0) # long time (2s on sage.math, 2011), needs sage.plot sage: abs(sum([CDF(a[i])/125-CDF(b[i]) for i in range(125)])) < 2**-16 True @@ -471,7 +474,7 @@ cdef class FastFourierTransform_complex(FastFourierTransform_base): sage: for i in range(1, 60): b[i]=1 sage: a.forward_transform() sage: a.backward_transform() - sage: (a.plot() + b.plot()).show(ymin=0) + sage: (a.plot() + b.plot()).show(ymin=0) # needs sage.plot """ cdef gsl_fft_complex_wavetable * wt cdef gsl_fft_complex_workspace * mem diff --git a/src/sage/calculus/wester.py b/src/sage/calculus/wester.py index 766ce586eb2..e33409a49ac 100644 --- a/src/sage/calculus/wester.py +++ b/src/sage/calculus/wester.py @@ -1,3 +1,4 @@ +# sage.doctest: needs sage.symbolic r""" Further examples from Wester's paper diff --git a/src/sage/structure/nonexact.py b/src/sage/structure/nonexact.py index f5894af43ac..fa5157b3a65 100644 --- a/src/sage/structure/nonexact.py +++ b/src/sage/structure/nonexact.py @@ -9,7 +9,7 @@ sage: R. = PowerSeriesRing(QQ) sage: R.default_prec() 20 - sage: cos(x) # optional - sage.symbolic + sage: cos(x) # needs sage.symbolic 1 - 1/2*x^2 + 1/24*x^4 - 1/720*x^6 + 1/40320*x^8 - 1/3628800*x^10 + 1/479001600*x^12 - 1/87178291200*x^14 + 1/20922789888000*x^16 - 1/6402373705728000*x^18 + O(x^20) @@ -19,7 +19,7 @@ sage: R. = PowerSeriesRing(QQ, default_prec=10) sage: R.default_prec() 10 - sage: cos(x) + sage: cos(x) # needs sage.symbolic 1 - 1/2*x^2 + 1/24*x^4 - 1/720*x^6 + 1/40320*x^8 + O(x^10) .. NOTE::