1
1
``` julia
2
- using FiniteDifferences, ArbNumerics, DelimitedFiles, SpecialFunctions
3
-
4
- if ! (@isdefined vx)
5
- vgrid = sort (vcat (range (0.05 , 15.0 , length= 20 ), [0.5 , 1.5 , 2.5 , 3.5 ]))
6
- xgrid = range (15.0 , 30.0 , length= 20 )
7
- const vx = vec (collect (Iterators. product (vgrid, xgrid)))
8
- end
2
+ using ArbNumerics, DelimitedFiles, SpecialFunctions
9
3
4
+ # Because you can't just use FiniteDifferences due to the "numerical noise".
10
5
simplefd (f,x,h= ArbReal (1e-40 )) = (f (x+ h)- f (x))/ h
11
6
12
7
ArbNumerics. setprecision (ArbReal, digits= 100 )
8
+ arb_besselk (v,x) = ArbNumerics. besselk (ArbReal (v), ArbReal (x))
9
+ arb_besselkx (v,x) = arb_besselk (v,x)* ArbNumerics. exp (ArbReal (x))
13
10
14
- function arb_besselkx (v,x)
15
- (av, ax) = ArbReal .((v,x))
16
- ArbNumerics. besselk (av, ax)* ArbNumerics. exp (ax)
17
- end
11
+ #
12
+ # besselkx_levin test:
13
+ #
14
+ vgrid = sort (vcat (range (0.05 , 15.0 , length= 20 ), [0.5 , 1.5 , 2.5 , 3.5 ]))
15
+ xgrid = range (15.0 , 30.0 , length= 20 )
16
+ vx = vec (collect (Iterators. product (vgrid, xgrid)))
18
17
19
18
ref_values = map (vx) do vxj
20
19
(v,x) = vxj
@@ -29,4 +28,25 @@ out_matrix = hcat(getindex.(vx, 1), # test v argument
29
28
getindex .(ref_values, 2 )) # test d/dx value
30
29
31
30
writedlm (" besselkx_levin_enzyme_tests.csv" , out_matrix)
31
+
32
+ #
33
+ # besselk_power_series test:
34
+ #
35
+ vgrid = [- 1e-8 , 1e-8 , 1 - 1e-8 , 1.0 , 1 + 1e-8 , 2 - 1e-8 , 2.0 , 2 + 1e-8 , 3 - 1e-8 , 3.0 , 3 + 1e-8 ]
36
+ xgrid = range (1e-5 , 1.5 , length= 20 )
37
+ vx = vec (collect (Iterators. product (vgrid, xgrid)))
38
+
39
+ ref_values = map (vx) do vxj
40
+ (v,x) = vxj
41
+ dx = simplefd (_x-> arb_besselk (v, _x), x) # NOT besselkx!
42
+ dv = simplefd (_v-> arb_besselk (_v, x), v) # NOT besselkx!
43
+ Float64 .((dv, dx))
44
+ end
45
+
46
+ out_matrix = hcat (getindex .(vx, 1 ), # test v argument
47
+ getindex .(vx, 2 ), # test x argument
48
+ getindex .(ref_values, 1 ), # test d/dv value
49
+ getindex .(ref_values, 2 )) # test d/dx value
50
+
51
+ writedlm (" besselk_power_series_enzyme_tests.csv" , out_matrix)
32
52
```
0 commit comments