more rules and docs

mathics/autoload/rules/Bessel.m

@@ -1,26 +1,77 @@
-(* Adapted from symja_android_library/symja_android_library/rules/Bessel{I,J}.m
-   Note: These are not currently covered by SymPy.
- *)
+(*Extended rules for handling expressions with Bessel functions*)
+
+
+
+Unprotect[HankelH1];
+(*HankelH1[x_Integer?NegativeQ, z_]:=-HankelH1[-x, z];*)
+(*Limit cases*)
+HankelH1[nu_, 0] := DirectedInfinity[];
+Protect[HankelH1];
+
+
+Unprotect[HankelH2];
+(*HankelH2[x_Integer?NegativeQ, z_]:=-HankelH2[-x, z];*)
+(*Limit cases*)
+HankelH2[nu_,0] := DirectedInfinity[];
+Protect[HankelH2];
 
 
 Unprotect[BesselI]
-BesselI[nu_/;(nu>0 && IntegerQ[2*nu]),z_]:=Module[{u,f,k= nu-1/2},f=Sinh[u]/u;While[k>0, k=k-1;f = (-D[f, u]/u)]; (Sqrt[2/Pi z] * ((-u)^(nu-1/2)*f))/.u->z];
-BesselI[nu_/;(nu<0 && IntegerQ[2*nu]),z_]:=Module[{u,f,k=-nu-1/2},f=Cosh[u]/u;While[k>0, k=k-1;f = (-D[f, u]/u)]; (Sqrt[2/Pi z] * ((-u)^(-nu-1/2)*f))/.u->z];
+(*Rayleight's formulas for half-integer indices*)
+BesselI[nu_/;(nu>0 && IntegerQ[2*nu]), z_]:=Module[{u,f,k= nu-1/2},f=Sinh[u]/u;While[k>0, k=k-1;f = (-D[f, u]/u)]; (Sqrt[2/Pi z] * ((-u)^(nu-1/2)*f))/.u->z];
+BesselI[nu_/;(nu<0 && IntegerQ[2*nu]), z_]:=Module[{u,f,k=-nu-1/2},f=Cosh[u]/u;While[k>0, k=k-1;f = (-D[f, u]/u)]; (Sqrt[2/Pi z] * ((-u)^(-nu-1/2)*f))/.u->z];
+(*Limit cases*)
+BesselI[0, 0] := 1;
+BesselI[nu_Integer,0]:=0;
+BesselI[nu_Rational, 0] := If[nu>0, 0, DirectedInfinity[]];
+BesselI[nu_Real, 0] := If[nu>0, 0, DirectedInfinity[]];
+BesselI[nu_, DirectedInfinity[z___]] := 0;
 Protect[BesselI]
 
 Unprotect[BesselK]
-BesselK[nu_/;(nu>0 && IntegerQ[2*nu]),z_]:=Module[{u,f,k= nu-1/2},f=Exp[-u]/u;While[k>0, k=k-1;f = (D[f, u]/u)]; (Sqrt[Pi/2 z] * ((-u)^(nu-1/2)*f))/.u->z];
-BesselK[nu_/;(nu<0 && IntegerQ[2*nu]),z_]:=Module[{u,f,k=-nu-1/2},f=Exp[-u]/u;While[k>0, k=k-1;f = (D[f, u]/u)]; (Sqrt[Pi/2 z] * ((-u)^(-nu-1/2)*f))/.u->z];
+(*Rayleight's formulas for half-integer indices*)
+BesselK[nu_/;(nu>0 && IntegerQ[2*nu]), z_]:=Module[{u,f,k= nu-1/2},f=Exp[-u]/u;While[k>0, k=k-1;f = (D[f, u]/u)]; (Sqrt[Pi/2 z] * ((-u)^(nu-1/2)*f))/.u->z];
+BesselK[nu_/;(nu<0 && IntegerQ[2*nu]), z_]:=Module[{u,f,k=-nu-1/2},f=Exp[-u]/u;While[k>0, k=k-1;f = (D[f, u]/u)]; (Sqrt[Pi/2 z] * ((-u)^(-nu-1/2)*f))/.u->z];
+(*Limit cases*)
+BesselK[0, 0] = DirectedInfinity[-1];
+BesselK[nu_?NumericQ, 0] = DirectedInfinity[];
 Protect[BesselK]
 
 
 Unprotect[BesselJ]
-BesselJ[nu_/;(nu>0 && IntegerQ[2*nu]),z_]:=Module[{u,f,k= nu-1/2},f=Sin[u]/u;While[k>0, k=k-1;f = (D[f, u]/u)]; (Sqrt[2/Pi z] * ((-u)^(nu-1/2)*f))/.u->z];
-BesselJ[nu_/;(nu<0 && IntegerQ[2*nu]),z_]:=Module[{u,f,k=-nu-1/2},f=Cos[u]/u;While[k>0, k=k-1;f = (-D[f, u]/u)]; (Sqrt[2/Pi z] * ((-u)^(-nu-1/2)*f))/.u->z];
+(*Rayleight's formulas for half-integer indices*)
+BesselJ[nu_/;(nu>0 && IntegerQ[2*nu]), z_]:=Module[{u,f,k= nu-1/2},f=Sin[u]/u;While[k>0, k=k-1;f = (D[f, u]/u)]; (Sqrt[2/Pi z] * ((-u)^(nu-1/2)*f))/.u->z];
+BesselJ[nu_/;(nu<0 && IntegerQ[2*nu]), z_]:=Module[{u,f,k=-nu-1/2},f=Cos[u]/u;While[k>0, k=k-1;f = (-D[f, u]/u)]; (Sqrt[2/Pi z] * ((-u)^(-nu-1/2)*f))/.u->z];
+(*Limit cases*)
+BesselJ[0, 0] := 1;
+BesselJ[nu_Integer,0]:=0;
+BesselJ[nu_Rational, 0] := If[nu>0, 0, DirectedInfinity[]];
+BesselJ[nu_Real, 0] := If[nu>0, 0, DirectedInfinity[]];
+BesselJ[nu_, DirectedInfinity[z___]] := 0;
 Protect[BesselJ]
 
 
 Unprotect[BesselY]
-BesselY[nu_/;(nu>0 && IntegerQ[2*nu]),z_]:=Module[{u,f,k= nu-1/2},f=Cos[u]/u;While[k>0, k=k-1;f = (D[f, u]/u)]; (-Sqrt[2/Pi z] * ((-u)^(nu-1/2)*f))/.u->z];
-BesselY[nu_/;(nu<0 && IntegerQ[2*nu]),z_]:=Module[{u,f,k=-nu-1/2},f=Sin[u]/u;While[k>0, k=k-1;f = (D[f, u]/u)]; (Sqrt[2/Pi z] * ((u)^(-nu-1/2)*f))/.u->z];
+(*Rayleight's formulas for half-integer indices*)
+BesselY[nu_/;(nu>0 && IntegerQ[2*nu]), z_]:=Module[{u,f,k= nu-1/2},f=Cos[u]/u;While[k>0, k=k-1;f = (D[f, u]/u)]; (-Sqrt[2/Pi z] * ((-u)^(nu-1/2)*f))/.u->z];
+BesselY[nu_/;(nu<0 && IntegerQ[2*nu]), z_]:=Module[{u,f,k=-nu-1/2},f=Sin[u]/u;While[k>0, k=k-1;f = (D[f, u]/u)]; (Sqrt[2/Pi z] * ((u)^(-nu-1/2)*f))/.u->z];
+(*Limit cases*)
+BesselY[0, 0] = DirectedInfinity[-1];
+BesselY[nu_, 0] = DirectedInfinity[];
 Protect[BesselY]
+
+
+
+
+
+Unprotect[Integrate];
+(*  See https://dlmf.nist.gov/10.9 *)
+Integrate[Cos[z_Real Cos[Theta_]], {Theta_, 0, Pi}]:= Pi BesselJ[0, Abs[z]];
+
+
+(* This rule needs to implement Elements*)
+Integrate[Cos[z_ Cos[Theta_]], {Theta_, 0, Pi}]:= ConditionalExpression[Pi BesselJ[0, Abs[z]], Element[z, Reals]];
+
+Protect[Integrate];
+
+(*TODO: extend me with series expansions, integrals, etc*)

mathics/builtin/specialfns/bessel.py

@@ -344,6 +344,10 @@ class BesselI(_Bessel):
 
         >> Plot[BesselI[0, x], {x, 0, 5}]
          = -Graphics-
+
+        The special case of half-integer index is expanded using Rayleigh's formulas:
+        >> BesselI[3/2, x]
+         = Sqrt[2] Sqrt[x] (-Sinh[x] / x ^ 2 + Cosh[x] / x) / Sqrt[Pi]
     """
 
     mpmath_name = "besseli"
@@ -388,12 +392,10 @@ class BesselJ(_Bessel):
 
     >> Plot[BesselJ[0, x], {x, 0, 10}]
      = -Graphics-
-    """
 
-    # TODO: Sympy Backend is not as powerful as Mathematica
-    """
+    The special case of half-integer index is expanded using Rayleigh's formulas:
     >> BesselJ[1/2, x]
-     = Sqrt[2 / Pi] Sin[x] / Sqrt[x]
+     = Sqrt[2] Sin[x] / (Sqrt[x] Sqrt[Pi])
     """
 
     mpmath_name = "besselj"
@@ -428,6 +430,11 @@ class BesselK(_Bessel):
 
     >> Plot[BesselK[0, x], {x, 0, 5}]
      = -Graphics-
+
+    The special case of half-integer index is expanded using Rayleigh's formulas:
+    >> BesselK[-3/2, x]
+     = Sqrt[2] Sqrt[x] Sqrt[Pi] (E ^ (-x) / x ^ 2 + E ^ (-x) / x) / 2
+
     """
 
     mpmath_name = "besselk"
@@ -460,19 +467,20 @@ class BesselY(_Bessel):
     >> BesselY[1.5, 4]
      = 0.367112
 
-    ## Returns ComplexInfinity instead
-    ## #> BesselY[0., 0.]
-    ##  = -Infinity
+    >> BesselY[0., 0.]
+      = -Infinity
 
     >> Plot[BesselY[0, x], {x, 0, 10}]
      = -Graphics-
-    """
 
-    # TODO: Special Values
-    """
+    The special case of half-integer index is expanded using Rayleigh's formulas:
+    >> BesselY[-3/2, x]
+     =  Sqrt[2] Sqrt[x] (-Sin[x] / x ^ 2 + Cos[x] / x) / Sqrt[Pi]
+
     >> BesselY[0, 0]
      = -Infinity
     """
+
     rules = {
         "Derivative[0,1][BesselY]": "(BesselY[-1 + #1, #2] / 2 - BesselY[1 + #1, #2] / 2)&",
     }