Recursive function in pgfplots

Question

I'm trying to plot Legendre polynomials, which can be done very easily by calling a recursive function. For instance, in Python this would like like this:

def legendre(xi, p):
    if p==0:
        return 1.
    elif p==1:
        return xi
    elif p>=2:
        return (1/p)*((2*p-1)*xi* legendre(xi, p-1) + (1-p)*legendre(xi, p-2))

However, if I want to plot these in pgfplots and it has turned to be very hard. I tried to follow this answer, but because of the recursive nature to the function (which calls the function with two different polynomial values), I can't implement it the same way as in the posted answer for the Fibonacci numbers. This is what I have so far:

\documentclass[border=5pt]{standalone}
\usepackage{pgfplots}
    \usetikzlibrary{math}
\begin{document}
\begin{tikzpicture}[
    evaluate={
        function legendre(\x,\p) {
            if \p == 0 then {
                return 1;
            } else {
                return legendre2(\x, \p);
            };
        };
        function legendre2(\x, \p) {
            if \p == 1 then {
                return \x;
            } else {
                return (1/\p)*((2*\p-1)*\x* legendre2(\x, \p-1) + (1-\p)*legendre2(\x, \p-2))
            };
        };
    },
]
\begin{axis}
  \addplot+ [domain=-1:1] {legendre(x,0)};
  \addplot+ [domain=-1:1] {legendre(x,1)};
  \addplot+ [domain=-1:1] {legendre(x,2)};
\end{axis}
\end{tikzpicture}
\end{document}

which can plot the first 2 Legendre polynomials. For the 3rd I enter in an infinite loop because I'm not calling the right function. I was wondering if there is a workaround this (I did not find a way to put a simple or statement like if \p == 0 or \p == 1 then so I could discriminate in this case for the first two orders.

I prefer to do calculations using C++ and send the numbers in table form for plotting with \pgfplots. Unless you are using \pgfmathparseFPU, pgfmath is not accurate enough for serious computation. — John Kormylo, May 05 '21 at 15:06

Sergio Llorente · Accepted Answer · 2021-05-05T18:46:08.950

You are very close to a solution. The problem is that legendre2 has not to call legendre2 itself, but legendre:

\documentclass[border=5pt]{standalone}
\usepackage{pgfplots}
    \usetikzlibrary{math}
\begin{document}
\begin{tikzpicture}[
    evaluate={
        function legendre(\x,\p) {
            if \p == 0 then {
                return 1;
            } else {
                return legendre2(\x, \p);
            };
        };
        function legendre2(\x, \p) {
            if \p == 1 then {
                return \x;
            } else {
                return (1/\p)*((2*\p-1)*\x* legendre(\x, \p-1) + (1-\p)*legendre(\x, \p-2));
            };
        };
    },
]
\begin{axis}[no markers]
  \addplot+ [domain=-1:1] {legendre(x,0)};
  \addplot+ [domain=-1:1] {legendre(x,1)};
  \addplot+ [domain=-1:1] {legendre(x,2)};
  \addplot+ [domain=-1:1] {legendre(x,3)};
  \addplot+ [domain=-1:1] {legendre(x,4)};
  \addplot+ [domain=-1:1] {legendre(x,5)};
\end{axis}
\end{tikzpicture}
\end{document}

Anyway, I agree with John Kormylo. You can use lua, as here, or gnuplot:

\documentclass[border=5pt]{standalone}
\usepackage{pgfplots}
    \usetikzlibrary{math}
\begin{document}
\newcommand{\defgnuplotfunc}{
  leg(x,p) = 
  (p==0) ? 1 :
  (p==1) ? x : 
  (1./p)*( (2.*p -1.) * x *leg(x,p-1) + (1.-p)*leg(x,p-2));
}
\begin{tikzpicture}
\begin{axis}
  \addplot+  gnuplot[raw gnuplot] {\defgnuplotfunc plot [-1:1] leg(x,0)};
  \addplot+  gnuplot[raw gnuplot] {\defgnuplotfunc plot [-1:1] leg(x,1)};
  \addplot+  gnuplot[raw gnuplot] {\defgnuplotfunc plot [-1:1] leg(x,2)};
  \addplot+  gnuplot[raw gnuplot] {\defgnuplotfunc plot [-1:1] leg(x,3)};
\end{axis}
\end{tikzpicture}
\end{document}

i forgot to comment that fot using gnuplot from latex you must enable system calla with the -shell-escape option. — Sergio Llorente, May 05 '21 at 17:13
Your solution won’t be correct for higher orders since for these it has to call itself. — aaragon, May 05 '21 at 17:34
In fact, I don’t think what you wrote for p=3 is correct either. — aaragon, May 05 '21 at 17:35
Oh I see how it works! It just bounces back and forth until it gets the right p. Thanks! — aaragon, May 05 '21 at 19:14
The first solution takes a really long time. Do you know why that is? — aaragon, May 05 '21 at 19:28
I guess latex is not fast doing floating-point arithmetic. Note that it is also not accurate (may be an issue for higher orders). — Sergio Llorente, May 08 '21 at 16:37

Recursive function in pgfplots

1 Answers1