Plot Bessel functions with pgfplots

Question

I need to plot the Bessel functions of first and second kind (J and Y), and the modified Bessel functions of first and second kind (I and K), with integer order, from order 0 to order 5, using pgfplots. Several solutions have been already presented in other questions, but:

Gnuplot makes available for plots only the Bessel functions of the first and second kind and only of order 0 and 1, so this answer is not suitable. I would not like to use nor pstricks, neither LuaTeX, so also this answer and this answer are not suitable. This question is quite similar, but has no answer and no examples in the comments.

How to accomplish this with pgfplots? Also an external solution with numpy, scipy and Matplotlib as in this example would be ok (however, the code seems to work only for the J function).

---

I tried using the code in the last linked example, inserting as an \addplot the file example-04.txt. It works for the function J, but it prints unacceptable values (e+09) for function Y, which is very large and negative near 0. The code which generates example-04.txt is Python and here it would be almost certainly OT.

Answer 1

With LuaLaTeX and FFI (requires LuaJITTeX or LuaTeX ≥ 1.0.3) one can make use of the GNU Scientific Library (GSL) which implements all the special Bessel functions.

Typeset with --shell-escape. You have to have the GSL installed.

\documentclass{article}
\usepackage{pgfplots}
\usepackage{luacode}

\begin{luacode*}
local ffi = require("ffi")
gsl = ffi.load("gsl")

ffi.cdef[[
double gsl_sf_bessel_Jn(int n, double x);
double gsl_sf_bessel_Yn(int n, double x);
double gsl_sf_bessel_In(int n, double x);
double gsl_sf_bessel_Kn(int n, double x);
]]
\end{luacode*}

\newcommand\declarebesselfunction[1]{%
  \pgfmathdeclarefunction{Bessel#1}{2}{%
    \pgfmathfloatparsenumber{%
      \directlua{tex.print(gsl.gsl_sf_bessel_#1n(
        \pgfmathfloatvalueof{##1},\pgfmathfloatvalueof{##2}))}%
    }%
  }%
}

\declarebesselfunction{J}
\declarebesselfunction{Y}
\declarebesselfunction{I}
\declarebesselfunction{K}

\begin{document}

\begin{tikzpicture}
  \begin{axis}[samples=100,no marks,restrict y to domain=-3:3]
    \pgfplotsinvokeforeach{0,...,5}{
      \addplot+[domain=0:10] {BesselJ(#1,x)};
      \addplot+[domain=.001:10] {BesselY(#1,x)};
      \addplot+[domain=0:10] {BesselI(#1,x)};
      \addplot+[domain=.001:10] {BesselK(#1,x)};
    }
  \end{axis}
\end{tikzpicture}

\end{document}

Answer 2

This is what I had in mind for my comment here: remove some of the code and change to the function you want; Bessel functions are given here. For plotting multiple functions let the ith function have x coordinates xi and y coordinates yi, use the function that you want.

\documentclass{article}
\usepackage[usenames,dvipsnames]{xcolor}
\usepackage{pgfplots}
\usepackage{sagetex}
\pagestyle{empty}
\begin{document}
Plotting the Bessel function using pgfplots and sagetex.
\begin{sagesilent}
LowerX = 0
UpperX = 12
LowerY = -1
UpperY = 1.25
step = .01
Scale = 1.0
xscale=1.0
yscale=1.0
output = r""
output += r"\begin{tikzpicture}"
output += r"[line cap=round,line join=round,x=8.75cm,y=8cm]"
output += r"\begin{axis}["
output += r"grid = none,"
output += r"minor tick num=4,"
output += r"every major grid/.style={Red!30, opacity=1.0},"
output += r"every minor grid/.style={ForestGreen!30, opacity=1.0},"
output += r"height= %f\textwidth,"%(yscale)
output += r"width = %f\textwidth,"%(xscale)
output += r"thick,"
output += r"black,"
output += r"axis lines=center,"
output += r"domain=%f:%f,"%(LowerX,UpperX)
output += r"line join=bevel,"
output += r"xmin=%f,xmax=%f,ymin= %f,ymax=%f,"%(LowerX,UpperX,LowerY, UpperY)
#output += r"xticklabels=\empty,"
#output += r"yticklabels=\empty,"
output += r"major tick length=5pt,"
output += r"minor tick length=0pt,"
output += r"major x tick style={black,very thick},"
output += r"major y tick style={black,very thick},"
output += r"minor x tick style={black,thin},"
output += r"minor y tick style={black,thin},"
#output += r"xtick=\empty,"
#output += r"ytick=\empty"
output += r"]"
##FUNCTION 1
t1 =  var('t1')
x1_coords = srange(LowerX,UpperX,step)
y1_coords = [bessel_J(1, t1).n(digits=6) for t1 in x1_coords]
output += r"\addplot[thin, NavyBlue, unbounded coords=jump] coordinates {"
for i in range(0,len(x1_coords)):
    if (y1_coords[i])<LowerY or (y1_coords[i])>UpperY:
        output += r"(%f,inf) "%(x1_coords[i])
    else:
        output += r"(%f,%f) "%(x1_coords[i],y1_coords[i])
output += r"};"
##FUNCTION 2
t2 =  var('t2')
x2_coords = srange(LowerX,UpperX,step)
y2_coords = [bessel_J(0, t2).n(digits=6) for t2 in x2_coords]
output += r"\addplot[thin, red, unbounded coords=jump] coordinates {"
for i in range(0,len(x2_coords)):
    if (y2_coords[i])<LowerY or (y2_coords[i])>UpperY:
        output += r"(%f,inf) "%(x2_coords[i])
    else:
        output += r"(%f,%f) "%(x2_coords[i],y2_coords[i])
output += r"};"
##### COMMENT OUT A LINE OF SAGESILENT BY STARTING WITH #
output += r"\end{axis}"
output += r"\end{tikzpicture}"
\end{sagesilent}
\begin{center}
\sagestr{output}
\end{center}
\end{document}

The output running in Cocalc looks like this: