Jun

8

Getting started with gnuplot

Filed Under Applied Math, Software, Suggested reading, Tutorial | 4 Comments

gnuplot is an excellent scientific package for visualizing data and plotting functions. Despite its name, it has nothing to do with the GNU project, even though it’s Open Source and entirely free. This tool is very handy whenever you need to produce production quality graphics from a given data-set (or function). It is no wonder that in its 20+ years of existence, it’s been employed in all sorts of industries.

Flexible, powerful, and easy to use, gnuplot is not only handy as a stand-alone program but can be used successfully by a variety of programming languages, including but not limited to Ruby, Python and Perl. It has also been adopted as a plotting engine by Open Source programs like Maxima and GNU Octave. The minimal effort required to learn gnuplot is therefore a very worthwhile endeavor. Being able to visualize things as you analyze data and explore mathematics, is a very useful aid.

gnuplot is cross-platform, and you can download it from SourceForge.

Once you have installed gnuplot, you can start the tool by simply running gnuplot from the shell. Getting started with the program is straightforward, given that the most basic functionalities are fairly intuitive. For example, if you wanted to plot the sin(x) function, you could run:

gnuplot> plot sin(x)

gnuplot1.png

gnuplot will automatically decide for you what portions of the function should be visualized (in our case between -10 and 10 on the x axis). If you want to overwrite this, you can specify otherwise, as shown below:

gnuplot> plot [-pi:pi] sin(x)

gnuplot2.png

Similarly, you can customize the plotting, by specifying a range on the y axis:

gnuplot> plot [] [-0.5:0.5] sin(x)

gnuplot3.png

It is also possible to plot several functions at once:

gnuplot> plot [-2:2] x, x**2, x**3

gnuplot4.png

gnuplot can draw all sorts of functions and graphics. For example, the following plots two surfaces:

gnuplot> splot x**2+y**2, x**2-y**2

gnuplot5.png

What you see plotted above, is based on a sensible set of defaults, but gnuplot offers countless options to customize the appearance of your graphics.

gnuplot isn’t just useful for plotting mathematical functions. I often use it for plotting data that is stored in tabular format in simple text formats. For example, assume that you have collected the following data on Deaths by Major Causes in the US, from a website:

stats.txt

That’s a big chunk of data, but it’s not very easy to analyze. It’d be nice to be able to visualize it. Assuming you saved it in a stats.txt file (you can include or remove the first few lines of comments) from gnuplot you can run the following:

gnuplot> plot "stats.txt"

gnuplot6.png

That’s a meaningless mess. But gnuplot is very flexible, so we can easily do much better. Let’s use the first column for x values, and each of the remaining columns as a curve of its own by running:


gnuplot> plot "stats.txt" using 1:2 title "Heart disease" with lines, "stats.txt" using 1:3 title "Cancer" with lines, "stats.txt" using 1:4 title "Cerebro-vascular diseases" with lines, "stats.txt" using 1:5 title "Lower respiratory diseases" with lines, "stats.txt" using 1:6 title "Diabetes mellitus" with lines, "stats.txt" using 1:7 title "Influenza and pneumonia" with lines, "stats.txt" using 1:8 title "Chronic liver disease" with lines, "stats.txt" using 1:9 title "Accidents" with lines, "stats.txt" using 1:10 title "Suicide" with lines, "stats.txt" using 1:11 title "Homicide" with lines

Or it’s abbreviated version:

gnuplot> plot "stats.txt" u 1:2 title "Heart disease" w l, "" u 1:3 title "Cancer" w l, "" u 1:4 title "Cerebro-vascular diseases" w l, "" u 1:5 title "Lower respiratory diseases" w l, "" u 1:6 title "Diabetes mellitus" w l, "" u 1:7 title "Influenza and pneumonia" w l, "" u 1:8 title "Chronic liver disease" w l, "" u 1:9 title "Accidents" w l, "" u 1:10 title "Suicide" w l, "" u 1:11 title "Homicide" w l

gnuplot7.png

Much better! Of course, the command above has a bit of repetition in it, but it’s justified by the flexibility of being able to pull data from different files if required, and plot some data with lines and others without (or in different styles) if we want to. To improve the appearance and legibility of the chart, let’s add a grid:

gnuplot> set grid
gnuplot> replot

gnuplot8.png

Nice! And notice how the replot command runs the last plotting directive for us, saving us some typing (or scrolling through the previous commands). To improve this further, we can also set the labels for the axis and a title for the whole chart:

gnuplot> set xlabel "Years"
gnuplot> set ylabel "Death rates per 100,000 people"
gnuplot> set title "Deaths by Major Causes, 1960–2005"
gnuplot> replot

gnuplot9.png
(Click here for a larger version)

Effortlessly, we went from having some data in a text file, to a meaningful, and professional looking chart that could be used in an article or a book. Not bad. You can save it in your favorite format, by setting the output file name and the terminal type. In the example below, I specify the font type, size and its location on my Mac OS X system, but it’s entirely optional. If missing, the default fonts apply.

gnuplot> set output "mychart.png"
gnuplot> set terminal png font "/Library/Fonts/Times New Roman.ttf, 11"
gnuplot> replot

This very short introduction should you give you a glimpse into what gnuplot can do for you, and even get you started with the tool. But gnuplot is much more than this and infinitely customizable, so you really may want to consider learning more about it. Unfortunately, despite being widely used for more than a decade, gnuplot has never had a book published about it.

gnuplot.jpgThankfully, Manning Publications will be putting out a book called Gnuplot in Action in October. The good news is that you don’t have to wait for the dead tree version, you can purchase it today through the Manning Early Access Program (MEAP). 12 out of the 15 chapters are already available, so the book is pretty much complete. If you’re interested in making the best out of gnuplot, this book is a must have. You don’t have to be a mathematician or a programmer to follow along, and it’s so much more than a bunch of options for customizing your graphics. The book guides you through the best techniques for taking advantage of the tool in several common scenarios. By reading, this easy to follow book, you’ll be able to generate colorful, nifty plots and have good mastery of the tool in a very short amount of time.

Mar

6

Polynomial Root-finding with the Jenkins-Traub Algorithm

Filed Under Applied Math, Tutorial | 3 Comments

The Jenkins-Traub Algorithm is a standard in the field of numerical computation of polynomial roots, fundamentally developed as a numerical algorithm specifically for the task of computing polynomial roots. In other words, (i) because it was planned from the outset for numerical purposes rather than being simply an adaptation of an analytic formula, it is extremely robust, effectively minimizing the effects of computer round-off error, while (ii) also being extremely efficient compared to more general methods not written specifically for the task of computing polynomial roots; in fact, the algorithm converges to polynomial roots at a rate better than quadratic. Furthermore, since being introduced over thirty years ago, the algorithm has had time to be rigorously tested and has successfully proven its quality; as a result, it has gained a wide distribution as evidenced by its incorporation in commercial software products and the posting on the NETLIB website of source code for programs based on the algorithm.

An Introduction to the Field of Root-finding

Given a function which is defined in terms of one or more independent variables, the roots (also called zeros) of the equation are the values of the independent variable(s) for which the function equals 0:

\displaystyle f(z1, z2, z3, z4, z5, . . . ) = 0

Note that, in general, the z values are complex numbers, comprised of real and imaginary components (indicated by the i): z = x + iy.

Consider the following equation:

\displaystyle 2z_{1} + 0.5z_{2}z_{3} + 3.14z_{4}^{2}+1.414\ln(z_{5})=23

The values of z_{1}, z_{2}, z_{3}, z_{4} and z_{5} which satisfy this equation may not be possible by analytical methods, so the equation would be rearranged into the following form:

\displaystyle 2z_{1} + 0.5z_{2}z_{3} + 3.14z_{4}^{2}+1.414\ln(z_{5}) - 23 = 0

which is of the form

\displaystyle f(z1, z2, z3, z4, z5) = 0

Once in this form (the standard form), solving the original equation becomes a matter of finding the z_{1}, z_{2}, z_{3}, z_{4} and z_{5} values for which f equals 0, and well-developed fields of mathematics and computer science provide several root-finding techniques for solving such a problem. Note that since the theme of this article is polynomial root-finding, further examples will focus on single-variable equations, specifically, polynomials.

Analytical Root-finding Techniques

Consider the following quadratic equation:

\displaystyle f(z) = z^{2} + z - 2

This equation is a second-degree polynomial–the highest power applied to the independent variable is 2. Consequently, this equation has two roots; an n-degree polynomial equation has n roots.

Because the equation is a simple polynomial, a first approach to finding its zeros might be to put the equation into the form (z - \alpha_{1})(z - \alpha_{2}) and make a few educated guesses; a person could quickly determine that \alpha_{1} = 1 (a strictly real result) is one root of this equation. This root could then be divided out of the original polynomial

\displaystyle \frac{(z^{2} + z -2)}{(z-1)} = (z+2)

to yield the second root: \alpha_{2} = -2.

In this simple example, the two roots were found easily; the second root was found immediately after the first one was known. However, in most real-world problems the roots of a quadratic are not found so easily. Furthermore, in problems where the original polynomial is of degree greater than two, when one root is found, the other roots do not follow immediately — more sophisticated techniques are required.

Among the techniques available are analytical ones, in other words, techniques that yield explicit, algebraic results. For example, in the case of a quadratic equation, explicit expressions for the two roots are possible via the Quadratic Formula: once a quadratic equation is arranged into the standard form

\displaystyle \bold{a}z^{2} + \bold{b}z + \bold{c} = 0

(a, b and c are known constants), the two roots are found via the Quadratic Formula:

\displaystyle \alpha = \frac{{ - b \pm \sqrt {b^2 - 4ac} }}{{2a}}

The quantity within the square root sign (b^2 - 4ac) is called the discriminant and determines whether the solution is complex or strictly real. If the discriminant is less than 0, the numerator contains an imaginary component and the roots are, therefore, complex. If the discriminant is greater than or equal to 0, the solutions are real. In fact, if the discriminant is 0, the two roots are real and equal:

\displaystyle \alpha_{1} = \frac{{ - b + \sqrt {0} }}{{2a}} = \frac{-b}{2a}

\displaystyle \alpha_{2} = \frac{{ - b - \sqrt {0} }}{{2a}} = \frac{-b}{2a}

In addition to the Quadratic Formula for quadratic equations, analogous formulae exist for polynomials of degree three and four. However, for polynomials of degree five and higher, analytical solutions are not possible (except in special cases), only numerical solutions are possible.

Numerical Root-finding Techniques

For the multitude of real-world problems that are not amenable to an analytical solution, numerical root-finding is a well-developed field offering a wide variety of tools from which to choose; many computer programs written for the computation of roots are based on algorithms that make them applicable to the computation of roots of functions in addition to polynomials, and virtually all of them employ an iterative approach that terminates once the desired degree of tolerance has been achieved. For example, the bisection method is a simple and reliable method for computing roots of a function when they are known ahead of time to be real only. The algorithm starts with the assumption that a zero is somewhere on a user-supplied interval [a,b]. In other words, the function value f(a) is of the opposite sign of f(b) – in going from a to b, f either goes from being positive to being negative, or it goes from being negative to being positive. A point, m, in the middle of the interval [a,b] is then selected. The function is evaluated at point m: f(m). If the sign of f(m) is the same as the sign of f(a), the desired zero is not in the interval [a,m]. In this case, the interval is cut in half and the new interval becomes [m,b]. On the other hand, if the sign of f(m) is not the same as the sign of f(a), the desired zero is in the interval [a,m]. In this case, the interval is cut in half and the new interval becomes [a,m]. This process is repeated until either the interval converges to a desired tolerance or a value of f is found that is within an acceptable tolerance, in which case the value of z at this point is the value returned by the program as the root.

In addition to the bisection method, more elegant–and faster converging–techniques are available: the False Position method, Brent’s method, the Newton-Raphson method, the secant method, and others. The Newton-Raphson method uses the function and its derivative to quickly converge to a solution. This method is good for situations in which the derivatives are either known or can be calculated with a low computational cost. In other situations, the cost of computing derivatives may be too high. The secant method, by comparison, does not require the use of derivatives, but the calculation of each iterate requires the use of the two previous iterates and does not converge to a solution as quickly as the Newton-Raphson method. In some situations, this method, too, may be deemed impractical. Usually, when some thought is given to a problem, a particular technique can be applied to it that is more appropriate than another. In fact, for some problems, one technique may fail to find a solution at all, whereas another technique will succeed. For other problems, several techniques may, indeed, be able to solve the problem and the numerical analyst may select the one that is simply more computationally efficient.

The Jenkins-Traub Algorithm

The field of numerical root-finding is so well-developed that the use of a good-quality numerical program is recommended — when numerical results are sought — even when an analytical solution to a problem is known. In other words, even though an analytical solution to, say, the quartic equation exists, writing a computer program that simply implements the textbook formula is not recommended; computer round-off errors often render results of such programs meaningless. The use of a robust numerical program, based upon sound theory and an excellent algorithm, and coded to thoroughly deal with computer round-off errors, is the recommended action. The Jenkins-Traub Algorithm is such an algorithm; it is a three-stage, extremely effective, globally convergent algorithm designed specifically for computing the roots of polynomials.

Stage One is the “No Shift” stage; the main purpose of this stage is to accentuate the smaller zeros. The search for a zero is started by taking an initial guess of 0 for a fixed number, M, of iterations (M is usually assigned the value 5 on the basis of numerical experience(1)).

Stage Two is the “Fixed Shift” stage, the purpose of this stage is to separate zeros of equal or almost equal magnitude. As a starting point in this stage, the following value is used:

\displaystyle s = \beta e^{i\theta}

\beta is a lower bound on the magnitudes of the probable zeros in the cluster. \theta could be taken at random, since the cluster could be anywhere in the complex plane; however, in practice \theta is usually initialized to 49°, putting s near the middle of the first quadrant of the complex plane. After a certain number of iterations, if s does not converge to a root, s is assigned a new value by increasing \theta by 94°. Repeated attempts would have the search for a root start with points in all four quadrants of the complex plane until the search is returned to the first quadrant. Should the search, indeed, return to the first quadrant, successive cycles start at points 16° away from the starting point of the preceding cycle.

Stage Three is the “Variable Shift” stage, which is terminated when the computed value of the polynomial at a possible zero is less than or equal to a specified bound.

In addition to the three fundamental stages around which it was developed, the Jenkins-Traub Algorithm incorporates several other techniques for making it as effective as possible. One of those techniques is deflation of the polynomial by synthetic division each time a root is found. Consider the following monic polynomial:

\displaystyle z^n + \bold{a_{n-1}}z^{n-1}+ \bold{a_{n-2}}z^{n-2}+\bold{a_{n-3}}z^{n-3} + \ldots + \bold{a_{1}}z + \bold{a_0} = 0

The \bold{a_i} are known constants and, in general, are complex. Now say the root z_1 (z_1 = x + iy) has been found. Synthetic division would be employed to divide that root out of the original polynomial:

\displaystyle \frac{z^n + \bold{a_{n-1}}z^{n-1}+ \bold{a_{n-2}}z^{n-2}+\bold{a_{n-3}}z^{n-3} + \ldots + \bold{a_0}}{z-z_1}

to yield

\displaystyle \bold{b_{n-1}}z^{n-1}+ \bold{b_{n-2}}z^{n-2}+\bold{a_{b-3}}z^{n-3} + \ldots + \bold{b_0}

The \bold{b_i} are new constants. The root-finding process is then repeated on this new–simpler–polynomial. As each root is found, the polynomial becomes successively simpler and each successive iteration of the algorithm involves, in general, fewer computations.

For polynomials whose coefficients are real only, when a complex root is found (z_1 = x + iy), an additional benefit arises: that root’s complex conjugate is also a root (z_2 = x - iy). In other words, two roots are computed — and the polynomial can be deflated by two degrees — in a single iteration of the algorithm. Furthermore, this deflation involves real only, rather than complex, operations; the product of these two roots is a real quadratic equation:

\displaystyle (z-z_1)(z - z_2) = z^2 - 2xz + (x^2 + y^2)

In this case, synthetic division is employed on the polynomial as follows:

\displaystyle \frac{z^n + \bold{a_{n-1}}z^{n-1}+ \bold{a_{n-2}}z^{n-2}+\bold{a_{n-3}}z^{n-3} + \ldots + \bold{a_0}}{z^2 - 2xz + (x^2 + y^2)}

In fact, for computing roots of polynomials which have only real coefficients, a modified version of the Jenkins-Traub Algorithm has been written which incorporates several features that take advantage of the characteristics of real-only polynomials to yield significant decreases in program execution time; “If the complex and real algorithms are applied to the same real polynomial, the real algorithm is about four times as fast.”(4)

The Jenkins-Traub Algorithm not only deflates the polynomial as roots are computed, it computes the roots roughly in order of increasing magnitude. This approach is taken because deflation of a polynomial can be unstable unless done by factors in order of increasing magnitude. By starting the search for roots with an initial guess of 0, as is done in Stage One, roots are, indeed, computed roughly in order of increasing magnitude, the factors by which the polynomial is successively deflated are roughly in order of increasing magnitude — the deflation of the polynomial is made quite stable.

As a quick check of the effectiveness of the Jenkins-Traub Algorithm, consider a contrived numerical example:

\displaystyle f(z) = (z - 1.001)(z - 0.998)(z - 1.00002)(z - 0.99999)

The roots of this polynomial are very close and should test how effectively the algorithm discerns near-multiple roots:

1.001
0.998
1.00002
0.99999

In addition, the polynomial is real-only so that readers may confirm the results with a free, ready-to-use, online Javascript program (URL given below).

Expanded, the test polynomial becomes

\displaystyle \displaystyle z^{4} - 3.99901z^{3} + 5.9970279898z^{2} - 3.9970259795802z + 0.9990079897802004

and the solutions computed by the Javascript root-finder follow:

0.9980051486669109
0.9998618452719573
1.0001539207016332
1.0009890853594987

The program, in fact, did a very good job:

  1. no errors were returned,
  2. no erroneous imaginary components were found,
  3. four distinct roots were computed, and
  4. the computed roots are extremely close to the exact roots (round-off errors are unavoidable).

Conclusion

Since its introduction in 1969 (in Michael Jenkins’ PhD thesis) the Jenkins-Traub Algorithm has been rigorously tested and has proven itself as an excellent algorithm for the computation of the roots of polynomials. It is extremely robust, is globally convergent for any distribution of zeros, and converges at a better than quadratic rate. Attesting to its high quality is the fact that source code for programs (written in FORTRAN) based on the algorithm have been posted on the NETLIB site, a repository of high-quality numerical programs; the complex coefficient version is posted as CPOLY (http://www.netlib.org/toms/419) and a real coefficient only version has been posted as RPOLY (http://www.netlib.org/toms/493). Translations into C++ and Javascript are also available; for example, Polynomial Root-finder (http://www.akiti.ca/PolyRootRe.html) is a Javascript translation of RPOLY and has been successfully tested with data published in Jenkins’ PhD thesis. In addition, the algorithm is incorporated in commercial software products (i.e. Mathematica). The proven effectiveness of the Jenkins-Traub Algorithm and its wide-spread distribution make it an extremely important tool in the toolbox of numerical analysts.

References

  1. “A First Course in Numerical Analysis: Second Edition”
    Anthony Ralston and Philip Rabinowitz
    McGraw-Hill Book Company
    New York
    1978
  2. “Numerical Methods and Software”
    David Kahaner, Cleve Moler, Stephen Nash
    Prentice-Hall, Inc.
    Englewood Cliffs, New Jersey
    1989
  3. “Numerical Recipes. The Art of Scientific Computing.”
    William H. Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling
    Cambridge University Press
    New York
    1986
    (3rd edition now available.)
  4. Wikipedia
    Jenkins-Traub Algorithm for Polynomial Zeros
  5. Jenkins, M.A. Three-stage variable-shift iterations for the solution of polynomial equations with a posteriori error bounds for the zeros. Diss., Rep. CS 138, Comput. Sci.
    Dep., Stanford U., Stanford, Cal., 1969.

This article was written by David Binner. You can visit his math site at http://akiti.ca. If you’d like to write for Math-Blog.com, please email us at submissions@math-blog.com.