使用 pgfplots 绘制贝塞尔函数

Question

我需要使用 pgfplots 绘制第一类和第二类贝塞尔函数（J 和 Y）以及第一类和第二类修正贝塞尔函数（I 和 K），整数阶，从 0 阶到 5 阶。其他问题中已经提出了几种解决方案，但是：

Gnuplot 只允许绘制第一类和第二类贝塞尔函数，且只能绘制 0 阶和 1 阶贝塞尔函数，因此这个答案不适合。我既不想使用 pstricks，也不想使用 LuaTeX，所以这个答案和这个答案不适合。这个问题非常相似，但评论中没有答案也没有例子。

如何使用 pgfplots 实现此目的？也可以使用 numpy、scipy 和 Matplotlib 的外部解决方案，如下所示这个例子就可以了（不过该代码似乎只适用于 J 函数）。

---

我尝试使用最后一个链接中的代码例子，插入为文件\addplot。example-04.txt它适用于函数 J，但它会打印不可接受的值（e+09）函数 Y，该值非常大且接近 0 的负值。生成的代码example-04.txt是 Python，在这里几乎肯定是 OT。

Answer 1

使用 LuaLaTeX 和 FFI（需要 LuaJITTeX 或 LuaTeX ≥ 1.0.3），可以使用 GNU 科学库 (GSL)，它实现了所有的特殊贝塞尔函数。

用 排版--shell-escape。您必须安装 GSL。

\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

这是我的评论这里：删除部分代码，改为自己想要的函数；贝塞尔函数如下这里。为了绘制多个函数，让第i个函数具有x坐标xi和y坐标yi，使用您想要的函数。

\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}

输出运行可钙看起来像这样：