Instalacja

ze źródeł
skompilowane pakiety

pakiety

W Ubuntu 14.04 są pakiety:^[1]

gsl-bin: GNU Scientific Library (GSL) -- binary package
libgsl0-dbg: GNU Scientific Library (GSL) -- debug symbols package
libgsl0-dev: GNU Scientific Library (GSL) -- development package
libgsl0ldbl: GNU Scientific Library (GSL) -- library package

Instalacja

 sudo apt-get install libgsl0ldbl
 sudo apt-get install libgsl0-dev

Sprawdzamy biblioteki:

gsl-config --libs

Otrzymujemy informacje o katalogu i opcjach linkera:

-L/usr/lib -lgsl -lgslcblas -lm

Inna metoda: ^[2]

sudo apt-get install gsl-bin
sudo apt-get -y install libgsl-dev

Pierwszy program

Poniższy przykładowy program oblicza wartość funkcji Bessela^[3] dla argumentu 5^[4]:

#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>

int main(void)
{
  double x = 5.0;
  double y = gsl_sf_bessel_J0(x);
  printf("J0(%g) = %.18e\n", x, y);
  return 0;
}

Kompilacja

Compile GSL C Code

ręczna

gcc -L/usr/lib -lgsl -lgslcblas -lm example.c

z użyciem gsl-config

gcc $(gsl-config --cflags) example.c $(gsl-config --libs)

make

Opis ^[5]

Wynik pracy programu

Wynik powinien być poprawny dla podwójnej precyzji:

J0(5) = -1.775967713143382920e-01

Silnia

W formie naturalnej

#include <stdio.h>
#include <gsl/gsl_sf.h>


// gcc -I/usr/include/gsl/ -L/usr/local/lib/ -lgsl -lgslcblas -lm f.c
// http://www.gnu.org/software/gsl/manual/html_node/Special-Function-Usage.html#Special-Function-Usage
// The natural form returns only the value of the function and can be used directly in mathematical expressions.

int main(void)
{
  unsigned int n = 100;
  double result;

  result = gsl_sf_fact(n);


 
  printf("factorial( %d ) = %.0f\n", n,result);
  

  return 0;
}

Wynik działania:

./a.out

factorial( 100 ) = 933262154439441509656467047959538825784009703731840988310128895405822272385704312950661130
89288327277825849664006524270554535976289719382852181865895959724032

W formie wychwytującej błędy

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf.h>


// gcc -I/usr/include/gsl/ -L/usr/local/lib/ -lgsl -lgslcblas -lm f.c
// gcc -I/usr/include/gsl/ -L/usr/lib/ -lgsl -lgslcblas -lm f.c
// factorial
// http://www.gnu.org/software/gsl/manual/html_node/Special-Function-Usage.html#Special-Function-Usage
// The error-handling function
// http://www.gnu.org/software/gsl/manual/html_node/Special-Functions-Examples.html#Special-Functions-Examples
/*




*/
 
int main(void)
{
  unsigned int n = 100;
  int status; 
  gsl_sf_result result;

  status = gsl_sf_fact_e (n, &result);


  printf("status  = %s\n", gsl_strerror(status));
  printf("factorial( %2d ) = %f\n", n,result.val);
  printf("+/- % .18f\n",result.err); 

  return 0;
}

Wynik:

./a.out
status  = success
factorial( 100 ) =    933262154439441509656467047959538825784009703731840988310128895405822272385704312950661130
89288327277825849664006524270554535976289719382852181865895959724032
+/-  4144516527479786869998377376409676501474209436767641660238087976812259116180322817471642386
      5792481641279576393802186138583881092629521428905984.000000000000000000

Znajdowanie zer równań

Wielowymiarowe metody^[6]

Aproksymacja wielomianowa

/*
https://rosettacode.org/wiki/Polynomial_regression
Polynomial regression

Find an approximating polynomial of known degree for a given data.

Example: For input data:

x = {0,  1,  2,  3,  4,  5,  6,   7,   8,   9,   10};
y = {1,  6,  17, 34, 57, 86, 121, 162, 209, 262, 321};

The approximating polynomial is: 3x^2 + 2x + 1
Here, the polynomial's coefficients are (3, 2, 1).

one can check it online: https://arachnoid.com/polysolve/

=========================================

gcc -Wall -I/usr/local/include  -c p.c
gcc -L/usr/local/lib p.o -lgsl -lgslcblas -lm

./a.out
1.000000
2.000000
3.000000

*/

