juliabessel-functions

matrix of special (besselj) functions


I'm new to julia so I would welcome some advice to improve the following function,

using SpecialFunctions

function rb(x, nu_max)

  bj = Array{Complex64}(length(x), nu_max)

  nu = 0.5 + (0:nu_max)
  # somehow dot broadcast isn't happy
  # bj .= [ besselj(_nu,_x)*sqrt(pi/2*_x) for _nu in nu, _x in x]

  bj   = [ besselj(_nu,_x)*sqrt(pi/2*_x) for _nu in nu, _x in x]

end

rb(1.0:0.1:2.0, 500)

basically, I'm not quite sure what's the recommended way to get a matrix over these two parameters (x and nu). The documentation doesn't offer much information, but I understand that the underlying fortran routine internally loops over nu, so I'd rather not do it again in the interest of performance.


Edit: I'm asked about the goal; it's to compute the Riccati-Bessel functions $j_1(x,\nu),h_1(x,\nu)$ for multiple values of $x$ and $\nu$.

I've stripped down stylistic questions from the original version to focus on this core issue.


Solution

  • This is a great example where you can take full advantage of broadcasting. It looks like you want the cartesian product between x and nu, where the rows are populated by the values of nu and the columns are x. This is exactly what broadcasting can do — you just need to reshape x such that it's a single row across many columns:

    julia> using SpecialFunctions
    
    julia> x = 1.0:0.1:2.0
    1.0:0.1:2.0
    
    julia> nu = 0.5 + (0:500)
    0.5:1.0:500.5
    
     # this shows how broadcast works — these are the arguments and their location in the matrix
    julia> tuple.(nu, reshape(x, 1, :))
    501×11 Array{Tuple{Float64,Float64},2}:
     (0.5, 1.0)    (0.5, 1.1)    …  (0.5, 1.9)    (0.5, 2.0)
     (1.5, 1.0)    (1.5, 1.1)       (1.5, 1.9)    (1.5, 2.0)
     (2.5, 1.0)    (2.5, 1.1)       (2.5, 1.9)    (2.5, 2.0)
     (3.5, 1.0)    (3.5, 1.1)       (3.5, 1.9)    (3.5, 2.0)
     ⋮                           ⋱                ⋮
     (497.5, 1.0)  (497.5, 1.1)     (497.5, 1.9)  (497.5, 2.0)
     (498.5, 1.0)  (498.5, 1.1)     (498.5, 1.9)  (498.5, 2.0)
     (499.5, 1.0)  (499.5, 1.1)     (499.5, 1.9)  (499.5, 2.0)
     (500.5, 1.0)  (500.5, 1.1)  …  (500.5, 1.9)  (500.5, 2.0)
    
    julia> bj = besselj.(nu,reshape(x, 1, :)).*sqrt.(pi/2*reshape(x, 1, :))
    501×11 Array{Float64,2}:
     0.841471    0.891207     0.932039     …  0.9463       0.909297
     0.301169    0.356592     0.414341        0.821342     0.870796
     0.0620351   0.0813173    0.103815        0.350556     0.396896
     0.00900658  0.0130319    0.0182194       0.101174     0.121444
     ⋮                                     ⋱               ⋮
     0.0         0.0          0.0             0.0          0.0
     0.0         0.0          0.0             0.0          0.0
     0.0         0.0          0.0             0.0          0.0
     0.0         0.0          0.0          …  0.0          0.0