# Ticket #12455: trac_12455_newstyle_airy2.patch

File trac_12455_newstyle_airy2.patch, 41.4 KB (added by eviatarbach, 6 years ago)
• ## doc/en/reference/functions/index.rst

# HG changeset patch
# User Eviatar Bach <eviatarbach@gmail.com>
# Date 1371600084 25200
# Branch airy
# Node ID 7ae7060f2c1b674e08d62973eb44f939327e86f6
# Parent  a52dffb6d915258df25e61dd167fe46dd72e0f4e
Various fixes, including fixing an incorrect identity

diff --git a/doc/en/reference/functions/index.rst b/doc/en/reference/functions/index.rst
 a sage/functions/orthogonal_polys sage/functions/other sage/functions/special sage/functions/airy sage/functions/exp_integral sage/functions/wigner sage/functions/generalized
• ## sage/functions/airy.py

diff --git a/sage/functions/airy.py b/sage/functions/airy.py
 a r""" Airy Functions This module implements Airy functions and their generalized derivatives. It supports symbolic functionality through Maxima and numeric evaluation through mpmath and Maxima. Airy functions are solutions to the differential equation f''(z) +f(z)x=0. AUTHORS: - Oscar Gerardo Lazo Arjona (2010): initial version - Douglas McNeil (2012): rewrite EXAMPLES: Verify that the Airy functions are solutions to the differential equation:: sage: diff(airy_ai(x), x, 2) - x * airy_ai(x) 0 sage: diff(airy_bi(x), x, 2) - x * airy_bi(x) 0 """ #***************************************************************************** #       Copyright (C) 2010 Oscar Gerardo Lazo Arjona algebraicamente@gmail.com #      Copyright (C) 2010 Oscar Gerardo Lazo Arjona #      Copyright (C) 2012 Douglas McNeil # #  Distributed under the terms of the GNU General Public License (GPL) # #    but WITHOUT ANY WARRANTY; without even the implied warranty of #    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU #    General Public License for more details. # #  The full text of the GPL is available at: # #  as published by the Free Software Foundation; either version 2 of #  the License, or (at your option) any later version. #                  http://www.gnu.org/licenses/ #***************************************************************************** from sage.symbolic.expression import Expression from sage.symbolic.function import is_inexact from sage.structure.coerce import parent as sage_structure_coerce_parent from sage.functions.other import gamma from sage.rings.integer_ring import ZZ from sage.rings.real_double import RDF from sage.rings.rational import Rational as R from sage.functions.special import meval from sage.calculus.var import var from sage.calculus.functional import derivative # from sage.rings.real_mpfr import RR from sage.rings.rational import Rational as R class FunctionAiryAiGeneral(BuiltinFunction): def __init__(self): The generalized derivative of the Airy Ai function INPUT: - alpha.- Return the \alpha-th order fractional derivative with respect to z. For \alpha = n = 1,2,3,\ldots this gives the derivative \operatorname{Ai}^{(n)}(z), and for \alpha = -n = -1,-2,-3,\ldots this gives the n-fold iterated integral - alpha -- Return the \alpha-th order fractional derivative with respect to z. For \alpha = n = 1,2,3,\ldots this gives the derivative \operatorname{Ai}^{(n)}(z), and for \alpha = -n = -1,-2,-3,\ldots this gives the n-fold iterated integral. .. math :: f_0(z) = \operatorname{Ai}(z) f_n(z) = \int_0^z f_{n-1}(t) dt. - x.- The argument of the function. f_n(z) = \int_0^z f_{n-1}(t) dt - hold.- Whether or not to stop from returning higher derivatives in terms of \operatorname{Ai}(x) and \operatorname{Ai}'(x). - x -- The argument of the function - hold_derivative -- Whether or not to stop from returning higher derivatives in terms of \operatorname{Ai}(x) and \operatorname{Ai}'(x) EXAMPLES:: sage: x,n=var('x n') sage: airy_ai(-2,x) sage: x, n = var('x n') sage: airy_ai(-2, x) airy_ai(-2, x) sage: derivative(airy_ai(-2,x),x) sage: derivative(airy_ai(-2, x), x) airy_ai(-1, x) sage: airy_ai(n,x) sage: airy_ai(n, x) airy_ai(n, x) sage: derivative(airy_ai(n,x),x) sage: derivative(airy_ai(n, x), x) airy_ai(n + 1, x) sage: airy_ai(2,x,True) sage: airy_ai(2, x, hold_derivative=True) airy_ai(2, x) sage: derivative(airy_ai(2,x,hold_derivative=True),x) sage: derivative(airy_ai(2, x, hold_derivative=True), x) airy_ai(3, x) """ BuiltinFunction.__init__(self, "airy_ai", nargs=2, latex_name=r"\operatorname{Ai}") latex_name=r"\operatorname{Ai}") def _derivative_(self, alpha, *args, **kwds): """ EXAMPLES:: sage: x,n=var('x n') sage: derivative(airy_ai(n,x),x) # indirect doctest sage: x, n = var('x n') sage: derivative(airy_ai(n, x), x) # indirect doctest airy_ai(n + 1, x) """ x=args[0] return airy_ai_general(alpha+1,x) x = args[0] return airy_ai_general(alpha + 1, x) def _eval_(self, alpha, *args): """ EXAMPLES:: sage: x,n=var('x n') sage: airy_ai(-2,1.0) # indirect doctest sage: x, n = var('x n') sage: airy_ai(-2, 1.0) # indirect doctest 0.136645379421096 sage: airy_ai(n,1.0) # indirect doctest sage: airy_ai(n, 1.0) # indirect doctest airy_ai(n, 1.00000000000000) """ x = args[0] x=args[0] if not isinstance(x,Expression) and not isinstance(alpha,Expression): if not isinstance(x, Expression) and not isinstance(alpha, Expression): if is_inexact(x): return self._evalf_(alpha,x) return self._evalf_(alpha, x) else: return None """ EXAMPLES:: sage: airy_ai(-2,1.0) # indirect doctest 0.136645379421096 sage: airy_ai(-2, 1.0) # indirect doctest 0.136645379421096 """ algorithm = kwargs.get('algorithm', None) or 'mpmath' parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x) prec = parent.prec() if hasattr(parent, 'prec') else 53 parent = sage_structure_coerce_parent(x) if algorithm == 'mpmath': import mpmath from sage.libs.mpmath import utils as mpmath_utils return mpmath_utils.call(mpmath.airyai, x, derivative=alpha, parent=parent) return mpmath_utils.call(mpmath.airyai, x, derivative=alpha, parent=parent) elif algorithm == 'maxima': raise NotImplementedError("general case not available in maxima") raise NotImplementedError("general case not available in Maxima") else: raise ValueError("unknown algorithm") class FunctionAiryAiSimple(BuiltinFunction): def __init__(self): """ The class for the Airy Ai function. The class for the Airy Ai function EXAMPLES:: sage: from sage.functions.airy import FunctionAiryAiSimple sage: airy_ai_simple= FunctionAiryAiSimple() sage: f=airy_ai_simple(x); f sage: airy_ai_simple = FunctionAiryAiSimple() sage: f = airy_ai_simple(x); f airy_ai(x) """ BuiltinFunction.__init__(self, "airy_ai", latex_name=r'\operatorname{Ai}', conversions=dict(mathematica='AiryAi', maxima='airy_ai')) BuiltinFunction.__init__(self, "airy_ai", latex_name=r'\operatorname{Ai}', conversions=dict(mathematica='AiryAi', maxima='airy_ai')) def _derivative_(self, x, diff_param=None): """ EXAMPLES:: sage: derivative(airy_ai(x),x) # indirect doctest sage: derivative(airy_ai(x), x) # indirect doctest airy_ai_prime(x) """ return airy_ai_prime(x) def _eval_(self, x): """ EXAMPLES:: sage: airy_ai(0) # indirect doctest sage: airy_ai(0)  # indirect doctest 1/3*3^(1/3)/gamma(2/3) sage: airy_ai(0.0) # indirect doctest sage: airy_ai(0.0)  # indirect doctest 0.355028053887817 sage: airy_ai(I) # indirect doctest sage: airy_ai(I)  # indirect doctest airy_ai(I) sage: airy_ai(1.0*I) # indirect doctest sage: airy_ai(1.0 * I) # indirect doctest 0.331493305432141 - 0.317449858968444*I """ if not isinstance(x,Expression): if not isinstance(x, Expression): if is_inexact(x): return self._evalf_(x) elif x==0: r=ZZ(2)/3 return 1/(3**(r)*gamma(r)) elif x == 0: r = ZZ(2) / 3 return 1 / (3 ** (r) * gamma(r)) else: return None def _evalf_(self, x, **kwargs): """ EXAMPLES:: sage: airy_ai(0.0) # indirect doctest sage: airy_ai(0.0)  # indirect doctest 0.355028053887817 sage: airy_ai(1.0*I) # indirect doctest sage: airy_ai(1.0 * I) # indirect doctest 0.331493305432141 - 0.317449858968444*I We can use several methods for numerical evaluation:: 0.00659113935746 sage: airy_ai(3).n(algorithm='mpmath') 0.00659113935746072 sage: airy_ai(3).n(algorithm='mpmath',prec=100) sage: airy_ai(3).n(algorithm='mpmath', prec=100) 0.0065911393574607191442574484080 TESTS:: sage: airy_ai(3).n(algorithm='maxima', prec=70) Traceback (most recent call last): ... ValueError: for the maxima algorithm the precision must be 53 ValueError: for the Maxima algorithm the precision must be 53 """ algorithm = kwargs.get('algorithm', None) or 'mpmath' parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x) prec = parent.prec() if hasattr(parent, 'prec') else 53 parent = sage_structure_coerce_parent(x) if algorithm == 'mpmath': import mpmath from sage.libs.mpmath import utils as mpmath_utils return mpmath_utils.call(mpmath.airyai, x, parent=parent) elif algorithm == 'maxima': if prec != 53: raise ValueError, "for the maxima algorithm the precision must be 53" return RDF(meval("airy_ai(%s)" % RDF(x))) raise ValueError("for the Maxima algorithm the precision must " "be 53") return RDF(meval("airy_ai({})".format(RDF(x)))) else: raise ValueError("unknown algorithm") class FunctionAiryAiPrime(BuiltinFunction): def __init__(self): """ EXAMPLES:: sage: x,n=var('x n') sage: x, n = var('x n') sage: airy_ai_prime(x) airy_ai_prime(x) sage: airy_ai_prime(0) -1/3*3^(2/3)/gamma(1/3) """ BuiltinFunction.__init__(self, "airy_ai_prime", latex_name=r"\operatorname{Ai}'", conversions=dict(mathematica='AiryAiPrime', maxima='airy_dai')) latex_name=r"\operatorname{Ai}'", conversions=dict(mathematica='AiryAiPrime', maxima='airy_dai')) def _derivative_(self, x, diff_param=None): """ EXAMPLES:: sage: derivative(airy_ai_prime(x),x) # indirect doctest sage: derivative(airy_ai_prime(x), x) # indirect doctest x*airy_ai(x) """ return x*airy_ai_simple(x) return x * airy_ai_simple(x) def _eval_(self, x): """ EXAMPLES:: sage: airy_ai(1,0) # indirect doctest sage: airy_ai(1, 0) # indirect doctest -1/3*3^(2/3)/gamma(1/3) sage: airy_ai(1,0.0) # indirect doctest sage: airy_ai(1, 0.0) # indirect doctest -0.258819403792807 """ if not isinstance(x,Expression): if not isinstance(x, Expression): if is_inexact(x): return self._evalf_(x) elif x==0: r=ZZ(1)/3 return -1/(3**(r)*gamma(r)) elif x == 0: r = ZZ(1) / 3 return -1 / (3 ** (r) * gamma(r)) else: return None def _evalf_(self, x, **kwargs): """ EXAMPLES:: sage: airy_ai(1,0.0) # indirect doctest sage: airy_ai(1, 0.0) # indirect doctest -0.258819403792807 We can use several methods for numerical evaluation:: sage: airy_ai(1,4).n(algorithm='maxima') sage: airy_ai(1, 4).n(algorithm='maxima') -0.0019586409502 sage: airy_ai(1,4).n(algorithm='mpmath') sage: airy_ai(1, 4).n(algorithm='mpmath') -0.00195864095020418 sage: airy_ai(1,4).n(algorithm='mpmath', prec=100) sage: airy_ai(1, 4).n(algorithm='mpmath', prec=100) -0.0019586409502041789001381409184 TESTS:: sage: airy_ai(1, 4).n(algorithm='maxima', prec=70) Traceback (most recent call last): ... ValueError: for the maxima algorithm the precision must be 53 ValueError: for the Maxima algorithm the precision must be 53 """ algorithm = kwargs.get('algorithm', None) or 'mpmath' parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x) prec = parent.prec() if hasattr(parent, 'prec') else 53 parent = sage_structure_coerce_parent(x) if algorithm == 'mpmath': import mpmath from sage.libs.mpmath import utils as mpmath_utils return mpmath_utils.call(mpmath.airyai, x, derivative=1, parent=parent) return mpmath_utils.call(mpmath.airyai, x, derivative=1, parent=parent) elif algorithm == 'maxima': if prec != 53: raise ValueError, "for the maxima algorithm the precision must be 53" return RDF(meval("airy_dai(%s)" % RDF(x))) raise ValueError("for the Maxima algorithm the precision must " "be 53") return RDF(meval("airy_dai({})".format(RDF(x)))) else: raise ValueError("unknown algorithm") airy_ai_general = FunctionAiryAiGeneral() airy_ai_simple = FunctionAiryAiSimple() airy_ai_prime = FunctionAiryAiPrime() airy_ai_general=FunctionAiryAiGeneral() airy_ai_simple= FunctionAiryAiSimple() airy_ai_prime=  FunctionAiryAiPrime() def airy_ai(alpha,x=None, hold_derivative=False, *args, **kwds): r""" def airy_ai(alpha, x=None, hold_derivative=False, *args, **kwds): r""" The Airy Ai function \operatorname{Ai}(x) is one of the two linearly independent solutions to the Airy differental equation linearly independent solutions to the Airy differential equation f''(z) +f(z)x=0, defined by the initial conditions: .. math :: \operatorname{Ai}(0)=\frac{1}{2^{2/3} \Gamma(\frac{2}{3})}, \operatorname{Ai}(0)=\frac{1}{2^{2/3} \Gamma\left(\frac{2}{3}\right)}, \operatorname{Ai}'(0)=-\frac{1}{2^{1/3} \Gamma(\frac{1}{3})}. \operatorname{Ai}'(0)=-\frac{1}{2^{1/3}\Gamma\left(\frac{1}{3}\right)}. Another way to define the Airy Ai function is: \cos\left(\frac{1}{3}t^3+xt\right) dt. INPUT: - alpha.- Return the \alpha-th order fractional derivative with respect to z. For \alpha = n = 1,2,3,\ldots this gives the derivative \operatorname{Ai}^{(n)}(z), and for \alpha = -n = -1,-2,-3,\ldots this gives the n-fold iterated integral - alpha -- Return the \alpha-th order fractional derivative with respect to z. For \alpha = n = 1,2,3,\ldots this gives the derivative \operatorname{Ai}^{(n)}(z), and for \alpha = -n = -1,-2,-3,\ldots this gives the n-fold iterated integral. .. math :: f_0(z) = \operatorname{Ai}(z) f_n(z) = \int_0^z f_{n-1}(t) dt. - x.- The argument of the function. f_n(z) = \int_0^z f_{n-1}(t) dt - hold_derivative.- Whether or not to stop from returning higher derivatives in terms of \operatorname{Ai}(x) and \operatorname{Ai}'(x). - x -- The argument of the function - hold_derivative -- Whether or not to stop from returning higher derivatives in terms of \operatorname{Ai}(x) and \operatorname{Ai}'(x) EXAMPLES:: sage: n,x=var('n x') sage: n, x = var('n x') sage: airy_ai(x) airy_ai(x) It can return derivatives or integrals:: sage: airy_ai(1,x) sage: airy_ai(1, x) airy_ai_prime(x) sage: airy_ai(2,x) sage: airy_ai(2, x) x*airy_ai(x) sage: airy_ai(2,x,True) sage: airy_ai(2, x, hold_derivative=True) airy_ai(2, x) sage: airy_ai(-2,x) sage: airy_ai(-2, x) airy_ai(-2, x) sage: airy_ai(n, x) airy_ai(n, x) It can be evaluated symbolically or numerically for real or complex values:: It can be evaluated symbolically or numerically for real or complex values:: sage: airy_ai(0) 1/3*3^(1/3)/gamma(2/3) sage: airy_ai(1.0*I) 0.331493305432141 - 0.317449858968444*I The functions can be evaluated numerically using either mpmath (the default) or maxima, where mpmath can compute the values to arbitrary precision:: The functions can be evaluated numerically using either mpmath (the default) or Maxima, where mpmath can compute the values to arbitrary precision:: sage: airy_ai(2).n(prec=100) 0.034924130423274379135322080792 sage: airy_ai(2).n(algorithm='mpmath',prec=100) sage: airy_ai(2).n(algorithm='mpmath', prec=100) 0.034924130423274379135322080792 sage: airy_ai(2).n(algorithm='maxima') 0.0349241304233 And the derivatives can be evaluated:: sage: airy_ai(1,0) sage: airy_ai(1, 0) -1/3*3^(2/3)/gamma(1/3) sage: airy_ai(1,0.0) sage: airy_ai(1, 0.0) -0.258819403792807 Plots:: sage: plot(airy_ai(x),(x,-10,5))+plot(airy_ai_prime(x),(x,-10,5),color='red') sage: (plot(airy_ai(x), (x, -10, 5)) +\ plot(airy_ai_prime(x), (x, -10, 5), color='red')) **References** - Abramowitz, Milton; Stegun, Irene A., eds. (1965), "Chapter 10" - http://en.wikipedia.org/wiki/Airy_function """ # We catch the case with no alpha if x is None: x = alpha return airy_ai_simple(x, **kwds) #We catch the case with no alpha if x==None: x=alpha # We raise an error if there are too many arguments if len(args) > 0: raise TypeError(("symbolic function airy_ai takes at most 3 arguments " "({} given)").format(len(args) + 3)) # We take care of all other cases. if not alpha in ZZ and not isinstance(alpha, Expression): raise ValueError("alpha must be an integer") if hold_derivative: return airy_ai_general(alpha, x, **kwds) elif alpha == 0: return airy_ai_simple(x, **kwds) #We raise an error if there are too many arguments if len(args) > 0: raise TypeError("Symbolic function airy_ai takes at most 3 arguments (%s given)" % (len(args)+3)) #We take care of all other cases. if hold_derivative: return airy_ai_general(alpha,x,**kwds) elif alpha==0: return airy_ai_simple(x, **kwds) elif alpha==1: elif alpha == 1: return airy_ai_prime(x, **kwds) elif alpha>1: #we use a different variable here because if x is a #particular value, we would be differentiating a constant #which would return 0. What we want is the value of #the derivative at the value and not the derivative of #a particular value of the function. v=var('v') return derivative(airy_ai_simple(v,**kwds),v,alpha).subs(v=x) elif alpha > 1: # We use a different variable here because if x is a # particular value, we would be differentiating a constant # which would return 0. What we want is the value of # the derivative at the value and not the derivative of # a particular value of the function. v = var('v') return derivative(airy_ai_simple(v, **kwds), v, alpha).subs(v=x) else: return airy_ai_general(alpha,x,**kwds) return airy_ai_general(alpha, x, **kwds) ######################################################################## ######################################################################## class FunctionAiryBiGeneral(BuiltinFunction): def __init__(self): r""" The generalized derivative of the Airy Bi function INPUT: - alpha.- Return the \alpha-th order fractional derivative with respect to z. For \alpha = n = 1,2,3,\ldots this gives the derivative \operatorname{Bi}^{(n)}(z), and for \alpha = -n = -1,-2,-3,\ldots this gives the n-fold iterated integral - alpha -- Return the \alpha-th order fractional derivative with respect to z. For \alpha = n = 1,2,3,\ldots this gives the derivative \operatorname{Bi}^{(n)}(z), and for \alpha = -n = -1,-2,-3,\ldots this gives the n-fold iterated integral. .. math :: f_0(z) = \operatorname{Bi}(z) f_n(z) = \int_0^z f_{n-1}(t) dt. - x.- The argument of the function. f_n(z) = \int_0^z f_{n-1}(t) dt - hold.- Whether or not to stop from returning higher derivatives in terms of \operatorname{Bi}(x) and \operatorname{Bi}'(x). - x -- The argument of the function - hold_derivative -- Whether or not to stop from returning higher derivatives in terms of \operatorname{Bi}(x) and \operatorname{Bi}'(x) EXAMPLES:: sage: x,n=var('x n') sage: airy_bi(-2,x) sage: x, n = var('x n') sage: airy_bi(-2, x) airy_bi(-2, x) sage: derivative(airy_bi(-2,x),x) sage: derivative(airy_bi(-2, x), x) airy_bi(-1, x) sage: airy_bi(n,x) sage: airy_bi(n, x) airy_bi(n, x) sage: derivative(airy_bi(n,x),x) sage: derivative(airy_bi(n, x), x) airy_bi(n + 1, x) sage: airy_bi(2,x,True) sage: airy_bi(2, x, hold_derivative=True) airy_bi(2, x) sage: derivative(airy_bi(2,x,hold_derivative=True),x) sage: derivative(airy_bi(2, x, hold_derivative=True), x) airy_bi(3, x) """ BuiltinFunction.__init__(self, "airy_bi", nargs=2, latex_name=r"\operatorname{Bi}") latex_name=r"\operatorname{Bi}") def _derivative_(self, alpha, *args, **kwds): """ EXAMPLES:: sage: x,n=var('x n') sage: derivative(airy_bi(n,x),x) # indirect doctest sage: x, n = var('x n') sage: derivative(airy_bi(n, x), x) # indirect doctest airy_bi(n + 1, x) """ x=args[0] return airy_bi_general(alpha+1,x) x = args[0] return airy_bi_general(alpha + 1, x) def _eval_(self, alpha, *args): """ EXAMPLES:: sage: x,n=var('x n') sage: airy_bi(-2,1.0) # indirect doctest sage: x, n = var('x n') sage: airy_bi(-2, 1.0) # indirect doctest 0.388621540699059 sage: airy_bi(n,1.0) # indirect doctest sage: airy_bi(n, 1.0) # indirect doctest airy_bi(n, 1.00000000000000) """ x=args[0] if not isinstance(x,Expression) and not isinstance(alpha,Expression): x = args[0] if not isinstance(x, Expression) and not isinstance(alpha, Expression): if is_inexact(x): return self._evalf_(alpha,x) return self._evalf_(alpha, x) else: return None """ EXAMPLES:: sage: airy_bi(-2,1.0) # indirect doctest sage: airy_bi(-2, 1.0) # indirect doctest 0.388621540699059 """ algorithm = kwargs.get('algorithm', None) or 'mpmath' parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x) prec = parent.prec() if hasattr(parent, 'prec') else 53 parent = sage_structure_coerce_parent(x) if algorithm == 'mpmath': import mpmath from sage.libs.mpmath import utils as mpmath_utils return mpmath_utils.call(mpmath.airybi, x, derivative=alpha, parent=parent) return mpmath_utils.call(mpmath.airybi, x, derivative=alpha, parent=parent) elif algorithm == 'maxima': raise NotImplementedError("general case not available in maxima") raise NotImplementedError("general case not available in Maxima") else: raise ValueError("unknown algorithm") class FunctionAiryBiSimple(BuiltinFunction): def __init__(self): """ The class for the Airy Bi function. The class for the Airy Bi function EXAMPLES:: sage: from sage.functions.airy import FunctionAiryBiSimple sage: airy_bi_simple= FunctionAiryBiSimple() sage: f=airy_bi_simple(x); f sage: f = airy_bi_simple(x); f airy_bi(x) """ BuiltinFunction.__init__(self, "airy_bi", latex_name=r'\operatorname{Bi}', conversions=dict(mathematica='AiryBi', maxima='airy_bi')) BuiltinFunction.__init__(self, "airy_bi", latex_name=r'\operatorname{Bi}', conversions=dict(mathematica='AiryBi', maxima='airy_bi')) def _derivative_(self, x, diff_param=None): """ EXAMPLES:: sage: derivative(airy_bi(x),x) # indirect doctest sage: derivative(airy_bi(x), x) # indirect doctest airy_bi_prime(x) """ return airy_bi_prime(x) def _eval_(self, x): """ EXAMPLES:: sage: airy_bi(0) # indirect doctest 1/3*3^(5/6)/gamma(1/3) sage: airy_bi(0.0) # indirect doctest sage: airy_bi(0)  # indirect doctest 1/3*3^(5/6)/gamma(2/3) sage: airy_bi(0.0)  # indirect doctest 0.614926627446001 sage: airy_bi(I) # indirect doctest sage: airy_bi(0).n() == airy_bi(0.0)  # indirect doctest True sage: airy_bi(I)  # indirect doctest airy_bi(I) sage: airy_bi(1.0*I) # indirect doctest sage: airy_bi(1.0 * I) # indirect doctest 0.648858208330395 + 0.344958634768048*I """ if not isinstance(x,Expression): if not isinstance(x, Expression): if is_inexact(x): return self._evalf_(x) elif x==0: one_sixth = ZZ(1)/6 return 1/(3**(one_sixth)*gamma(2*one_sixth)) elif x == 0: one_sixth = ZZ(1) / 6 return 1 / (3 ** (one_sixth) * gamma(4 * one_sixth)) else: return None def _evalf_(self, x, **kwargs): """ EXAMPLES:: sage: airy_bi(0.0) # indirect doctest sage: airy_bi(0.0)  # indirect doctest 0.614926627446001 sage: airy_bi(1.0*I) # indirect doctest sage: airy_bi(1.0 * I) # indirect doctest 0.648858208330395 + 0.344958634768048*I We can use several methods for numerical evaluation:: sage: airy_bi(3).n(algorithm='maxima', prec=100) Traceback (most recent call last): ... ValueError: for the maxima algorithm the precision must be 53 ValueError: for the Maxima algorithm the precision must be 53 """ algorithm = kwargs.get('algorithm', None) or 'mpmath' parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x) prec = parent.prec() if hasattr(parent, 'prec') else 53 parent = sage_structure_coerce_parent(x) if algorithm == 'mpmath': import mpmath from sage.libs.mpmath import utils as mpmath_utils return mpmath_utils.call(mpmath.airybi, x, parent=parent) elif algorithm == 'maxima': if prec != 53: raise ValueError, "for the maxima algorithm the precision must be 53" return RDF(meval("airy_bi(%s)" % RDF(x))) raise ValueError("for the Maxima algorithm the precision must " "be 53") return RDF(meval("airy_bi({})".format(RDF(x)))) else: raise ValueError("unknown algorithm") class FunctionAiryBiPrime(BuiltinFunction): def __init__(self): """ EXAMPLES:: sage: x,n=var('x n') sage: airy_bi_prime(x) # indirect doctest sage: x, n = var('x n') sage: airy_bi_prime(x)  # indirect doctest airy_bi_prime(x) sage: airy_bi_prime(0) # indirect doctest sage: airy_bi_prime(0)  # indirect doctest 3^(1/6)/gamma(1/3) """ BuiltinFunction.__init__(self, "airy_bi_prime", latex_name=r"\operatorname{Bi}'", conversions=dict(mathematica='AiryBiPrime', maxima='airy_dbi')) latex_name=r"\operatorname{Bi}'", conversions=dict(mathematica='AiryBiPrime', maxima='airy_dbi')) def _derivative_(self, x, diff_param=None): """ EXAMPLES:: sage: derivative(airy_bi_prime(x),x) # indirect doctest sage: derivative(airy_bi_prime(x), x) # indirect doctest x*airy_bi(x) """ return x*airy_bi_simple(x) return x * airy_bi_simple(x) def _eval_(self, x): """ EXAMPLES:: sage: airy_bi(1,0) # indirect doctest sage: airy_bi(1, 0) # indirect doctest 3^(1/6)/gamma(1/3) sage: airy_bi(1,0.0) # indirect doctest sage: airy_bi(1, 0.0) # indirect doctest 0.448288357353826 """ if not isinstance(x,Expression): if not isinstance(x, Expression): if is_inexact(x): return self._evalf_(x) elif x==0: one_sixth = ZZ(1)/6 return 3**(one_sixth)/gamma(2*one_sixth) elif x == 0: one_sixth = ZZ(1) / 6 return 3 ** (one_sixth) / gamma(2 * one_sixth) else: return None def _evalf_(self, x, **kwargs): """ EXAMPLES:: sage: airy_bi(1,0.0) # indirect doctest sage: airy_bi(1, 0.0) # indirect doctest 0.448288357353826 We can use several methods for numerical evaluation:: sage: airy_bi(1,4).n(algorithm='maxima') sage: airy_bi(1, 4).n(algorithm='maxima') 161.926683505 sage: airy_bi(1,4).n(algorithm='mpmath') sage: airy_bi(1, 4).n(algorithm='mpmath') 161.926683504613 sage: airy_bi(1,4).n(algorithm='mpmath', prec=100) sage: airy_bi(1, 4).n(algorithm='mpmath', prec=100) 161.92668350461340184309492429 TESTS:: sage: airy_bi(1,4).n(algorithm='maxima', prec=70) sage: airy_bi(1, 4).n(algorithm='maxima', prec=70) Traceback (most recent call last): ... ValueError: for the maxima algorithm the precision must be 53 ValueError: for the Maxima algorithm the precision must be 53 """ algorithm = kwargs.get('algorithm', None) or 'mpmath' parent = kwargs.get('parent', None) or sage_structure_coerce_parent(x) prec = parent.prec() if hasattr(parent, 'prec') else 53 parent = sage_structure_coerce_parent(x) if algorithm == 'mpmath': import mpmath from sage.libs.mpmath import utils as mpmath_utils return mpmath_utils.call(mpmath.airybi, x, derivative=1, parent=parent) return mpmath_utils.call(mpmath.airybi, x, derivative=1, parent=parent) elif algorithm == 'maxima': if prec != 53: raise ValueError, "for the maxima algorithm the precision must be 53" return RDF(meval("airy_dbi(%s)" % RDF(x))) raise ValueError("for the Maxima algorithm the precision must " "be 53") return RDF(meval("airy_dbi({})".format(RDF(x)))) else: raise ValueError("unknown algorithm") airy_bi_general=FunctionAiryBiGeneral() airy_bi_simple= FunctionAiryBiSimple() airy_bi_prime= FunctionAiryBiPrime() airy_bi_general = FunctionAiryBiGeneral() airy_bi_simple = FunctionAiryBiSimple() airy_bi_prime = FunctionAiryBiPrime() def airy_bi(alpha,x=None, hold_derivative=False, *args, **kwds): def airy_bi(alpha, x=None, hold_derivative=False, *args, **kwds): r""" The Airy Bi function \operatorname{Bi}(x) is one of the two linearly independent solutions to the Airy differental equation linearly independent solutions to the Airy differential equation f''(z) +f(z)x=0, defined by the initial conditions: .. math :: \operatorname{Bi}(0)=\frac{1}{3^{1/6} \Gamma(\frac{2}{3})}, \operatorname{Bi}(0)=\frac{1}{3^{1/6} \Gamma\left(\frac{2}{3}\right)}, \operatorname{Bi}'(0)=\frac{3^{1/6}}{ \Gamma(\frac{1}{3})}. \operatorname{Bi}'(0)=\frac{3^{1/6}}{ \Gamma\left(\frac{1}{3}\right)}. Another way to define the Airy Bi function is: +\sin\left(xt + \frac{1}{3}t^3\right) \right ] dt. INPUT: - alpha.- Return the \alpha-th order fractional derivative with respect to z. For \alpha = n = 1,2,3,\ldots this gives the derivative \operatorname{Bi}^{(n)}(z), and for \alpha = -n = -1,-2,-3,\ldots this gives the n-fold iterated integral - alpha -- Return the \alpha-th order fractional derivative with respect to z. For \alpha = n = 1,2,3,\ldots this gives the derivative \operatorname{Bi}^{(n)}(z), and for \alpha = -n = -1,-2,-3,\ldots this gives the n-fold iterated integral. .. math :: f_0(z) = \operatorname{Bi}(z) f_n(z) = \int_0^z f_{n-1}(t) dt. - x.- The argument of the function. f_n(z) = \int_0^z f_{n-1}(t) dt - hold_derivative.- Whether or not to stop from returning higher derivatives in terms of \operatorname{Bi}(x) and \operatorname{Bi}'(x). - x -- The argument of the function - hold_derivative -- Whether or not to stop from returning higher derivatives in terms of \operatorname{Bi}(x) and \operatorname{Bi}'(x) EXAMPLES:: sage: n,x=var('n x') sage: n, x = var('n x') sage: airy_bi(x) airy_bi(x) It can return derivatives or integrals:: sage: airy_bi(1,x) sage: airy_bi(1, x) airy_bi_prime(x) sage: airy_bi(2,x) sage: airy_bi(2, x) x*airy_bi(x) sage: airy_bi(2,x,True) sage: airy_bi(2, x, hold_derivative=True) airy_bi(2, x) sage: airy_bi(-2,x) sage: airy_bi(-2, x) airy_bi(-2, x) sage: airy_bi(n, x) airy_bi(n, x) It can be evaluated symbolically or numerically for real or complex values:: It can be evaluated symbolically or numerically for real or complex values:: sage: airy_bi(0) 1/3*3^(5/6)/gamma(1/3) 1/3*3^(5/6)/gamma(2/3) sage: airy_bi(0.0) 0.614926627446001 sage: airy_bi(I) sage: airy_bi(1.0*I) 0.648858208330395 + 0.344958634768048*I The functions can be evaluated numerically using either mpmath (the default) or maxima, where mpmath can compute the values to arbitrary precision:: The functions can be evaluated numerically using either mpmath (the default) or Maxima, where mpmath can compute the values to arbitrary precision:: sage: airy_bi(2).n(prec=100) 3.2980949999782147102806044252 sage: airy_bi(2).n(algorithm='mpmath',prec=100) sage: airy_bi(2).n(algorithm='mpmath', prec=100) 3.2980949999782147102806044252 sage: airy_bi(2).n(algorithm='maxima') 3.29809499998 And the derivatives can be evaluated:: sage: airy_bi(1,0) sage: airy_bi(1, 0) 3^(1/6)/gamma(1/3) sage: airy_bi(1,0.0) sage: airy_bi(1, 0.0) 0.448288357353826 Plots:: sage: plot(airy_bi(x),(x,-10,5))+plot(airy_bi_prime(x),(x,-10,5),color='red') sage: (plot(airy_bi(x), (x, -10, 5)) +\ plot(airy_bi_prime(x), (x, -10, 5), color='red')) **References** - Abramowitz, Milton; Stegun, Irene A., eds. (1965), "Chapter 10" - http://en.wikipedia.org/wiki/Airy_function """ # We catch the case with no alpha if x is None: x = alpha return airy_bi_simple(x, **kwds) #We catch the case with no alpha if x==None: x=alpha # We raise an error if there are too many arguments if len(args) > 0: raise TypeError(("symbolic function airy_ai takes at most 3 arguments " "({} given)").format(len(args) + 3)) # We take care of all other cases. if hold_derivative: return airy_bi_general(alpha, x, **kwds) elif alpha == 0: return airy_bi_simple(x, **kwds) #We raise an error if there are too many arguments if len(args) > 0: raise TypeError("Symbolic function airy_ai takes at most 3 arguments (%s given)" % (len(args)+3)) #We take care of all other cases. if hold_derivative: return airy_bi_general(alpha,x,**kwds) elif alpha==0: return airy_bi_simple(x, **kwds) elif alpha==1: elif alpha == 1: return airy_bi_prime(x, **kwds) elif alpha>1: #we use a different variable here because if x is a #particular value, we would be differentiating a constant #which would return 0. What we want is the value of #the derivative at the value and not the derivative of #a particular value of the function. v=var('v') return derivative(airy_bi_simple(v,**kwds),v,alpha).subs(v=x) elif alpha > 1: # We use a different variable here because if x is a # particular value, we would be differentiating a constant # which would return 0. What we want is the value of # the derivative at the value and not the derivative of # a particular value of the function. v = var('v') return derivative(airy_bi_simple(v, **kwds), v, alpha).subs(v=x) else: return airy_bi_general(alpha,x,**kwds) return airy_bi_general(alpha, x, **kwds)