#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit.h> // 
#include <stdbool.h>
#include <math.h>


// http://www.gnu.org/software/gsl/manual/html_node/Fitting-Examples.html here] helped alot to code this quickly):



#define NP 11
double x[] = {0,  1,  2,  3,  4,  5,  6,   7,   8,   9,   10};
double y[] = {1,  6,  17, 34, 57, 86, 121, 162, 209, 262, 321};

#define DEGREE 3
double coeff[DEGREE];


bool polynomialfit(int obs, int degree, 
		   double *dx, double *dy, double *store) /* n, p */
{
  gsl_multifit_linear_workspace *ws;
  gsl_matrix *cov, *X;
  gsl_vector *y, *c;
  double chisq;

  int i, j;

  X = gsl_matrix_alloc(obs, degree);
  y = gsl_vector_alloc(obs);
  c = gsl_vector_alloc(degree);
  cov = gsl_matrix_alloc(degree, degree);

  for(i=0; i < obs; i++) {
    for(j=0; j < degree; j++) {
      gsl_matrix_set(X, i, j, pow(dx[i], j));
    }
    gsl_vector_set(y, i, dy[i]);
  }

  ws = gsl_multifit_linear_alloc(obs, degree);
  gsl_multifit_linear(X, y, c, cov, &chisq, ws);

  /* store result ... */
  for(i=0; i < degree; i++)
  {
    store[i] = gsl_vector_get(c, i);
  }

  gsl_multifit_linear_free(ws);
  gsl_matrix_free(X);
  gsl_matrix_free(cov);
  gsl_vector_free(y);
  gsl_vector_free(c);
  return true; /* we do not "analyse" the result (cov matrix mainly)
		  to know if the fit is "good" */
}

int main()
{
  int i;

  polynomialfit(NP, DEGREE, x, y, coeff);
  for(i=0; i < DEGREE; i++) {
    printf("%lf\n", coeff[i]);
  }
  return 0;
}

interpolacja

Interpolacja = znajdowanie wzoru funkcji przechodzącej przez n punktów

Interpolacja typy

wg wymiaru (1D / 2D)
wg metody

Metody interpolacji 1D:

gsl_interp_linear = interpolacja liniowa
gsl_interp_polynomial = Interpolacja wielomianowa ( ang. polynomial interpolation) = interpolacja za pomocą wielmianów
gsl_interp_cspline = Cubic splines = Interpolacja funkcjami sklejanymi ( 3 stopnia = cubic )
gsl_interp_cspline_periodic
gsl_interp_akima = Akima splines
gsl_interp_akima_periodic
gsl_interp_steffen = Steffen splines

Poniższy program demonstruje użycie funkcji interpolacji funkcjami sklejanymi. Oblicza interpolację splajnu sześciennego 10-punktowego zbioru danych $(x_{i},y_{i})$ gdzie

$x_{i}=i+\sin(i)/2$
$y_{i}=i+\cos(i^{2})$
$i=0\dots 9$

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

int
main (void)
{
  int i;
  double xi, yi, x[10], y[10];

  printf ("#m=0,S=17\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;
}

Kompilacja

gcc -Wall s.c -lm -lgsl -lgslcblas

Uruchomienie

./a.out > interp.dat

Utworzenie grafiki SVG za pomocą programu graph

graph -T svg < interp.dat > interp.svg

ODE

    
/*

gcc o.c -Wall -Wextra -lm -lplplot -lgsl -lgslcblas


*/

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_statistics.h>
#include <plplot/plplot.h>






int ode_func(double t, const double y[], double f[], void *params) {
    double k = *((double*)params);
    (void)t;
    f[0] = y[1];
    f[1] = k * y[1] * (1.0 - y[0] * y[0]) - y[0];
    return GSL_SUCCESS;
}

/* Some methods ask for you to give the Jacobian of it too
 * for better performance. You have to calculate this by hand.
 */
int ode_jac(double t, const double y[], double *dfdy, double dfdt[], void *params) {
    double k = *((double*)params);
    (void)t;
    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, -(1.0 + 2.0 * k * y[0] * y[1]));
    gsl_matrix_set(m, 1, 1, k * (1.0 * y[0] * y[0]));
    /* Autonomous. */
    dfdt[0] = 0.0;
    dfdt[1] = 0.0;
    return GSL_SUCCESS;
}

