Recursion == Looping
The Scheme programming language does not have an explicit loop construct. Iteration in Scheme is done by recursion. Recursion requires only the ability to define a function and a conditional operator. Gnuplot's language is quite limited but it contains both those features.
In the following examples I show how to use recursion in Gnuplot to do
some things that might seem to require a second scripting
language. The more involved examples use recursion to build a command
string and then use eval to evaluate the command string. To debug
these examples, one can use print in place of eval.
A Simple Recursion
This is the simplest example of a recursive function:
1. factorial(n) = (n < 2 ? 1 : n*factorial(n-1)); 2. print factorial(5);
120
A Simple Tail Recursion
The first example does not use Tail recursion. Here is an example that does
1. factorial(n, s) = (n < 2 ? s : factorial(n-1, n*s)); 2. print factorial(10, 1);
3628800
Tail Recursion is Better!
Gnuplot has a built-in stack limit of 250. Let us compare a
tail-recursive tr versus a non tail-recursive t function that computes
$\displaystyle\sum\limits_{n=0}^x n$.
non Tail-Recursive
The non tail-recursive function t will give us a stack-overflow error
for a very small stack size.
1. t(x) = (x < 1 ? 0 : x+t(x-1)) 2. print t(98) 3. print t(99)
4851
line 0: stack overflow
Tail Recursive
The tail-recursive function tr will not result in a stack overflow.
1. tr(x, s) = (x < 1 ? s : tr(x-1, s+x)) 2. print tr(20000,0)
200010000
Ok, that is true, but it may bust other things.
1. tr(x, s) = (x < 1 ? s : tr(x-1, s+x)) 2. looptr(a, b, cmd) = (b < a ? cmd : looptr(a+1, b, sprintf("%s;\n print 'a=%g ', tr(%g,0) ", cmd, a, a))) 3. eval looptr(20135,20155,"")
a=20135 202719180
a=20136 202739316
a=20137 202759453
a=20138 202779591
a=20139 202799730
a=20140 202819870
a=20141 202840011
a=20142 202860153
a=20143 202880296
a=20144 /bin/bash: line 6: 3668 Segmentation fault gnuplot <<EOF
tr(x, s) = (x < 1 ? s : tr(x-1, s+x))
looptr(a, b, cmd) = (b < a ? cmd : looptr(a+1, b, sprintf("%s;\n print 'a=%g ', tr(%g,0) ", cmd, a, a)))
eval looptr(20135,20155,"")
EOF
The segmentation fault is non-deterministic. It is caused by tr, not
looptr, as can be verified by hand.
A Tail-Recursive Plot Command
In this example, more akin to many Scheme examples, we keep track of the state ($s$) at each step, in addition to the step itself ($n$). We count backward from $n$ to 0, and return the accumulated state $s$ at that point. At each step, we print part of a plot command line, which is passed as the current state to plotcmd with $n$ reduced by 1.
The output shows that plotcmd recurses till $n=0$, at which point the state $s$ is the initial value passed, which is returned and printed on the left for the $n=1$ recursion, followed by $\sin(1*x)$ and the plot options, which is saved in the state $s$ for the $n=2$ recursion, and so on.
| Step ($n$) | Input State ($s$) |
|---|---|
| 3 | "plot [0:2*pi] 0" |
| 2 | "plot [0:2*pi] 0, sin(3*x) linewidth 3" |
| 1 | "plot [0:2*pi] 0, sin(3*x) linewidth 3, sin(2*x) linewidth 3" |
| 0 | "plot [0:2*pi] 0, sin(3*x) linewidth 3, sin(2*x) linewidth 3, sin(1*x) linewidth 3" |
1. plotcmd(s, n, options) = (n < 1 ? s : \ 2. plotcmd( \ 3. sprintf("%s, sin(%g*x) %s", s, n, options), \ 4. n-1, options)) 5. 6. set key outside 7. set terminal svg 8. set output 'mathjax-html/plotcmd-tr.svg' 9. options = "linewidth 3" 10. plotcmd = plotcmd("plot [0:2*pi] 0", 3, options) 11. eval plotcmd 12. print plotcmd
plot [0:2*pi] 0, sin(3*x) linewidth 3, sin(2*x) linewidth 3, sin(1*x) linewidth 3
 |
| Multiple plots of sin. |
Remark
I would be remiss not to mention that in version 4.4 of Gnuplot, there is a limited looping ability with plot for.
1. set key outside 2. set terminal svg 3. set output 'mathjax-html/plot-for.svg' 4. set xrange [0:2*pi] 5. plot for [a=0:3] sin(a*x) linewidth 3 title sprintf("sin(%g*x)", a)
 |
