update besselj

Author: heltonmc

src/U_polynomials.jl

@@ -40,10 +40,17 @@ function Uk_poly_In(p, v, p2, ::Type{Float32})
     u10 = 1 / 88580102706155225088000 * evalpoly(p2, (9745329584487361980740625, -1230031256571145165088463750, 27299183373230345667273718125, -246750339886026017414509498824, 1177120360439828012193658602930, -3327704366990695147540934069220, 5876803711285273203043452095250, -6564241639632418015173104205000, 4513386761946134740461797128125, -1745632061522350031610173343750, 290938676920391671935028890625))
     u11 = -1 / 27636992044320430227456000 * evalpoly(p2, (15237265774872558064250728125, -2321657500166464779660536015625, 62011003282542082472252466220875, -676389476843440422173605288596087, 3926191452593448964331218647028194, -13704902022868787041097596217578170, 30589806122850866110941936529264950, -44801790321820682384740638703320750, 42936745153513007436411401865860625, -25963913760458280822131901737878125, 8997860461116953237934638500359375, -1363312191078326248171914924296875))
     u12 = 1 / 39797268543821419527536640000 * evalpoly(p2, (120907703923613748239829527671875, -21882222767154197351962677311437500, 692277766674325563109617997687563750, -8958590476947726766450604043798559500, 62055079517573388459132793029571876461, -261201165596865827608687437905740780920, 714528665351965363868467348116538170900, -1314368459332124683504418467809129387000, 1642838631056253395899341314188134178125, -1378260730939829908037976894025628887500, 743739612850105971846901081692858843750, -233469939346545526651936774615688437500, 32426380464797989812768996474401171875))
+    u13 = -1 / 955134445051714068660879360000 * evalpoly(p2, (17438611142828905996129258798828125, -3698121486504259988094897605296209375, 136735019134677724428035696765082813750, -2069933923586966756183324291232117362250, 16843538631795795357786827345307534156063, -83924867223075156862785921508524155665245, 274983827478138958623041508409195988431140, -616410216242554698436702353237166008656700, 962926533925253374359704288384340809260875, -1049095945162229046324321461816274931715625, 782463969315283937856703223178540650343750, -381190503845282445953419057314665534156250, 109361210755577700442544717509565392265625, -14020668045586884672121117629431460546875))
+    u14 = 1 / 45846453362482275295722209280000 * evalpoly(p2, (5448320367052402487647812713291015625, -1338184074771428116079233795614103631250, 57170953417612444837142230812990944671875, -1000503839668383458999731491914516564625300, 9440449669103391509091075981237243128469201, -54857705817658080981995319669299096598482382, 211477117385619365164298957821904115470535555, -564850830044980230366140582063618983657685400, 1070439683260179398514645952824098811025619475, -1451823699927947602004297385351260623500372750, 1401302601668131482630653233972052729323190625, -940627071986145750405920450097257805227812500, 417630985812040040477569028872405152769921875, -110320224449895843354117801955418504167031250, 13133360053559028970728309756597440972265625))   
+    u15 = -1 / 9820310310243703368343697227776000000 * evalpoly(p2, (8178936810213560828419581728001773291015625, -2303431987527333128955769182299845911305390625, 112597271053778779753048514469995937998172890625, -2254933495791765108580529087615802250458013685625, 24403480234538299231733883413666768614198435948125, -163359140754958502896104062604202934925448173291477, 730367145705123976114617970888707594104468177381925, -2284251621937242886581917667066122422330060024456125, 5136561256208409671660362298619778869859994724706875, -8420533422834140468835467666391400380550043688871875, 10085018700249896522602267572484630409640764997271875, -8735135969643867540297524795790262235822823374296875, 5329871927856528282113994744458999865006055974609375, -2173722139119126857509156976742601742357422740234375, 532039967451707060045861997017872377315039521484375, -59115551939078562227317999668652486368337724609375))
 
-    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12))
+    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15))
 end
 
+
+
+    
+   
 #=
 function Uk_poly_Kn(p, v, p2, ::Type{T}) where T <: BigFloat
     u0 = one(T)

src/besselj.jl

@@ -159,12 +159,21 @@ Nu must be real.
 function _besselj(nu, x)
     nu == 0 && return besselj0(x)
     nu == 1 && return besselj1(x)
-    if x < 20.0
+    if x < 4.0
         if nu > 60.0
             return log_besselj_small_arguments_orders(nu, x)
         else
             return besselj_small_arguments_orders(nu, x)
         end
+    elseif x < 20.0
+        if nu > 45
+            return besselj_debye(nu, x)
+        else
+            v = nu + 45
+            jnu = besselj_debye(v, x)
+            jnup1 = besselj_debye(v+1, x)
+            return besselj_down_recurrence(x, jnu, jnup1, v, nu)[2]
+        end
     elseif x > 1.85*nu
         return besselj_large_argument(nu, x)
     end
@@ -190,8 +199,18 @@ function _besselj(nu, x)
             jnum1 = besselj_large_argument(v2 -1, x)
             return besselj_up_recurrence(x, jnu, jnum1, v2, nu)[2]
         end
-    elseif 2*x < nu
+    elseif 1.5*x < nu
         return besselj_debye(nu, x)
+    elseif nu < 1000
+        v = nu + ceil(nu*1.3)
+        jnu = besselj_debye(v, x)
+        jnup1 = besselj_debye(v+1, x)
+        return besselj_down_recurrence(x, jnu, jnup1, v, nu)[2]
+    else
+        v = nu + 500
+        jnu = besselj_debye(v, x)
+        jnup1 = besselj_debye(v+1, x)
+        return besselj_down_recurrence(x, jnu, jnup1, v, nu)[2]
     end
 end
 
@@ -204,12 +223,13 @@ function besselj_large_argument(v, x::T) where T
     return cos_sum(α, -xn)*b
 end
 
+# generally can only use for x < 4.0
 # this needs a better way to sum these as it produces large errors
 # only valid in non-oscillatory regime (v>1/2, 0<t<sqrt(v^2 - 0.25))
 # power series has premature underflow for large orders
 # gives wrong answers for x > 20.0 (might want to fix this)
 function besselj_small_arguments_orders(v, x::T) where T
-    MaxIter = 1000
+    MaxIter = 2000
     out = zero(T)
     a = (x/2)^v / gamma(v + one(T))
     t2 = (x/2)^2

test/besselj_test.jl

@@ -100,11 +100,10 @@ for (v, x) in vnu
 end
 
 ## test all numbers and orders for 0<nu<100
-x = 0.1:0.5:100.0
-nu = 2:100
+x = 0.1:3.0:800.0
+nu = 2:4:500
 for v in nu, xx in x
-    @show v, xx
-    @test isapprox(Bessels._besselj(BigFloat(v), BigFloat(xx)), SpecialFunctions.besselj(BigFloat(v), BigFloat(xx)), rtol=1e-12)
+    @test isapprox(Bessels._besselj(BigFloat(v), BigFloat(xx)), SpecialFunctions.besselj(BigFloat(v), BigFloat(xx)), rtol=1e-14)
 end
 
 # test half orders (SpecialFunctions does not give big float precision)