diff options
Diffstat (limited to 'gsl-1.9/doc/poly.texi')
| -rw-r--r-- | gsl-1.9/doc/poly.texi | 302 |
1 files changed, 302 insertions, 0 deletions
diff --git a/gsl-1.9/doc/poly.texi b/gsl-1.9/doc/poly.texi new file mode 100644 index 0000000..531f77a --- /dev/null +++ b/gsl-1.9/doc/poly.texi @@ -0,0 +1,302 @@ +@cindex polynomials, roots of + +This chapter describes functions for evaluating and solving polynomials. +There are routines for finding real and complex roots of quadratic and +cubic equations using analytic methods. An iterative polynomial solver +is also available for finding the roots of general polynomials with real +coefficients (of any order). The functions are declared in the header +file @code{gsl_poly.h}. + +@menu +* Polynomial Evaluation:: +* Divided Difference Representation of Polynomials:: +* Quadratic Equations:: +* Cubic Equations:: +* General Polynomial Equations:: +* Roots of Polynomials Examples:: +* Roots of Polynomials References and Further Reading:: +@end menu + +@node Polynomial Evaluation +@section Polynomial Evaluation +@cindex polynomial evaluation +@cindex evaluation of polynomials + +@deftypefun double gsl_poly_eval (const double @var{c}[], const int @var{len}, const double @var{x}) +This function evaluates the polynomial +@c{$c[0] + c[1] x + c[2] x^2 + \dots + c[len-1] x^{len-1}$} +@math{c[0] + c[1] x + c[2] x^2 + \dots + c[len-1] x^@{len-1@}} using +Horner's method for stability. The function is inlined when possible. +@end deftypefun + +@node Divided Difference Representation of Polynomials +@section Divided Difference Representation of Polynomials +@cindex divided differences, polynomials +@cindex evaluation of polynomials, in divided difference form + +The functions described here manipulate polynomials stored in Newton's +divided-difference representation. The use of divided-differences is +described in Abramowitz & Stegun sections 25.1.4 and 25.2.26. + +@deftypefun int gsl_poly_dd_init (double @var{dd}[], const double @var{xa}[], const double @var{ya}[], size_t @var{size}) +This function computes a divided-difference representation of the +interpolating polynomial for the points (@var{xa}, @var{ya}) stored in +the arrays @var{xa} and @var{ya} of length @var{size}. On output the +divided-differences of (@var{xa},@var{ya}) are stored in the array +@var{dd}, also of length @var{size}. +@end deftypefun + +@deftypefun double gsl_poly_dd_eval (const double @var{dd}[], const double @var{xa}[], const size_t @var{size}, const double @var{x}) +This function evaluates the polynomial stored in divided-difference form +in the arrays @var{dd} and @var{xa} of length @var{size} at the point +@var{x}. +@end deftypefun + +@deftypefun int gsl_poly_dd_taylor (double @var{c}[], double @var{xp}, const double @var{dd}[], const double @var{xa}[], size_t @var{size}, double @var{w}[]) +This function converts the divided-difference representation of a +polynomial to a Taylor expansion. The divided-difference representation +is supplied in the arrays @var{dd} and @var{xa} of length @var{size}. +On output the Taylor coefficients of the polynomial expanded about the +point @var{xp} are stored in the array @var{c} also of length +@var{size}. A workspace of length @var{size} must be provided in the +array @var{w}. +@end deftypefun + +@node Quadratic Equations +@section Quadratic Equations +@cindex quadratic equation, solving + +@deftypefun int gsl_poly_solve_quadratic (double @var{a}, double @var{b}, double @var{c}, double * @var{x0}, double * @var{x1}) +This function finds the real roots of the quadratic equation, +@tex +\beforedisplay +$$ +a x^2 + b x + c = 0 +$$ +\afterdisplay +@end tex +@ifinfo + +@example +a x^2 + b x + c = 0 +@end example + +@end ifinfo +@noindent +The number of real roots (either zero, one or two) is returned, and +their locations are stored in @var{x0} and @var{x1}. If no real roots +are found then @var{x0} and @var{x1} are not modified. If one real root +is found (i.e. if @math{a=0}) then it is stored in @var{x0}. When two +real roots are found they are stored in @var{x0} and @var{x1} in +ascending order. The case of coincident roots is not considered +special. For example @math{(x-1)^2=0} will have two roots, which happen +to have exactly equal values. + +The number of roots found depends on the sign of the discriminant +@math{b^2 - 4 a c}. This will be subject to rounding and cancellation +errors when computed in double precision, and will also be subject to +errors if the coefficients of the polynomial are inexact. These errors +may cause a discrete change in the number of roots. However, for +polynomials with small integer coefficients the discriminant can always +be computed exactly. + +@end deftypefun + +@deftypefun int gsl_poly_complex_solve_quadratic (double @var{a}, double @var{b}, double @var{c}, gsl_complex * @var{z0}, gsl_complex * @var{z1}) + +This function finds the complex roots of the quadratic equation, +@tex +\beforedisplay +$$ +a z^2 + b z + c = 0 +$$ +\afterdisplay +@end tex +@ifinfo + +@example +a z^2 + b z + c = 0 +@end example + +@end ifinfo +@noindent +The number of complex roots is returned (either one or two) and the +locations of the roots are stored in @var{z0} and @var{z1}. The roots +are returned in ascending order, sorted first by their real components +and then by their imaginary components. If only one real root is found +(i.e. if @math{a=0}) then it is stored in @var{z0}. + +@end deftypefun + + +@node Cubic Equations +@section Cubic Equations +@cindex cubic equation, solving + +@deftypefun int gsl_poly_solve_cubic (double @var{a}, double @var{b}, double @var{c}, double * @var{x0}, double * @var{x1}, double * @var{x2}) + +This function finds the real roots of the cubic equation, +@tex +\beforedisplay +$$ +x^3 + a x^2 + b x + c = 0 +$$ +\afterdisplay +@end tex +@ifinfo + +@example +x^3 + a x^2 + b x + c = 0 +@end example + +@end ifinfo +@noindent +with a leading coefficient of unity. The number of real roots (either +one or three) is returned, and their locations are stored in @var{x0}, +@var{x1} and @var{x2}. If one real root is found then only @var{x0} is +modified. When three real roots are found they are stored in @var{x0}, +@var{x1} and @var{x2} in ascending order. The case of coincident roots +is not considered special. For example, the equation @math{(x-1)^3=0} +will have three roots with exactly equal values. + +@end deftypefun + +@deftypefun int gsl_poly_complex_solve_cubic (double @var{a}, double @var{b}, double @var{c}, gsl_complex * @var{z0}, gsl_complex * @var{z1}, gsl_complex * @var{z2}) + +This function finds the complex roots of the cubic equation, +@tex +\beforedisplay +$$ +z^3 + a z^2 + b z + c = 0 +$$ +\afterdisplay +@end tex +@ifinfo + +@example +z^3 + a z^2 + b z + c = 0 +@end example + +@end ifinfo +@noindent +The number of complex roots is returned (always three) and the locations +of the roots are stored in @var{z0}, @var{z1} and @var{z2}. The roots +are returned in ascending order, sorted first by their real components +and then by their imaginary components. + +@end deftypefun + + +@node General Polynomial Equations +@section General Polynomial Equations +@cindex general polynomial equations, solving + +The roots of polynomial equations cannot be found analytically beyond +the special cases of the quadratic, cubic and quartic equation. The +algorithm described in this section uses an iterative method to find the +approximate locations of roots of higher order polynomials. + +@deftypefun {gsl_poly_complex_workspace *} gsl_poly_complex_workspace_alloc (size_t @var{n}) +This function allocates space for a @code{gsl_poly_complex_workspace} +struct and a workspace suitable for solving a polynomial with @var{n} +coefficients using the routine @code{gsl_poly_complex_solve}. + +The function returns a pointer to the newly allocated +@code{gsl_poly_complex_workspace} if no errors were detected, and a null +pointer in the case of error. +@end deftypefun + +@deftypefun void gsl_poly_complex_workspace_free (gsl_poly_complex_workspace * @var{w}) +This function frees all the memory associated with the workspace +@var{w}. +@end deftypefun + +@deftypefun int gsl_poly_complex_solve (const double * @var{a}, size_t @var{n}, gsl_poly_complex_workspace * @var{w}, gsl_complex_packed_ptr @var{z}) +This function computes the roots of the general polynomial +@c{$P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_{n-1} x^{n-1}$} +@math{P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_@{n-1@} x^@{n-1@}} using +balanced-QR reduction of the companion matrix. The parameter @var{n} +specifies the length of the coefficient array. The coefficient of the +highest order term must be non-zero. The function requires a workspace +@var{w} of the appropriate size. The @math{n-1} roots are returned in +the packed complex array @var{z} of length @math{2(n-1)}, alternating +real and imaginary parts. + +The function returns @code{GSL_SUCCESS} if all the roots are found. If +the QR reduction does not converge, the error handler is invoked with +an error code of @code{GSL_EFAILED}. Note that due to finite precision, +roots of higher multiplicity are returned as a cluster of simple roots +with reduced accuracy. The solution of polynomials with higher-order +roots requires specialized algorithms that take the multiplicity +structure into account (see e.g. Z. Zeng, Algorithm 835, ACM +Transactions on Mathematical Software, Volume 30, Issue 2 (2004), pp +218--236). +@end deftypefun + +@node Roots of Polynomials Examples +@section Examples + +To demonstrate the use of the general polynomial solver we will take the +polynomial @math{P(x) = x^5 - 1} which has the following roots, +@tex +\beforedisplay +$$ +1, e^{2\pi i /5}, e^{4\pi i /5}, e^{6\pi i /5}, e^{8\pi i /5} +$$ +\afterdisplay +@end tex +@ifinfo + +@example +1, e^@{2\pi i /5@}, e^@{4\pi i /5@}, e^@{6\pi i /5@}, e^@{8\pi i /5@} +@end example + +@end ifinfo +@noindent +The following program will find these roots. + +@example +@verbatiminclude examples/polyroots.c +@end example + +@noindent +The output of the program is, + +@example +$ ./a.out +@verbatiminclude examples/polyroots.out +@end example + +@noindent +which agrees with the analytic result, @math{z_n = \exp(2 \pi n i/5)}. + +@node Roots of Polynomials References and Further Reading +@section References and Further Reading + +The balanced-QR method and its error analysis are described in the +following papers, + +@itemize @asis +@item +R.S. Martin, G. Peters and J.H. Wilkinson, ``The QR Algorithm for Real +Hessenberg Matrices'', @cite{Numerische Mathematik}, 14 (1970), 219--231. + +@item +B.N. Parlett and C. Reinsch, ``Balancing a Matrix for Calculation of +Eigenvalues and Eigenvectors'', @cite{Numerische Mathematik}, 13 (1969), +293--304. + +@item +A. Edelman and H. Murakami, ``Polynomial roots from companion matrix +eigenvalues'', @cite{Mathematics of Computation}, Vol.@: 64, No.@: 210 +(1995), 763--776. +@end itemize + +@noindent +The formulas for divided differences are given in Abramowitz and Stegun, + +@itemize @asis +@item +Abramowitz and Stegun, @cite{Handbook of Mathematical Functions}, +Sections 25.1.4 and 25.2.26. +@end itemize |