Programowanie w systemie UNIX/GSL

Z Wikibooks, biblioteki wolnych podręczników.

Instalacja[edytuj]

  • ze źródeł
  • skompilowane pakiety

pakiety[edytuj]

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[edytuj]

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[edytuj]

ręczna[edytuj]

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

z użyciem gsl-config[edytuj]

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

make[edytuj]

Opis [5]

Wynik pracy programu[edytuj]

Wynik powinien być poprawny dla podwójnej precyzji:


J0(5) = -1.775967713143382920e-01

Silnia[edytuj]

W formie naturalnej[edytuj]

#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[edytuj]

#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ń[edytuj]

Wielowymiarowe metody[6]

Aproksymacja wielomianowa[edytuj]

/*
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[edytuj]

Interpolaja funckjami sklejanymi 1 D. Wynik działana programu

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 gdzie

#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[edytuj]

Van der Pol oscylator


    
/*

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[edytuj]

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[edytuj]

Funcje Bessela

  • the Cylindrical Bessel functions ,
  • Modified cylindrical Bessel functions ,
  • Spherical Bessel functions ,
  • Modified Spherical Bessel functions ,


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

#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[edytuj]

  1. Ask Ubuntu : Install GNU Scientific library (GSL) on Ubuntu 14.04 via terminal
  2. askubuntu 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. gsl: Example-programs-for-Multidimensional-Root-finding