int main(int argc, char **argv) {
    enum Constexpr {n_points = 1000};
    double mu = 1.0;
    gsl_odeiv2_system sys = {ode_func, ode_jac, 2, &mu};
    gsl_odeiv2_driver * d= gsl_odeiv2_driver_alloc_y_new(
        &sys,
        gsl_odeiv2_step_rk8pd,
        1e-6,
        1e-6,
        0.0
    );
    int i;
    double t = 0.0;
    double t1 = 100.0;
    /* Initial condition: f = 0 with speed 0. */
    double y[2] = {1.0, 0.0};
    double dt = t1 / n_points;
    double datax[n_points];
    double datay[n_points];

    for (i = 0; i < n_points; i++) {
        double ti = i * dt;
        int status = gsl_odeiv2_driver_apply(d, &t, ti, y);
        if (status != GSL_SUCCESS) {
            fprintf(stderr, "error, return value=%d\n", status);
            break;
        }

        /* Get output. */
        printf("%.5e %.5e %.5e\n", t, y[0], y[1]);
        datax[i] = y[0];
        datay[i] = y[1];
    }

    /* Cleanup. */
    gsl_odeiv2_driver_free(d);

    /* Plot. */
    if (argc > 1 && argv[1][0] == '1') {
        plsdev("xwin");
        plinit();
        plenv(
            gsl_stats_min(datax, 1, n_points),
            gsl_stats_max(datax, 1, n_points),
            gsl_stats_min(datay, 1, n_points),
            gsl_stats_max(datay, 1, n_points),
            0,
            1
        );
        plstring(n_points, datax, datay, "*");
        plend();
    }

    return EXIT_SUCCESS;
}

Funkcje specjalne

Zestaw wszystkich funkcji specjalnych ( ang special functions) biblioteki GSL jest dostępny poprzez plik gsl_sf.h

Lista funkcji:

Airy functions, plik gsl_sf_airy.h
Bessel functions, plik gsl_sf_bessel.h
Clausen functions
Coulomb wave functions
Coupling coefficients
Dawson function
Debye functions
Dilogarithms
Elliptic integrals
Jacobi elliptic functions
Error functions
Exponential integrals
Fermi-Dirac functions
Gamma functions
Gegenbauer functions
Hermite polynomials and functions
Hypergeometric functions
Laguerre functions
Legendre functions and Spherical Harmonics
Psi (Digamma) Function
Synchrotron functions
Transport functions
Trigonometric functions
Zeta functions

Formy użycia

natural form: funkcja zwraca tylko wartość funkcji
an error-handling form: funkcja zwraca wartość i kod błędu

Przykłady:

double y = gsl_sf_bessel_J0 (x); // natural form

// an error-handling form
gsl_sf_result result;
int status = gsl_sf_bessel_J0_e (x, &result);

Bessel functions

Funcje Bessela

the Cylindrical Bessel functions $J_{n}(x)$ , $Y_{n}(x)$
Modified cylindrical Bessel functions $I_{n}(x)$ , $K_{n}(x)$
Spherical Bessel functions $j_{l}(x)$ , $y_{l}(x)$
Modified Spherical Bessel functions $i_{l}(x)$ , $k_{l}(x)$

Funkcja Bessela pierwszego rodzaju rzędu 0 (ang. the regular cylindrical Bessel function of zeroth order) $J_{0}(x)$

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_sf_bessel.h>

int main(void) {
    double x = 5.0;
    printf("J0(%g) = %.18e\n", x, gsl_sf_bessel_J0(x));
    return EXIT_SUCCESS;
}

Kompilacja:

gcc -Wall b.c -lm -lgsl -lgslcblas

Uruchomienie:

./a.out

Wynik:

J0(5) = -1.775967713143382642e-01

Odnośniki

[1] Ask Ubuntu : Install GNU Scientific library (GSL) on Ubuntu 14.04 via terminal

[2] skubuntu question :how-can-i-run-an-c-c-program-that-use-gnu-scientific-library-gsl

[3] Funkcje Bessela w wikipedii

[4] An Example Program – GNU Scientific Library – Reference Manual

[5] Getting started with GSL by Darren Wilkinson

[6] sl: Example-programs-for-Multidimensional-Root-finding

[1]

[2]

[3]

[4]

[5]

[6]