From a215d3bc25cf6a9008e9807ed71f79d663c00236 Mon Sep 17 00:00:00 2001 From: Ralf Stephan Date: Thu, 5 Mar 2015 08:19:34 +0100 Subject: [PATCH] 17678: special values of Bessel functions --- src/sage/functions/bessel.py | 115 ++++++++++++++++++---- src/sage/rings/infinity.py | 2 +- src/sage/symbolic/expression.pyx | 2 +- src/sage/symbolic/ring.pyx | 9 +- src/sage/tests/french_book/recequadiff.py | 15 +-- 5 files changed, 114 insertions(+), 29 deletions(-) diff --git a/src/sage/functions/bessel.py b/src/sage/functions/bessel.py index e0472302827..20fb2911a0f 100644 --- a/src/sage/functions/bessel.py +++ b/src/sage/functions/bessel.py @@ -97,7 +97,7 @@ sage: bessel_J(0, x) bessel_J(0, x) sage: bessel_J(0, 0) - bessel_J(0, 0) + 1 sage: bessel_J(0, x).diff(x) -1/2*bessel_J(1, x) + 1/2*bessel_J(-1, x) @@ -193,11 +193,12 @@ # remove after deprecation period from sage.calculus.calculus import maxima -from sage.functions.other import real, imag +from sage.functions.trig import sin, cos +from sage.functions.other import real, imag, sqrt from sage.misc.sage_eval import sage_eval from sage.rings.real_mpfr import RealField from sage.plot.plot import plot -from sage.rings.all import ZZ +from sage.rings.all import ZZ, QQ class Function_Bessel_J(BuiltinFunction): @@ -302,6 +303,37 @@ def __init__(self): maxima='bessel_j', sympy='besselj')) + def _eval_(self, n, x): + """ + EXAMPLES:: + + sage: n = var('n') + sage: bessel_J(0, 0) + 1 + sage: bessel_J(5/2, 0) + 0 + sage: bessel_J(-5/2, 0) + Infinity + sage: bessel_J(1/2, x) + sqrt(2)*sqrt(1/(pi*x))*sin(x) + sage: bessel_J(-1/2, x) + sqrt(2)*sqrt(1/(pi*x))*cos(x) + sage: bessel_J(n, 0) + bessel_J(n, 0) + """ + from sage.rings.infinity import unsigned_infinity + if x == 0: + if n == 0: + return ZZ(1) + elif n > 0 or n in ZZ: + return ZZ(0) + elif not isinstance(n, Expression): + return unsigned_infinity + elif n == QQ(1)/2: + return sqrt(2/pi/x) * sin(x) + elif n == QQ(-1)/2: + return sqrt(2/pi/x) * cos(x) + def _evalf_(self, n, x, parent=None, algorithm=None): """ EXAMPLES:: @@ -474,6 +506,29 @@ def __init__(self): maxima='bessel_y', sympy='bessely')) + def _eval_(self, n, x): + """ + EXAMPLES:: + + sage: n = var('n') + sage: bessel_Y(1, 0) + Infinity + sage: bessel_Y(1/2, x) + -sqrt(2)*sqrt(1/(pi*x))*cos(x) + sage: bessel_Y(-1/2, x) + sqrt(2)*sqrt(1/(pi*x))*sin(x) + """ + from sage.rings.infinity import infinity, unsigned_infinity + if x == 0: + if n == 0: + return -infinity + elif not isinstance(n, Expression): + return unsigned_infinity + elif n == QQ(1)/2: + return -sqrt(2/pi/x) * cos(x) + elif n == QQ(-1)/2: + return sqrt(2/pi/x) * sin(x) + def _evalf_(self, n, x, parent=None, algorithm=None): """ EXAMPLES:: @@ -643,26 +698,43 @@ def _eval_(self, n, x): """ EXAMPLES:: - sage: y=var('y') - sage: bessel_I(y,x) + sage: n,y = var('n,y') + sage: bessel_I(y, x) bessel_I(y, x) - sage: bessel_I(0.0, 1.0) - 1.26606587775201 + sage: bessel_I(0, 0) + 1 + sage: bessel_I(7/2, 0) + 0 + sage: bessel_I(-7/2, 0) + Infinity sage: bessel_I(1/2, 1) sqrt(2)*sinh(1)/sqrt(pi) sage: bessel_I(-1/2, pi) sqrt(2)*cosh(pi)/pi + sage: bessel_I(n, 0) + bessel_I(n, 0) """ + from sage.rings.infinity import unsigned_infinity # special identities - if n == Integer(1) / Integer(2): + if x == 0: + if n == 0: + return ZZ(1) + elif n > 0 or n in ZZ: + return ZZ(0) + elif not isinstance(n, Expression): + return unsigned_infinity + if n == QQ(1)/2: return sqrt(2 / (pi * x)) * sinh(x) - elif n == -Integer(1) / Integer(2): + elif n == -QQ(1)/2: return sqrt(2 / (pi * x)) * cosh(x) + def _evalf_(self, n, x, parent=None, algorithm=None): """ EXAMPLES:: + sage: bessel_I(0.0, 1.0) + 1.26606587775201 sage: bessel_I(1,3).n(digits=20) 3.9533702174026093965 """ @@ -738,9 +810,9 @@ class Function_Bessel_K(BuiltinFunction): -1/2*bessel_K(3, x) - 1/2*bessel_K(1, x) sage: bessel_K(1/2, x) - bessel_K(1/2, x) + sqrt(1/2)*sqrt(pi)*e^(-x)/sqrt(x) sage: bessel_K(1/2, -1) - bessel_K(1/2, -1) + -I*sqrt(1/2)*sqrt(pi)*e sage: bessel_K(1/2, 1) sqrt(1/2)*sqrt(pi)*e^(-1) @@ -816,23 +888,30 @@ def _eval_(self, n, x): """ EXAMPLES:: - sage: bessel_K(1,0) - bessel_K(1, 0) - sage: bessel_K(1.0, 0.0) - +infinity - sage: bessel_K(-1, 1).n(128) - 0.60190723019723457473754000153561733926 + sage: n = var('n') + sage: bessel_K(1, 0) + Infinity + sage: bessel_K(1/2, x) + sqrt(1/2)*sqrt(pi)*e^(-x)/sqrt(x) + sage: bessel_K(n, 0) + bessel_K(n, 0) """ + from sage.rings.infinity import unsigned_infinity # special identity - if n == Integer(1) / Integer(2) and x > 0: + if x == 0 and not isinstance(n, Expression): + return unsigned_infinity + if n == QQ(1)/2 or n == -QQ(1)/2 and x > 0: return sqrt(pi / 2) * exp(-x) * x ** (-Integer(1) / Integer(2)) + def _evalf_(self, n, x, parent=None, algorithm=None): """ EXAMPLES:: sage: bessel_K(0.0, 1.0) 0.421024438240708 + sage: bessel_K(-1, 1).n(128) + 0.60190723019723457473754000153561733926 sage: bessel_K(0, RealField(128)(1)) 0.42102443824070833333562737921260903614 """ diff --git a/src/sage/rings/infinity.py b/src/sage/rings/infinity.py index 3a9577c53be..da1d585cc2e 100644 --- a/src/sage/rings/infinity.py +++ b/src/sage/rings/infinity.py @@ -469,7 +469,7 @@ def _div_(self, other): sage: SR(infinity) / unsigned_infinity Traceback (most recent call last): ... - ValueError: unsigned oo times smaller number not defined + RuntimeError: indeterminate expression: 0 * infinity encountered. """ return self * ~other diff --git a/src/sage/symbolic/expression.pyx b/src/sage/symbolic/expression.pyx index 77aea5ad9df..a531a63619d 100644 --- a/src/sage/symbolic/expression.pyx +++ b/src/sage/symbolic/expression.pyx @@ -2710,7 +2710,7 @@ cdef class Expression(CommutativeRingElement): sage: x*unsigned_infinity Traceback (most recent call last): ... - ValueError: oo times number < oo not defined + RuntimeError: indeterminate expression: infinity * f(x) encountered. sage: SR(oo)*SR(oo) +Infinity diff --git a/src/sage/symbolic/ring.pyx b/src/sage/symbolic/ring.pyx index d883c622d91..9512b984633 100644 --- a/src/sage/symbolic/ring.pyx +++ b/src/sage/symbolic/ring.pyx @@ -103,6 +103,8 @@ cdef class SymbolicRing(CommutativeRing): 2 sage: SR.coerce(-infinity) -Infinity + sage: SR.coerce(unsigned_infinity) + Infinity sage: SR.has_coerce_map_from(ZZ['t']) True sage: SR.has_coerce_map_from(ZZ['t,u,v']) @@ -117,6 +119,8 @@ cdef class SymbolicRing(CommutativeRing): True sage: SR.has_coerce_map_from(GF(9, 'a')) True + sage: SR.has_coerce_map_from(UnsignedInfinityRing) + True TESTS: @@ -163,7 +167,8 @@ cdef class SymbolicRing(CommutativeRing): from sage.rings.polynomial.multi_polynomial_ring import is_MPolynomialRing from sage.rings.all import (ComplexField, - RLF, CLF, AA, QQbar, InfinityRing) + RLF, CLF, AA, QQbar, InfinityRing, + UnsignedInfinityRing) from sage.rings.finite_rings.finite_field_base import is_FiniteField from sage.interfaces.maxima import Maxima @@ -180,7 +185,7 @@ cdef class SymbolicRing(CommutativeRing): elif is_PolynomialRing(R) or is_MPolynomialRing(R) or is_FractionField(R): base = R.base_ring() return base is not self and self.has_coerce_map_from(base) - elif (R is InfinityRing + elif (R is InfinityRing or R is UnsignedInfinityRing or is_RealIntervalField(R) or is_ComplexIntervalField(R) or is_IntegerModRing(R) or is_FiniteField(R)): return True diff --git a/src/sage/tests/french_book/recequadiff.py b/src/sage/tests/french_book/recequadiff.py index 6deec70d9e8..b94e6c23980 100755 --- a/src/sage/tests/french_book/recequadiff.py +++ b/src/sage/tests/french_book/recequadiff.py @@ -12,6 +12,9 @@ sage: x = var('x') sage: y = function('y', x) + sage: _C = SR.var("_C") + sage: _K1 = SR.var("_K1") + sage: _K2 = SR.var("_K2") Sage example in ./recequadiff.tex, line 179:: @@ -104,10 +107,8 @@ Sage example in ./recequadiff.tex, line 367:: - sage: solve(ed, y)[0].subs_expr(_C==5).rhs() - Traceback (most recent call last): - ... - NameError: name '_C' is not defined + sage: solve(ed, y)[0].subs_expr(_C == 5).rhs() + e^(-sqrt(-2*cos(x) + 10)) Sage example in ./recequadiff.tex, line 377:: @@ -129,16 +130,16 @@ sage: P = Graphics() sage: for k in range(1,20,2): - ... P += plot(solve(ed, y)[0].subs_expr(c==1+k/4).rhs(), x, -3, 3) + ... P += plot(solve(ed, y)[0].subs_expr(c == 1+k/4).rhs(), x, -3, 3) sage: P - Graphics object consisting of 11 graphics primitives + Graphics object consisting of 10 graphics primitives Sage example in ./recequadiff.tex, line 426:: sage: P = Graphics() sage: for j in [0,1]: ... for k in range(1,10,2): - ... f = solve(ed,y)[j].subs_expr(c==2+0.25*k).rhs() + ... f = solve(ed,y)[j].subs_expr(c == 2+0.25*k).rhs() ... P += plot(f, x, -3, 3) sage: P Graphics object consisting of 10 graphics primitives