| The plot for command. |
A non Tail-Recursive Variant
In this example, we directly keep track of the state. The output shows that plotcmd recurses till $n=0$, at which point the state $s$ is the initial value passed, which is returned and printed on the left for the $n=1$ recursion, followed by $\sin(1*x)$ and the plot options, which is saved in the state $s$ for the $n=2$ recursion, and so on.
| Step ($n$) | State ($s$) |
|---|---|
| 0 | "plot [0:2*pi] 0" |
| 1 | "plot [0:2*pi] 0, sin(1*x) linewidth 3" |
| 2 | "plot [0:2*pi] 0, sin(1*x) linewidth 3, sin(2*x) linewidth 3" |
| 3 | "plot [0:2*pi] 0, sin(1*x) linewidth 3, sin(2*x) linewidth 3, sin(3*x) linewidth 3" |
1. plotcmd(s, n, options) = (n < 1 ? s : \ 2. sprintf("%s, sin(%g*x) %s", \ 3. s=plotcmd(s, n-1, options), \ 4. n, options)); 5. 6. set key outside 7. set terminal svg 8. set output 'mathjax-html/plotcmd.svg' 9. options = "linewidth 3" 10. plotcmd = plotcmd("plot [0:2*pi] 0", 3, options) 11. eval plotcmd 12. print plotcmd
plot [0:2*pi] 0, sin(1*x) linewidth 3, sin(2*x) linewidth 3, sin(3*x) linewidth 3
 |
| Multiple plots of sin. |
An Example using Riemann Sums
The final example comes from calculus, where we want to illustrate the
Riemann integral. At each iteration of the counter $c$, rectangles
appends the code to construct a single rectangle with bottom left
vertex at $(x, 0)$ and upper right vertex at $(x+w, y)$. We use the
counter $c$ to index the rectangles so that Gnuplot does not try to
put more rectangles than needed on the graph.
1. rectangles(a, w, c, cmd) = (c < 1 ? cmd : \ 2. rectangles(a+w, w, c-1, \ 3. sprintf("set object %g rectangle from %g, 0 to %g, %g;\n%s", \ 4. c, a, a+w, f(a), cmd))) 5. 6. set terminal svg 7. set output 'mathjax-html/riemann-sum.svg' 8. set style rectangle fillcolor rgb 'blue' back 9. f(x) = x**3 10. eval rectangles(0, 1.0/5, 5, "") 11. plot [0:1] f(x) notitle linewidth 3
On line 10, the function rectangles accumulates the commands to place
5 rectangles in the unit interval, and eval evaluates the accumulated
commands.
We could have put the command on line 11 into the final argument of
rectangles on line 10.
 |
| A plot of $f(x)=x^3$ and the left-hand rule for approximating $\displaystyle\int_0^1 f(x) \D{x}$. |
A More Elaborate Example
The previous example is somewhat spartan. It does not compute the Riemann sum, and the left endpoint rule is hard-coded. This example fixes those shortcomings. Just for flair, we also show how the color of the rectangles can be programmatically determined.
1. leftrule(a, b) = f(a); 2. rightrule(a, b) = f(b); 3. midrule(a, b) = 0.5*(f(a)+f(b)); 4. rule(a, b) = (rule==0 ? leftrule(a, b) : (rule==1 ? rightrule(a, b) : midrule(a, b))); 5. 6. rectoptions(y) = sprintf("fillcolor rgb '%s' back;", (y < 0.5 ? 'blue' : 'green')) 7. 8. rectangle(x, y, w, c) = sprintf("%s \ 9. set object %g rectangle from %g, 0 to %g, %g %s; \ 10. area = area+%g*%g;", \ 11. (c<2 ? "area=0;" : ""), c, x, x+w, y, rectoptions(y), w, y); 12. 13. rectangles(a, w, c, cmd) = (c < 1 ? cmd : \ 14. rectangles(a+w, w, c-1, \ 15. sprintf("%s\n%s", \ 16. rectangle(a, rule(a, a+w), w, c), cmd))) 17. 18. f(x) = x**3; N = 10; T = 1.0; W = T/N; 19. 20. set macros 21. set xtics 0, W, T; set ytics 0, 0.2, f(T) 22. set xrange [0:T]; set yrange [0:f(T)] 23. set terminal svg size 1500, 480 24. set output 'mathjax-html/riemann-sums.svg' 25. set multiplot layout 1, 3 26. label = "set label 1 sprintf('colored\narea=%.3f', area) at 0.5,0.9 center front;" 27. plotf = "@label;set title title;plot f(x) notitle linewidth 3, 0.5 notitle linewidth 3;" 28. rule = 0;title='left rule' ; eval rectangles(0, W, N, plotf) 29. rule = 1;title='right rule'; eval rectangles(0, W, N, plotf) 30. rule = 2;title='mid rule'; eval rectangles(0, W, N, plotf) 31. unset multiplot
 |
| The left-, right-, and mid-point rules for the Riemann integral of $f(x)=x^3$ on the unit interval. |
![[FSF Associate Member 10686]](http://static.fsf.org/nosvn/associate/fsf-10686.png)
