This is gsl-ref.info, produced by makeinfo version 4.0 from gsl-ref.texi. INFO-DIR-SECTION Scientific software START-INFO-DIR-ENTRY * gsl-ref: (gsl-ref). GNU Scientific Library - Reference END-INFO-DIR-ENTRY This file documents the GNU Scientific Library. Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001 The GSL Team. Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.1 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A copy of the license is included in the section entitled "GNU Free Documentation License".  File: gsl-ref.info, Node: Traveling Salesman Problem, Prev: Trivial example, Up: Examples with Simulated Annealing Traveling Salesman Problem -------------------------- The TSP ("Traveling Salesman Problem") is the classic combinatorial optimization problem. I have provided a very simple version of it, based on the coordinates of twelve cities in the southwestern United States. This should maybe be called the "Flying Salesman Problem", since I am using the great-circle distance between cities, rather than the driving distance. Also: I assume the earth is a sphere, so I don't use geoid distances. The `gsl_siman_solve()' routine finds a route which is 3490.62 Kilometers long; this is confirmed by an exhaustive search of all possible routes with the same initial city. The full code can be found in `siman/siman_tsp.c', but I include here some plots generated with in the following way: ./siman_tsp > tsp.output grep -v "^#" tsp.output | xyplot -xyil -d "x................y" -lx "generation" -ly "distance" -lt "TSP -- 12 southwest cities" | xyps -d > 12-cities.eps grep initial_city_coord tsp.output | awk '{print $2, $3, $4, $5}' | xyplot -xyil -lb0 -cs 0.8 -lx "longitude (- means west)" -ly "latitude" -lt "TSP -- initial-order" | xyps -d > initial-route.eps grep final_city_coord tsp.output | awk '{print $2, $3, $4, $5}' | xyplot -xyil -lb0 -cs 0.8 -lx "longitude (- means west)" -ly "latitude" -lt "TSP -- final-order" | xyps -d > final-route.eps This is the output showing the initial order of the cities; longitude is negative, since it is west and I want the plot to look like a map. # initial coordinates of cities (longitude and latitude) ###initial_city_coord: -105.95 35.68 Santa Fe ###initial_city_coord: -112.07 33.54 Phoenix ###initial_city_coord: -106.62 35.12 Albuquerque ###initial_city_coord: -103.2 34.41 Clovis ###initial_city_coord: -107.87 37.29 Durango ###initial_city_coord: -96.77 32.79 Dallas ###initial_city_coord: -105.92 35.77 Tesuque ###initial_city_coord: -107.84 35.15 Grants ###initial_city_coord: -106.28 35.89 Los Alamos ###initial_city_coord: -106.76 32.34 Las Cruces ###initial_city_coord: -108.58 37.35 Cortez ###initial_city_coord: -108.74 35.52 Gallup ###initial_city_coord: -105.95 35.68 Santa Fe The optimal route turns out to be: # final coordinates of cities (longitude and latitude) ###final_city_coord: -105.95 35.68 Santa Fe ###final_city_coord: -106.28 35.89 Los Alamos ###final_city_coord: -106.62 35.12 Albuquerque ###final_city_coord: -107.84 35.15 Grants ###final_city_coord: -107.87 37.29 Durango ###final_city_coord: -108.58 37.35 Cortez ###final_city_coord: -108.74 35.52 Gallup ###final_city_coord: -112.07 33.54 Phoenix ###final_city_coord: -106.76 32.34 Las Cruces ###final_city_coord: -96.77 32.79 Dallas ###final_city_coord: -103.2 34.41 Clovis ###final_city_coord: -105.92 35.77 Tesuque ###final_city_coord: -105.95 35.68 Santa Fe Here's a plot of the cost function (energy) versus generation (point in the calculation at which a new temperature is set) for this problem:  File: gsl-ref.info, Node: Ordinary Differential Equations, Next: Interpolation, Prev: Simulated Annealing, Up: Top Ordinary Differential Equations ******************************* This chapter describes functions for solving ordinary differential equation (ODE) initial value problems. The library provides a variety of low-level methods, such as Runge-Kutta and Bulirsch-Stoer routines, and higher-level components for adaptive step-size control. The components can be combined by the user to achieve the desired solution, with full access to any intermediate steps. These functions are declared in the header file `gsl_odeiv.h'. * Menu: * Defining the ODE System:: * Stepping Functions:: * Adaptive Step-size Control:: * Evolution:: * ODE Example programs:: * ODE References and Further Reading::  File: gsl-ref.info, Node: Defining the ODE System, Next: Stepping Functions, Up: Ordinary Differential Equations Defining the ODE System ======================= The routines solve the general n-dimensional first-order system, dy_i(t)/dt = f_i(t, y_1(t), ..., y_n(t)) for i = 1, \dots, n. The stepping functions rely on the vector of derivatives f_i and the Jacobian matrix, J_{ij} = df_i(t,y(t)) / dy_j. A system of equations is defined using the `gsl_odeiv_system' datatype. - Data Type: gsl_odeiv_system This data type defines a general ODE system with arbitrary parameters. `int (* FUNCTION) (double t, const double y[], double dydt[], void * params)' This function should store the vector elements f_i(t,y,params) in the array DYDT, for arguments (T,Y) and parameters PARAMS `int (* JACOBIAN) (double t, const double y[], double * dfdy, double dfdt[], void * params);' This function should store the vector elements df_i(t,y,params)/dt in the array DFDT and the Jacobian matrix J_{ij} in the the array DFDY regarded as a row-ordered matrix `J(i,j) = dfdy[i * dim + j]' where `dim' is the dimension of the system. Some of the simpler solver algorithms do not make use of the Jacobian matrix, so it is not always strictly necessary to provide it (this element of the struct can be replace by a null pointer). However, it is useful to provide the Jacobian to allow the solver algorithms to be interchanged - the best algorithms make use of the Jacobian. `size_t DIMENSION;' This is the dimension of the system of equations `void * PARAMS' This is a pointer to the arbitrary parameters of the system.  File: gsl-ref.info, Node: Stepping Functions, Next: Adaptive Step-size Control, Prev: Defining the ODE System, Up: Ordinary Differential Equations Stepping Functions ================== The lowest level components are the "stepping functions" which advance a solution from time t to t+h for a fixed step-size h and estimate the resulting local error. - Function: gsl_odeiv_step * gsl_odeiv_step_alloc (const gsl_odeiv_step_type * T, size_t DIM) This function returns a pointer to a newly allocated instance of a stepping function of type T for a system of DIM dimensions. - Function: int gsl_odeiv_step_reset (gsl_odeiv_step * S) This function resets the stepping function S. It should be used whenever the next use of S will not be a continuation of a previous step. - Function: void gsl_odeiv_step_free (gsl_odeiv_step * S) This function frees all the memory associated with the stepping function S. - Function: const char * gsl_odeiv_step_name (const gsl_odeiv_step * S) This function returns a pointer to the name of the stepping function. For example, printf("step method is '%s'\n", gsl_odeiv_step_name (s)); would print something like `step method is 'rk4''. - Function: unsigned int gsl_odeiv_step_order (const gsl_odeiv_step * S) This function returns the order of the stepping function on the previous step. This order can vary if the stepping function itself is adaptive. - Function: int gsl_odeiv_step_apply (gsl_odeiv_step * S, double T, double H, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv_system * DYDT) This function applies the stepping function S to the system of equations defined by DYDT, using the step size H to advance the system from time T and state Y to time T+H. The new state of the system is stored in Y on output, with an estimate of the absolute error in each component stored in YERR. If the argument DYDT_IN is not null it should point an array containing the derivatives for the system at time T on input. This is optional as the derivatives will be computed internally if they are not provided, but allows the reuse of existing derivative information. On output the new derivatives of the system at time T+H will be stored in DYDT_OUT if it is not null. The following algorithms are available, - Step Type: gsl_odeiv_step_rk2 Embedded 2nd order Runge-Kutta with 3rd order error estimate. - Step Type: gsl_odeiv_step_rk4 4th order (classical) Runge-Kutta. - Step Type: gsl_odeiv_step_rkf45 Embedded 4th order Runge-Kutta-Fehlberg method with 5th order error estimate. This method is a good general-purpose integrator. - Step Type: gsl_odeiv_step_rkck Embedded 4th order Runge-Kutta Cash-Karp method with 5th order error estimate. - Step Type: gsl_odeiv_step_rk8pd Embedded 8th order Runge-Kutta Prince-Dormand method with 9th order error estimate. - Step Type: gsl_odeiv_step_rk2imp Implicit 2nd order Runge-Kutta at Gaussian points - Step Type: gsl_odeiv_step_rk4imp Implicit 4th order Runge-Kutta at Gaussian points - Step Type: gsl_odeiv_step_bsimp Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm requires the Jacobian. - Step Type: gsl_odeiv_step_gear1 M=1 implicit Gear method - Step Type: gsl_odeiv_step_gear2 M=2 implicit Gear method  File: gsl-ref.info, Node: Adaptive Step-size Control, Next: Evolution, Prev: Stepping Functions, Up: Ordinary Differential Equations Adaptive Step-size Control ========================== The control function examines the proposed change to the solution and its error estimate produced by a stepping function and attempts to determine the optimal step-size for a user-specified level of error. - Function: gsl_odeiv_control * gsl_odeiv_control_standard_new (double EPS_ABS, double EPS_REL, double A_Y, double A_DYDT) The standard control object is a four parameter heuristic based on absolute and relative errors EPS_ABS and EPS_REL, and scaling factors A_Y and A_DYDT for the system state y(t) and derivatives y'(t) respectively. The step-size adjustment procedure for this method begins by computing the desired error level D_i for each component, D_i = eps_abs + eps_rel * (a_y |y_i| + a_dydt h |y'_i|) and comparing it with the observed error E_i = |yerr_i|. If the observed error E exceeds the desired error level D by more than 10% for any component then the method reduces the step-size by an appropriate factor, h_new = h_old * S * (D/E)^(1/q) where q is the consistency order of method (e.g. q=4 for 4(5) embedded RK), and S is a safety factor of 0.9. The ratio D/E is taken to be the maximum of the ratios D_i/E_i. If the observed error E is less than 50% of the desired error level D for the maximum ratio D_i/E_i then the algorithm takes the opportunity to increase the step-size to bring the error in line with the desired level, h_new = h_old * S * (E/D)^(1/(q+1)) This encompasses all the standard error scaling methods. - Function: gsl_odeiv_control * gsl_odeiv_control_y_new (double EPS_ABS, double EPS_REL) This function creates a new control object which will keep the local error on each step within an absolute error of EPS_ABS and relative error of EPS_REL with respect to the solution y_i(t). This is equivalent to the standard control object with A_Y=1 and A_DYDT=0. - Function: gsl_odeiv_control * gsl_odeiv_control_yp_new (double EPS_ABS, double EPS_REL) This function creates a new control object which will keep the local error on each step within an absolute error of EPS_ABS and relative error of EPS_REL with respect to the derivatives of the solution y'_i(t) . This is equivalent to the standard control object with A_Y=0 and A_DYDT=1. - Function: gsl_odeiv_control * gsl_odeiv_control_alloc (const gsl_odeiv_control_type * T) This function returns a pointer to a newly allocated instance of a control function of type T. This function is only needed for defining new types of control functions. For most purposes the standard control functions described above should be sufficient. - Function: int gsl_odeiv_control_init (gsl_odeiv_control * C, double EPS_ABS, double EPS_REL, double A_Y, double A_DYDT) This function initializes the control function C with the parameters EPS_ABS (absolute error), EPS_REL (relative error), A_Y (scaling factor for y) and A_DYDT (scaling factor for derivatives). - Function: void gsl_odeiv_control_free (gsl_odeiv_control * C) This function frees all the memory associated with the control function C. - Function: int gsl_odeiv_control_hadjust (gsl_odeiv_control * C, gsl_odeiv_step * S, const double y0[], const double yerr[], const double dydt[], double * H) This function adjusts the step-size H using the control function C, and the current values of Y, YERR and DYDT. The stepping function STEP is also needed to determine the order of the method. If the error in the y-values YERR is found to be too large then the step-size H is reduced and the function returns `GSL_ODEIV_HADJ_DEC'. If the error is sufficiently small then H may be increased and `GSL_ODEIV_HADJ_INC' is returned. The function returns `GSL_ODEIV_HADJ_NIL' if the step-size is unchanged. The goal of the function is to estimate the largest step-size which satisfies the user-specified accuracy requirements for the current point. - Function: const char * gsl_odeiv_control_name (const gsl_odeiv_control * C) This function returns a pointer to the name of the control function. For example, printf("control method is '%s'\n", gsl_odeiv_control_name (c)); would print something like `control method is 'standard''  File: gsl-ref.info, Node: Evolution, Next: ODE Example programs, Prev: Adaptive Step-size Control, Up: Ordinary Differential Equations Evolution ========= The highest level of the system is the evolution function which combines the results of a stepping function and control function to reliably advance the solution forward over an interval (t_0, t_1). If the control function signals that the step-size should be decreased the evolution function backs out of the current step and tries the proposed smaller step-size. This is process is continued until an acceptable step-size is found. - Function: gsl_odeiv_evolve * gsl_odeiv_evolve_alloc (size_t DIM) This function returns a pointer to a newly allocated instance of an evolution function for a system of DIM dimensions. - Function: int gsl_odeiv_evolve_apply (gsl_odeiv_evolve * E, gsl_odeiv_control * CON, gsl_odeiv_step * STEP, const gsl_odeiv_system * DYDT, double * T, double T1, double * H, double y[]) This function advances the system (E, DYDT) from time T and position Y using the stepping function STEP. The new time and position are stored in T and Y on output. The initial step-size is taken as H, but this will be modified using the control function C to achieve the appropriate error bound if necessary. The routine may make several calls to STEP in order to determine the optimum step-size. If the step-size has been changed the value of H will be modified on output. The maximum time T1 is guaranteed not to be exceeded by the time-step. On the final time-step the value of T will be set to T1 exactly. - Function: int gsl_odeiv_evolve_reset (gsl_odeiv_evolve * E) This function resets the evolution function E. It should be used whenever the next use of E will not be a continuation of a previous step. - Function: void gsl_odeiv_evolve_free (gsl_odeiv_evolve * E) This function frees all the memory associated with the evolution function E.  File: gsl-ref.info, Node: ODE Example programs, Next: ODE References and Further Reading, Prev: Evolution, Up: Ordinary Differential Equations Examples ======== The following program solves the second-order nonlinear Van der Pol oscillator equation, x''(t) + \mu x'(t) (x(t)^2 - 1) + x(t) = 0 This can be converted into a first order system suitable for use with the library by introducing a separate variable for the velocity, y = x'(t), x' = y y' = -x + \mu y (1-x^2) The program begins by defining functions for these derivatives and their Jacobian, #include #include #include #include int func (double t, const double y[], double f[], void *params) { double mu = *(double *)params; f[0] = y[1]; f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1); return GSL_SUCCESS; } int jac (double t, const double y[], double *dfdy, double dfdt[], void *params) { double mu = *(double *)params; gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2); gsl_matrix * m = &dfdy_mat.matrix; gsl_matrix_set (m, 0, 0, 0.0); gsl_matrix_set (m, 0, 1, 1.0); gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0); gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0)); dfdt[0] = 0.0; dfdt[1] = 0.0; return GSL_SUCCESS; } int main (void) { const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd; gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 2); gsl_odeiv_control * c = gsl_odeiv_control_y_new (1e-6, 0.0); gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (2); double mu = 10; gsl_odeiv_system sys = {func, jac, 2, &mu}; double t = 0.0, t1 = 100.0; double h = 1e-6; double y[2] = { 1.0, 0.0 }; while (t < t1) { int status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, t1, &h, y); if (status != GSL_SUCCESS) break; printf("%.5e %.5e %.5e\n", t, y[0], y[1]); } gsl_odeiv_evolve_free(e); gsl_odeiv_control_free(c); gsl_odeiv_step_free(s); return 0; } The main loop of the program evolves the solution from (y, y') = (1, 0) at t=0 to t=100. The step-size h is automatically adjusted by the controller to maintain an absolute accuracy of 10^{-6} in the function values Y. To obtain the values at regular intervals, rather than the variable spacings chosen by the control function, the main loop can be modified to advance the solution from one point to the next. For example, the following main loop prints the solution at the fixed points t = 0, 1, 2, \dots, 100, for (i = 1; i <= 100; i++) { double ti = i * t1 / 100.0; while (t < ti) { gsl_odeiv_evolve_apply (e, c, s, &sys, &t, ti, &h, y); } printf("%.5e %.5e %.5e\n", t, y[0], y[1]); } It is also possible to work with a non-adaptive integrator, using only the stepping function itself. The following program uses the `rk4' fourth-order Runge-Kutta stepping function with a fixed stepsize of 0.01, int main (void) { const gsl_odeiv_step_type * T = gsl_odeiv_step_rk4; gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 2); double mu = 10; gsl_odeiv_system sys = {func, jac, 2, &mu}; double t = 0.0, t1 = 100.0; double h = 1e-2; double y[2] = { 1.0, 0.0 }, y_err[2], dfdy[4], dydt_in[2], dydt_out[2]; /* initialise dydt_in */ GSL_ODEIV_JA_EVAL(&sys, t, y, dfdy, dydt_in); while (t < t1) { int status = gsl_odeiv_step_apply (s, t, h, y, y_err, dydt_in, dydt_out, &sys); if (status != GSL_SUCCESS) break; dydt_in[0] = dydt_out[0]; dydt_in[1] = dydt_out[1]; t += h; printf("%.5e %.5e %.5e\n", t, y[0], y[1]); } gsl_odeiv_step_free(s); return 0; } The derivatives and jacobian must be initialised for the starting point t=0 before the first step is taken. Subsequent steps use the output derivatives DYDT_OUT as inputs to the next step by copying their values into DYDT_IN.  File: gsl-ref.info, Node: ODE References and Further Reading, Prev: ODE Example programs, Up: Ordinary Differential Equations References and Further Reading ============================== Many of the the basic Runge-Kutta formulas can be found in the Handbook of Mathematical Functions, Abramowitz & Stegun (eds.), `Handbook of Mathematical Functions', Section 25.5. The implicit Bulirsch-Stoer algorithm `bsimp' is described in the following paper, G. Bader and P. Deuflhard, "A Semi-Implicit Mid-Point Rule for Stiff Systems of Ordinary Differential Equations.", Numer. Math. 41, 373-398, 1983.  File: gsl-ref.info, Node: Interpolation, Next: Numerical Differentiation, Prev: Ordinary Differential Equations, Up: Top Interpolation ************* This chapter describes functions for performing interpolation. The library provides a variety of interpolation methods, including Cubic splines and Akima splines. The interpolation types are interchangeable, allowing different methods to be used without recompiling. Interpolations can be defined for both normal and periodic boundary conditions. Additional functions are available for computing derivatives and integrals of interpolating functions. The functions described in this section are declared in the header files `gsl_interp.h' and `gsl_spline.h'. * Menu: * Introduction to Interpolation:: * Interpolation Functions:: * Interpolation Types:: * Index Look-up and Acceleration:: * Evaluation of interpolating functions:: * Higher-level interface:: * Interpolation Example programs:: * Interpolation References and Further Reading::  File: gsl-ref.info, Node: Introduction to Interpolation, Next: Interpolation Functions, Up: Interpolation Introduction ============ Given a set of data points (x_1, y_1) \dots (x_n, y_n) the routines described in this section compute a continuous interpolating function y(x) such that y_i = y(x_i). The interpolation is piecewise smooth, and its behavior at the points is determined by the type of interpolation used.  File: gsl-ref.info, Node: Interpolation Functions, Next: Interpolation Types, Prev: Introduction to Interpolation, Up: Interpolation Interpolation Functions ======================= The interpolation function for a given dataset is stored in a `gsl_interp' object. These are created by the following functions. - Function: gsl_interp * gsl_interp_alloc (const gsl_interp_type * T, size_t SIZE) This function returns a pointer to a newly allocated interpolation object of type T for SIZE data-points. - Function: int gsl_interp_init (gsl_interp * INTERP, const double XA[], const double YA[], size_t SIZE) This function initializes the interpolation object INTERP for the data (XA,YA) where XA and YA are arrays of size SIZE. The interpolation object (`gsl_interp') does not save the data arrays XA and YA and only stores the static state computed from the data. The XA data array is always assumed to be strictly ordered; the behavior for other arrangements is not defined. - Function: void gsl_interp_free (gsl_interp * INTERP) This function frees the interpolation object INTERP.  File: gsl-ref.info, Node: Interpolation Types, Next: Index Look-up and Acceleration, Prev: Interpolation Functions, Up: Interpolation Interpolation Types =================== The interpolation library provides five interpolation types: - Interpolation Type: gsl_interp_linear Linear interpolation. This interpolation method does not require any additional memory. - Interpolation Type: gsl_interp_polynomial Polynomial interpolation. This method should only be used for interpolating small numbers of points because polynomial interpolation introduces large oscillations, even for well-behaved datasets. The number of terms in the interpolating polynomial is equal to the number of points. - Interpolation Type: gsl_interp_cspline Cubic spline with natural boundary conditions - Interpolation Type: gsl_interp_cspline_periodic Cubic spline with periodic boundary conditions - Interpolation Type: gsl_interp_akima Akima spline with natural boundary conditions - Interpolation Type: gsl_interp_akima_periodic Akima spline with periodic boundary conditions The following related functions are available, - Function: const char * gsl_interp_name (const gsl_interp * INTERP) This function returns the name of the interpolation type used by INTERP. For example, printf("interp uses '%s' interpolation\n", gsl_interp_name (interp)); would print something like, interp uses 'cspline' interpolation. - Function: unsigned int gsl_interp_min_size (const gsl_interp * INTERP) This function returns the minimum number of points required by the interpolation type of INTERP. For example, Akima spline interpolation requires a minimum of 5 points.  File: gsl-ref.info, Node: Index Look-up and Acceleration, Next: Evaluation of interpolating functions, Prev: Interpolation Types, Up: Interpolation Index Look-up and Acceleration ============================== The state of searches can be stored in a `gsl_interp_accel' object, which is a kind of iterator for interpolation lookups. It caches the previous value of an index lookup. When the subsequent interpolation point falls in the same interval its index value can be returned immediately. - Function: size_t gsl_interp_bsearch (const double x_array[], double X, size_t INDEX_LO, size_t INDEX_HI) This function returns the index i of the array X_ARRAY such that `x_array[i] <= x < x_array[i+1]'. The index is searched for in the range [INDEX_LO,INDEX_HI]. - Function: gsl_interp_accel * gsl_interp_accel_alloc (void) This function returns a pointer to an accelerator object, which is a kind of iterator for interpolation lookups. It tracks the state of lookups, thus allowing for application of various acceleration strategies. - Function: size_t gsl_interp_accel_find (gsl_interp_accel * A, const double x_array[], size_t SIZE, double X) This function performs a lookup action on the data array X_ARRAY of size SIZE, using the given accelerator A. This is how lookups are performed during evaluation of an interpolation. The function returns an index i such that `xarray[i] <= x < xarray[i+1]'. - Function: void gsl_interp_accel_free (gsl_interp_accel* A) This function frees the accelerator object A.  File: gsl-ref.info, Node: Evaluation of interpolating functions, Next: Higher-level interface, Prev: Index Look-up and Acceleration, Up: Interpolation Evaluation of interpolating functions ===================================== - Function: double gsl_interp_eval (const gsl_interp * INTERP, const double XA[], const double YA[], double X, gsl_interp_accel * A) - Function: int gsl_interp_eval_e (const gsl_interp * INTERP, const double XA[], const double YA[], double X, gsl_interp_accel * A, double * Y) These functions return the interpolated value of Y for a given point X, using the interpolation object INTERP, data arrays XA and YA and the accelerator A. - Function: double gsl_interp_eval_deriv (const gsl_interp * INTERP, const double XA[], const double YA[], double X, gsl_interp_accel * A) - Function: int gsl_interp_eval_deriv_e (const gsl_interp * INTERP, const double XA[], const double YA[], double X, gsl_interp_accel * A, double * D) These functions return the derivative D of an interpolated function for a given point X, using the interpolation object INTERP, data arrays XA and YA and the accelerator A. - Function: double gsl_interp_eval_deriv2 (const gsl_interp * INTERP, const double XA[], const double YA[], double X, gsl_interp_accel * A) - Function: int gsl_interp_eval_deriv2_e (const gsl_interp * INTERP, const double XA[], const double YA[], double X, gsl_interp_accel * A, double * D2) These functions return the second derivative D2 of an interpolated function for a given point X, using the interpolation object INTERP, data arrays XA and YA and the accelerator A. - Function: double gsl_interp_eval_integ (const gsl_interp * INTERP, const double XA[], const double YA[], double A, double Bdouble X, gsl_interp_accel * A) - Function: int gsl_interp_eval_integ_e (const gsl_interp * INTERP, const double XA[], const double YA[], , double A, double B, gsl_interp_accel * A, double * RESULT) These functions return the numerical integral RESULT of an interpolated function over the range [A, B], using the interpolation object INTERP, data arrays XA and YA and the accelerator A.  File: gsl-ref.info, Node: Higher-level interface, Next: Interpolation Example programs, Prev: Evaluation of interpolating functions, Up: Interpolation Higher-level interface ====================== The functions described in the previous sections required the user to supply pointers to the x and y arrays on each call. The following functions are equivalent to the corresponding `gsl_interp' functions but maintain a copy of this data in the `gsl_spline' object. This removes the need to pass both XA and YA as arguments on each evaluation. These functions are defined in the header file `gsl_spline.h'. - Function: gsl_spline * gsl_spline_alloc (const gsl_interp_type * T, size_t N) - Function: int gsl_spline_init (gsl_spline * SPLINE, const double xa[], const double YA[], size_t SIZE) - Function: void gsl_spline_free (gsl_spline * SPLINE) - Function: double gsl_spline_eval (const gsl_spline * SPLINE, double X, gsl_interp_accel * A) - Function: int gsl_spline_eval_e (const gsl_spline * SPLINE, double X, gsl_interp_accel * A, double * Y) - Function: double gsl_spline_eval_deriv (const gsl_spline * SPLINE, double X, gsl_interp_accel * A) - Function: int gsl_spline_eval_deriv_e (const gsl_spline * SPLINE, double X, gsl_interp_accel * A, double * D) - Function: double gsl_spline_eval_deriv2 (const gsl_spline * SPLINE, double X, gsl_interp_accel * A) - Function: int gsl_spline_eval_deriv2_e (const gsl_spline * SPLINE, double X, gsl_interp_accel * A, double * D2) - Function: double gsl_spline_eval_integ (const gsl_spline * SPLINE, double A, double B, gsl_interp_accel * ACC) - Function: int gsl_spline_eval_integ_e (const gsl_spline * SPLINE, double A, double B, gsl_interp_accel * ACC, double * RESULT)  File: gsl-ref.info, Node: Interpolation Example programs, Next: Interpolation References and Further Reading, Prev: Higher-level interface, Up: Interpolation Examples ======== The following program demonstrates the use of the interpolation and spline functions. It computes a cubic spline interpolation of the 10-point dataset (x_i, y_i) where x_i = i + \sin(i)/2 and y_i = i + \cos(i^2) for i = 0 \dots 9. #include #include #include #include #include #include int main (void) { int i; double xi, yi, x[10], y[10]; printf ("#m=0,S=2\n"); for (i = 0; i < 10; i++) { x[i] = i + 0.5 * sin (i); y[i] = i + cos (i * i); printf ("%g %g\n", x[i], y[i]); } printf ("#m=1,S=0\n"); { gsl_interp_accel *acc = gsl_interp_accel_alloc (); gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, 10); gsl_spline_init (spline, x, y, 10); for (xi = x[0]; xi < x[9]; xi += 0.01) { yi = gsl_spline_eval (spline, xi, acc); printf ("%g %g\n", xi, yi); } gsl_spline_free (spline); gsl_interp_accel_free(acc); } return 0; } The output is designed to be used with the GNU plotutils `graph' program, $ ./a.out > interp.dat $ graph -T ps < interp.dat > interp.ps The result shows a smooth interpolation of the original points. The interpolation method can changed simply by varying the first argument of `gsl_spline_alloc'.  File: gsl-ref.info, Node: Interpolation References and Further Reading, Prev: Interpolation Example programs, Up: Interpolation References and Further Reading ============================== Descriptions of the interpolation algorithms and further references can be found in the following book, C.W. Ueberhuber, `Numerical Computation (Volume 1), Chapter 9 "Interpolation"', Springer (1997), ISBN 3-540-62058-3.  File: gsl-ref.info, Node: Numerical Differentiation, Next: Chebyshev Approximations, Prev: Interpolation, Up: Top Numerical Differentiation ************************* The functions described in this chapter compute numerical derivatives by finite differencing. An adaptive algorithm is used to find the best choice of finite difference and to estimate the error in the derivative. These functions are declared in the header file `gsl_diff.h' * Menu: * Numerical Differentiation functions:: * Numerical Differentiation Example:: * Numerical Differentiation References::  File: gsl-ref.info, Node: Numerical Differentiation functions, Next: Numerical Differentiation Example, Up: Numerical Differentiation Functions ========= - Function: int gsl_diff_central (const gsl_function *F, double X, double *RESULT, double *ABSERR) This function computes the numerical derivative of the function F at the point X using an adaptive central difference algorithm. The derivative is returned in RESULT and an estimate of its absolute error is returned in ABSERR. - Function: int gsl_diff_forward (const gsl_function *F, double X, double *RESULT, double *ABSERR) This function computes the numerical derivative of the function F at the point X using an adaptive forward difference algorithm. The function is evaluated only at points greater than X and at X itself. The derivative is returned in RESULT and an estimate of its absolute error is returned in ABSERR. This function should be used if f(x) has a singularity or is undefined for values less than X. - Function: int gsl_diff_backward (const gsl_function *F, double X, double *RESULT, double *ABSERR) This function computes the numerical derivative of the function F at the point X using an adaptive backward difference algorithm. The function is evaluated only at points less than X and at X itself. The derivative is returned in RESULT and an estimate of its absolute error is returned in ABSERR. This function should be used if f(x) has a singularity or is undefined for values greater than X.  File: gsl-ref.info, Node: Numerical Differentiation Example, Next: Numerical Differentiation References, Prev: Numerical Differentiation functions, Up: Numerical Differentiation Example ======= The following code estimates the derivative of the function f(x) = x^{3/2} at x=2 and at x=0. The function f(x) is undefined for x<0 so the derivative at x=0 is computed using `gsl_diff_forward'. #include #include #include double f (double x, void * params) { return pow (x, 1.5); } int main (void) { gsl_function F; double result, abserr; F.function = &f; F.params = 0; printf("f(x) = x^(3/2)\n"); gsl_diff_central (&F, 2.0, &result, &abserr); printf("x = 2.0\n"); printf("f'(x) = %.10f +/- %.5f\n", result, abserr); printf("exact = %.10f\n\n", 1.5 * sqrt(2.0)); gsl_diff_forward (&F, 0.0, &result, &abserr); printf("x = 0.0\n"); printf("f'(x) = %.10f +/- %.5f\n", result, abserr); printf("exact = %.10f\n", 0.0); return 0; } Here is the output of the program, $ ./demo f(x) = x^(3/2) x = 2.0 f'(x) = 2.1213203435 +/- 0.01490 exact = 2.1213203436 x = 0.0 f'(x) = 0.0012172897 +/- 0.05028 exact = 0.0000000000  File: gsl-ref.info, Node: Numerical Differentiation References, Prev: Numerical Differentiation Example, Up: Numerical Differentiation References and Further Reading ============================== The algorithms used by these functions are described in the following book, S.D. Conte and Carl de Boor, `Elementary Numerical Analysis: An Algorithmic Approach', McGraw-Hill, 1972.  File: gsl-ref.info, Node: Chebyshev Approximations, Next: Series Acceleration, Prev: Numerical Differentiation, Up: Top Chebyshev Approximations ************************ This chapter describes routines for computing Chebyshev approximations to univariate functions. A Chebyshev approximation is a truncation of the series f(x) = \sum c_n T_n(x), where the Chebyshev polynomials T_n(x) = \cos(n \arccos x) provide an orthogonal basis of polynomials on the interval [-1,1] with the weight function 1 / \sqrt{1-x^2}. The first few Chebyshev polynomials are, T_0(x) = 1, T_1(x) = x, T_2(x) = 2 x^2 - 1. The functions described in this chapter are declared in the header file `gsl_chebyshev.h'. * Menu: * The gsl_cheb_series struct:: * Creation and Calculation of Chebyshev Series:: * Chebyshev Series Evaluation:: * Derivatives and Integrals:: * Chebyshev Approximation examples:: * Chebyshev Approximation References and Further Reading::  File: gsl-ref.info, Node: The gsl_cheb_series struct, Next: Creation and Calculation of Chebyshev Series, Up: Chebyshev Approximations The gsl_cheb_series struct ========================== A Chebyshev series is stored using the following structure, typedef struct { double * c; /* coefficients c[0] .. c[order] */ int order; /* order of expansion */ double a; /* lower interval point */ double b; /* upper interval point */ } gsl_cheb_struct The approximation is made over the range [a,b] using ORDER+1 terms, including the coefficient c[0].  File: gsl-ref.info, Node: Creation and Calculation of Chebyshev Series, Next: Chebyshev Series Evaluation, Prev: The gsl_cheb_series struct, Up: Chebyshev Approximations Creation and Calculation of Chebyshev Series ============================================ - Function: gsl_cheb_series * gsl_cheb_alloc (const size_t N) This function allocates space for a Chebyshev series of order N and returns a pointer to a new `gsl_cheb_series' struct. - Function: void gsl_cheb_free (gsl_cheb_series * CS) This function frees a previously allocated Chebyshev series CS. - Function: int gsl_cheb_init (gsl_cheb_series * CS, const gsl_function * F, const double A, const double B) This function computes the Chebyshev approximation CS for the function F over the range (a,b) to the previously specified order. The computation of the Chebyshev approximation is an O(n^2) process, and requires n function evaluations.  File: gsl-ref.info, Node: Chebyshev Series Evaluation, Next: Derivatives and Integrals, Prev: Creation and Calculation of Chebyshev Series, Up: Chebyshev Approximations Chebyshev Series Evaluation =========================== - Function: double gsl_cheb_eval (const gsl_cheb_series * CS, double X) This function evaluates the Chebyshev series CS at a given point X. - Function: int gsl_cheb_eval_err (const gsl_cheb_series * CS, const double X, double * RESULT, double * ABSERR) This function computes the Chebyshev series CS at a given point X, estimating both the series RESULT and its absolute error ABSERR. The error estimate is made from the first neglected term in the series. - Function: double gsl_cheb_eval_n (const gsl_cheb_series * CS, size_t ORDER, double X) This function evaluates the Chebyshev series CS at a given point N, to (at most) the given order ORDER. - Function: int gsl_cheb_eval_n_err (const gsl_cheb_series * CS, const size_t ORDER, const double X, double * RESULT, double * ABSERR) This function evaluates a Chebyshev series CS at a given point X, estimating both the series RESULT and its absolute error ABSERR, to (at most) the given order ORDER. The error estimate is made from the first neglected term in the series.  File: gsl-ref.info, Node: Derivatives and Integrals, Next: Chebyshev Approximation examples, Prev: Chebyshev Series Evaluation, Up: Chebyshev Approximations Derivatives and Integrals ========================= The following functions allow a Chebyshev series to be differentiated or integrated, producing a new Chebyshev series. Note that the error estimate produced by evaluating the derivative series will be underestimated due to the contribution of higher order terms being neglected. - Function: int gsl_cheb_calc_deriv (gsl_cheb_series * DERIV, const gsl_cheb_series * CS) This function computes the derivative of the series CS, storing the derivative coefficients in the previously allocated DERIV. The two series CS and DERIV must have been allocated with the same order. - Function: int gsl_cheb_calc_integ (gsl_cheb_series * INTEG, const gsl_cheb_series * CS) This function computes the integral of the series CS, storing the integral coefficients in the previously allocated INTEG. The two series CS and INTEG must have been allocated with the same order. The lower limit of the integration is taken to be the left hand end of the range A.  File: gsl-ref.info, Node: Chebyshev Approximation examples, Next: Chebyshev Approximation References and Further Reading, Prev: Derivatives and Integrals, Up: Chebyshev Approximations Examples ======== The following example program computes Chebyshev approximations to a step function. This is an extremely difficult approximation to make, due to the discontinuity, and was chosen as an example where approximation error is visible. For smooth functions the Chebyshev approximation converges extremely rapidly and errors would not be visible. #include #include #include double f (double x, void *p) { if (x < 0.5) return 0.25; else return 0.75; } int main (void) { int i, n = 10000; gsl_cheb_series *cs = gsl_cheb_alloc (40); gsl_function F; F.function = f; F.params = 0; gsl_cheb_init (cs, &F, 0.0, 1.0); for (i = 0; i < n; i++) { double x = i / (double)n; double r10 = gsl_cheb_eval_n (cs, 10, x); double r40 = gsl_cheb_eval (cs, x); printf ("%g %g %g %g\n", x, GSL_FN_EVAL (&F, x), r10, r40); } gsl_cheb_free (cs); return 0; } The output from the program gives the original function, 10-th order approximation and 40-th order approximation, all sampled at intervals of 0.001 in x.  File: gsl-ref.info, Node: Chebyshev Approximation References and Further Reading, Prev: Chebyshev Approximation examples, Up: Chebyshev Approximations References and Further Reading ============================== The following paper describes the use of Chebyshev series, R. Broucke, "Ten Subroutines for the Manipulation of Chebyshev Series [C1] (Algorithm 446)". `Communications of the ACM' 16(4), 254-256 (1973)  File: gsl-ref.info, Node: Series Acceleration, Next: Discrete Hankel Transforms, Prev: Chebyshev Approximations, Up: Top Series Acceleration ******************* The functions described in this chapter accelerate the convergence of a series using the Levin u-transform. This method takes a small number of terms from the start of a series and uses a systematic approximation to compute an extrapolated value and an estimate of its error. The u-transform works for both convergent and divergent series, including asymptotic series. These functions are declared in the header file `gsl_sum.h'. * Menu: * Acceleration functions:: * Acceleration functions without error estimation:: * Example of accelerating a series:: * Series Acceleration References::  File: gsl-ref.info, Node: Acceleration functions, Next: Acceleration functions without error estimation, Up: Series Acceleration Acceleration functions ====================== The following functions compute the full Levin u-transform of a series with its error estimate. The error estimate is computed by propagating rounding errors from each term through to the final extrapolation. These functions are intended for summing analytic series where each term is known to high accuracy, and the rounding errors are assumed to originate from finite precision. They are taken to be relative errors of order `GSL_DBL_EPSILON' for each term. The calculation of the error in the extrapolated value is an O(N^2) process, which is expensive in time and memory. A faster but less reliable method which estimates the error from the convergence of the extrapolated value is described in the next section For the method described here a full table of intermediate values and derivatives through to O(N) must be computed and stored, but this does give a reliable error estimate. . - Function: gsl_sum_levin_u_workspace * gsl_sum_levin_u_alloc (size_t N) This function allocates a workspace for a Levin u-transform of N terms. The size of the workspace is O(2n^2 + 3n). - Function: int gsl_sum_levin_u_free (gsl_sum_levin_u_workspace * W) This function frees the memory associated with the workspace W. - Function: int gsl_sum_levin_u_accel (const double * ARRAY, size_t ARRAY_SIZE, gsl_sum_levin_u_workspace * W, double * SUM_ACCEL, double * ABSERR) This function takes the terms of a series in ARRAY of size ARRAY_SIZE and computes the extrapolated limit of the series using a Levin u-transform. Additional working space must be provided in W. The extrapolated sum is stored in SUM_ACCEL, with an estimate of the absolute error stored in ABSERR. The actual term-by-term sum is returned in `w->sum_plain'. The algorithm calculates the truncation error (the difference between two successive extrapolations) and round-off error (propagated from the individual terms) to choose an optimal number of terms for the extrapolation.