Tekst źródłowy strony Programowanie w systemie UNIX/GSL

=Instalacja=
* ze źródeł
* skompilowane pakiety

==pakiety==
W Ubuntu 14.04 są pakiety:<ref>[http://askubuntu.com/questions/490465/install-gnu-scientific-library-gsl-on-ubuntu-14-04-via-terminal Ask Ubuntu : Install GNU Scientific library (GSL) on Ubuntu 14.04 via terminal]</ref>
* 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: <ref>[http://askubuntu.com/questions/766076/how-can-i-run-an-c-c-program-that-use-gnu-scientific-library-gsl askubuntu question :how-can-i-run-an-c-c-program-that-use-gnu-scientific-library-gsl]</ref>

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

=Pierwszy program = 

Poniższy przykładowy program oblicza wartość funkcji Bessela<ref>[[w:Funkcje Bessela|Funkcje Bessela w wikipedii]]</ref> dla argumentu 5<ref>[http://www.gnu.org/software/gsl/manual/html_node/An-Example-Program.html An Example Program – GNU Scientific Library – Reference Manual<!-- Tytuł wygenerowany przez bota -->]</ref>:

<syntaxhighlight lang="c">
#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;
}
</syntaxhighlight>

==Kompilacja == 
* [https://doc.cocalc.com/howto/gsl.html 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 <ref>[http://www.mas.ncl.ac.uk/~ndjw1/teaching/sim/gsl.html Getting started with GSL by Darren Wilkinson ]</ref>

==Wynik pracy programu== 
Wynik powinien być  poprawny dla podwójnej precyzji:


 J0(5) = -1.775967713143382920e-01

=Silnia=

==W formie naturalnej ==

<syntaxhighlight lang=c>
#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;
}
</syntaxhighlight>


Wynik działania:

 ./a.out

 factorial( 100 ) = 933262154439441509656467047959538825784009703731840988310128895405822272385704312950661130
 89288327277825849664006524270554535976289719382852181865895959724032


==W formie wychwytującej błędy==

<syntaxhighlight lang=c>
#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;
}
</syntaxhighlight>


Wynik:

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

=Znajdowanie zer równań=
Wielowymiarowe metody<ref>[http://www.gnu.org/software/gsl/manual/html_node/Example-programs-for-Multidimensional-Root-finding.html gsl: Example-programs-for-Multidimensional-Root-finding]</ref>

= Aproksymacja wielomianowa=

<syntaxhighlight lang=c>
/*
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;
}
</syntaxhighlight>
=interpolacja=
[[Plik:Cubic spline 1D interpolation.svg|mały|prawo|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 <math>(x_i, y_i)</math> gdzie 
* <math>x_i = i + \sin(i)/2 </math> 
* <math>  y_i = i + \cos(i^2)</math> 
* <math> i = 0 \dots 9</math>

<syntaxhighlight lang="c">
#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;
}
</syntaxhighlight>

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=
* [https://www.gnu.org/software/gsl/doc/html/ode-initval.html dokumentacja] 
* [https://people.sc.fsu.edu/~jburkardt/c%20src/midpoint%20gsl%20test/midpoint%20gsl%20test.html midpoint_gsl_test, a C code which calls the Gnu Scientific Library (GSL) implicit midpoint method solver for ordinary differential equation (ODE), and uses gnuplot() to plot the resulting solution by   j burkardt]

[[Plik:Van der Pol oscillator solution and phase space generated with GNU Scientific Library.png|mały|prawo|Van der Pol oscylator]]


<syntaxhighlight lang= c>    
/*

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;
}

</syntaxhighlight >

=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 <math>J_n(x)</math> , <math>Y_n(x)</math>
* Modified cylindrical Bessel functions <math>I_n(x)</math> ,  <math>K_n(x)</math>
* Spherical Bessel functions <math>j_l(x)</math>, <math>y_l(x)</math>
* Modified Spherical Bessel functions <math>i_l(x)</math>, <math>k_l(x)</math>


Funkcja Bessela pierwszego rodzaju rzędu 0 (ang. the regular cylindrical Bessel function of zeroth order) <math>J_0(x)</math>

<syntaxhighlight lang="c">
#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;
}
</syntaxhighlight>


Kompilacja:

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

Uruchomienie:

 ./a.out

Wynik:

 J0(5) = -1.775967713143382642e-01

=Odnośniki=
<references/>