Programowanie w systemie UNIX/GSL
Instalacja
edytuj- ze źródeł
- skompilowane pakiety
pakiety
edytujW 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
edytujPoniż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
edytujręczna
edytujgcc -L/usr/lib -lgsl -lgslcblas -lm example.c
z użyciem gsl-config
edytujgcc $(gsl-config --cflags) example.c $(gsl-config --libs)
make
edytujOpis [5]
Wynik pracy programu
edytujWynik powinien być poprawny dla podwójnej precyzji:
J0(5) = -1.775967713143382920e-01
Silnia
edytujW 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ń
edytujWielowymiarowe 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
edytujInterpolacja = 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- dokumentacja
- 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
/*
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
edytujZestaw 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
edytujFuncje 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- ↑ Ask Ubuntu : Install GNU Scientific library (GSL) on Ubuntu 14.04 via terminal
- ↑ askubuntu question :how-can-i-run-an-c-c-program-that-use-gnu-scientific-library-gsl
- ↑ Funkcje Bessela w wikipedii
- ↑ An Example Program – GNU Scientific Library – Reference Manual
- ↑ Getting started with GSL by Darren Wilkinson
- ↑ gsl: Example-programs-for-Multidimensional-Root-finding