# Plot Bessel functions with pgfplots

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.

• One thing could be to define your own functions for this purpose if you know the mathematical expression. Pure curiosity, is there a specfic reason for you to discard luatex ? Oct 25, 2018 at 14:57
• @BambOo The mathematical expressions are quite complicated, using (at least in one definition) sums and the Gamma function. The only reason about luatex is to not complicate the compilation with another tool. However, I just tried the code in the linked luatex answer (adding the "magic commen" % !TeX program = lualatex) and it doesn't compile. Oct 25, 2018 at 17:33
• @HenriMenke Updated the question. Oct 25, 2018 at 17:37
• Giving the option --shell-escape to LuaLaTeX instead of pdflatex makes the code in the linked luatex answer work. Oct 25, 2018 at 22:25

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")

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{%
\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}{
}
\end{axis}
\end{tikzpicture}

\end{document} • This generates an error: Package PGF Math Error: Could not parse input '' as a floating point number, sorry. The unreadable part was near ''. (in 'BesselJ(0,x)'). }. I am using TeXstudio 2.12.6 with Qt 5.9.5. It is your code with % !TeX program = lualatex as first line. Oct 25, 2018 at 17:42
• I'm not sure what you mean by “with LuaTeX and FFI”; what should one do with the error could not load library gsl? Oct 25, 2018 at 22:44
• @egreg That means you don't have the GNU Scientific Library (GSL) installed. If you are using Homebrew or MacPorts you can probably get it using brew install gsl or port install gsl. If it still does not work, you may try ffi.load("/path/to/libgsl.dylib") instead of ffi.load("gsl"). Oct 26, 2018 at 9:30
• @HenriMenke Thanks, it works, but I think the dependency should be better explained in your answer. Oct 26, 2018 at 11:40

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: 