From 3942e7de88f31e4b211424f1f0f21ce4924b244d Mon Sep 17 00:00:00 2001 From: mmatera Date: Sun, 19 Mar 2023 12:41:40 -0300 Subject: [PATCH] more rules and docs --- mathics/autoload/rules/Bessel.m | 73 +++++++++++++++++++++++----- mathics/builtin/specialfns/bessel.py | 28 +++++++---- 2 files changed, 80 insertions(+), 21 deletions(-) diff --git a/mathics/autoload/rules/Bessel.m b/mathics/autoload/rules/Bessel.m index 84c0163f0..ac9be96fb 100644 --- a/mathics/autoload/rules/Bessel.m +++ b/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*) diff --git a/mathics/builtin/specialfns/bessel.py b/mathics/builtin/specialfns/bessel.py index 6d0ef5bb1..10017b09c 100644 --- a/mathics/builtin/specialfns/bessel.py +++ b/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)&", }