[Help-gsl] Gsl-1.8 compiled and running under MSYS-1.10+second order derivative

[Rene Girard]

Hi,

  First I do not know if I used the proper e-mail
address for the following e-mail. If it is not
please let me know which one I should use. Also, I
apologize if you find this
e-mail long. 

First I would like to inform you that I was successful
in compiling, linking and using 
gsl-1.8 under MSYS 1.10 (available from Mingw.org)
running under Windows XP, 
Windows XP Pro and Windows 2000. The Web site for
gsl-1.8 indicates only Cygwin
for running gsl-1.8 that under windows.      

In using the functions in gsl-1.8 for numerical
differentiation in Chapter 27, under Fedora Core 5 
using gcc 4.1.0 on AMD64 Athlon 3700, I modified the
example in the manual where I define and 
use the pointer *params. The modified example now
reads:

- Begin program t_dn1.c - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - - 

/*
 * Programme t_dn1.c: Le programme t_dn1.c donne une
methode d'utilisation 
 *                    des fonctions de la librarie GSL
pour le calcul de la 
 *                    derivee d'une fonction a une
variable y = f(x). Dans le
 *                    present on a : 
 * 
 *                          f(x) = x^(3/2)
 *
 *                    avec f'(x) = 3*sqrt(x)/2
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>


double g1(double, double *);

int main(void)
 {
  gsl_function f;
  double a,r1,er1;
  double p1[2];

  p1[0] = 1.0; p1[1] = 1.0e-04;

  f.function = (double (*) ()) &g1;
  f.params = &p1[0];

  printf("f(x) = x^(3/2)\n");


  a = p1[0];

  gsl_deriv_central(&f,2.0,p1[1],&r1,&er1);
  printf("Pour x = 2.0\n");
  printf("f'(x)         = %14.10f +/-
%14.10f\n",r1,er1);
  printf("Valeur exacte = %14.10f\n",a*1.5*sqrt(2.0));

  gsl_deriv_forward(&f,0.0,1.0e-10,&r1,&er1);
  printf("Pour x = 0.0\n");
  printf("f'(x)         = %14.10f +/-
%14.10f\n",r1,er1);
  printf("Valeur exacte = %14.10f\n",0.0);

  return 0;

 } /* fin du programme t_dn1.c */       

double g1(double x, double *params)
 {
  double a,y;

  a = params[0];
  y = a*pow(x,1.5);
  return y;  
 }
 
- End program t_dn1c - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - - - 

I would like to bring to your attention that when I
compile that simple example and
I set the pointer to the function from we want to have
the derivative with the following
statement:
 
  f.function = &g1; (like in the example in Chap. 27)

I get the following warning:

"warning: assignment from incompatible pointer type"

to remove the warning I had to cast the pointer as
follow:

 f.function = (double (*) ()) &g1;

Note that I found the pointer type using gdb. I do not
think this is a bug in 
gsl-1.8 and I would like to suggest that the gsl-1.8
users should be informed 
of this fact.

I was also able to use the numerical differentiation
function "gsl_deriv_central" to compute
first, second order partial derivatives of a function
with two variables. Also, I was able to compute 
the second order cross derivative. The programs are
listed below and they can be used as 
useful examples for the other users of gsl-1.8 on how
to use the numerical differentiation functions in 
gsl-1.8 to compute these derivatives. The test
function used in these examples is:

             f(x,y) = (sin(pi*x/2)*sin(pi*y/2))^2


Program t_dn2.c compute first order partial derivative
with respect to x of
f(x,y):

- Begin program t_dn2.c - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - - 

/*
 * Programme t_dn2.c: Le programme t_dn2.c donne une
methode d'utilisation 
 *                    des fonctions de la librarie GSL
pour le calcul de la 
 *                    derivee partielle du premier
ordre d'une fonction a deux variables 
 *                    z = f(x,y) par rapport a x. Dans
ce programme on considere 
 *                    la fonction test:
 * 
 *                          f(x,y) =
(sin(pi*x/2)*sin(pi*y/2))^2
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>


double g_x(double, double *);

int main(void)
 {
  gsl_function f_x;
  const double pi = 4.0*atan(1.0);     
  double ro,r1,dx,er1,x,y;
  double p1[2];
  int i;

  y = 0.65; 

  dx = 0.1;

  p1[0] = 1.0e-04;
  p1[1] = y;

  f_x.function = (double (*) ()) &g_x;
  f_x.params = &p1[0];

  x = -dx;
  printf("f(x,y) = (sin(pi*x/2)*sin(pi*y/2))^2\n");
  for(i = 1; i <= 11; i++)
   {
    x += dx;
    gsl_deriv_central(&f_x,x,p1[0],&r1,&er1);
    ro = pi*(sin(pi*p1[1]/2.0)*sin(pi*p1[1]/2.0))*
         cos(pi*x/2.0)*sin(pi*x/2.0);     
    printf("Pour x = %11.6f et y = %5.2f\n",x,y);
    printf("f_x           = %14.10f +/-
%14.10f\n",r1,er1);
    printf("Valeur exacte = %14.10f\n",ro); 
   }

  return 0;

 } /* fin du programme t_dn2.c */       

double g_x(double x, double *params)
 {
  const double pi = 4.0*atan(1.0);    
  double yo,y1;

  yo = (sin(pi*x/2.0)*sin(pi*params[1]/2));
  y1 = yo*yo;

  return y1;  
 }

- End program t_dn2.c - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - 

Program t_dn3.c compute second order partial
derivative with respect to x of
f(x,y):

- Begin program t_dn3.c - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - - 

/*
 * Programme t_dn3.c: Le programme t_dn3.c donne une
methode d'utilisation 
 *                    des fonctions de la librarie GSL
pour le calcul de la 
 *                    derivee du second ordre d'une
fonction a deux variables 
 *                    z = f(x,y). Dans ce programme on
considere la fonction test:
 * 
 *                          f(x,y) =
(sin(pi*x/2)*sin(pi*y/2))^2
 *
 *                    avec f_xx =
((pi^2)/2)*(sin(pi*y/2))^2*
 *                                ((cos(pi*x/2))^2 -
(sin(pi*x/2))^2)
 *
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>


double g_x(double, double *);
double g_xx(double, double *);

int main(void)
 {
  gsl_function f_xx;
  const double pi = 4.0*atan(1.0);     
  double ro,r1,ero,er1,pero,per1,x,y;
  double p1[2];

  x = 0.3; y = 0.65; 

  f_xx.function = (double (*) ()) &g_xx;
  f_xx.params = &p1[0];

  p1[0] = 1.0e-04;
  p1[1] = y;

  printf("f(x,y) = (sin(pi*x/2)*sin(pi*y/2))^2\n");
  gsl_deriv_central(&f_xx,x,p1[0],&r1,&er1);
  ro = (pi*pi/2.0)*(sin(pi*y/2.0)*sin(pi*y/2.0))*
       (cos(pi*x/2.0)*cos(pi*x/2.0) -
sin(pi*x/2.0)*sin(pi*x/2.0));     
  printf("Pour x = %5.2f et y = %5.2f\n",x,y);
  if (fabs(r1) >= 1.0e-06) 
   {
    per1 = er1*100.0/fabs(r1);
    printf("f_xx          = %14.10f +/-  %14.10f  per1
= %14.6f%%\n",r1,er1,per1);
   }
  else
   {
    printf("f_xx          = %14.10f +/- 
%14.10f\n",r1,er1);
   }

  ero = fabs(ro-r1);
  if (fabs(ro) >= 1.0e-06) 
   {
    pero = ero*100.0/fabs(ro);
    printf("Valeur exacte = %14.10f ero= %14.10f  pero
= %14.6f%%\n",ro,ero,pero);
   }
  else
   {
    printf("Valeur exacte = %14.10f ero=
%14.10f\n",ro,ero);
   }

  return 0;

 } /* fin du programme t_dn3.c */       

double g_xx(double x, double *params)
 {
  gsl_function f_x;
  double ero,y,p_x[2];
 
  p_x[0] = params[0];
  p_x[1] = params[1];

  f_x.function = (double (*) ()) &g_x;
  f_x.params = &p_x[0]; 

  gsl_deriv_central(&f_x,x,p_x[0],&y,&ero);

  return y;  
 }

double g_x(double x, double *params)
 {
  const double pi = 4.0*atan(1.0);    
  double yo,y1;

  yo = (sin(pi*x/2.0)*sin(pi*params[1]/2));
  y1 = yo*yo;

  return y1;  
 }

- End program t_dn3.c - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - 

Program t_dn4.c compute second order partial cross
derivative with respect to x and y of
f(x,y):

- Begin program t_dn4.c - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - - 

/*
 * Programme t_dn4.c: Le programme t_dn4.c donne une
methode d'utilisation 
 *                    des fonctions de la librarie GSL
pour le calcul de la 
 *                    derivee croisee du second ordre
d'une fonction a deux variables 
 *                    z = f(x,y). Dans ce programme on
considere la fonction test:
 * 
 *                          f(x,y) =
(sin(pi*x/2)*sin(pi*y/2))^2
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <adm/adm.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>


double g_y(double, double *);
double g_xy(double, double *);

int main(void)
 {
  gsl_function f_xy;
  const double pi = 4.0*atan(1.0);     
  double ro,r1,ero,er1,pero,per1,x,y;
  double p1[2];

  x = 0.3; y = 0.65; 

  f_xy.function = (double (*) ()) &g_xy;
  f_xy.params = &p1[0];

  p1[0] = 1.0e-04;
  p1[1] = y;

  printf("f(x,y) = (sin(pi*x/2)*sin(pi*y/2))^2\n");
  gsl_deriv_central(&f_xy,x,p1[0],&r1,&er1);
  ro =
(pi*pi)*(sin(pi*x/2.0)*cos(pi*x/2.0))*(sin(pi*y/2.0)*cos(pi*y/2.0));
  printf("Pour x = %5.2f et y = %5.2f\n",x,y);
  if (fabs(r1) >= 1.0e-06) 
   {
    per1 = er1*100.0/fabs(r1);
    printf("f_xy          = %14.10f +/-  %14.10f  per1
= %14.6f%%\n",r1,er1,per1);
   }
  else
   {
    printf("f_xy          = %14.10f +/- 
%14.10f\n",r1,er1);
   }

  ero = fabs(ro-r1);
  if (fabs(ro) >= 1.0e-06) 
   {
    pero = ero*100.0/fabs(ro);
    printf("Valeur exacte = %14.10f ero= %14.10f  pero
= %14.6f%%\n",ro,ero,pero);
   }
  else
   {
    printf("Valeur exacte = %14.10f ero=
%14.10f\n",ro,ero);
   }

  return 0;

 } /* fin du programme t_dn4.c */       

double g_xy(double x, double *params)
 {
  gsl_function f_y;
  double ero,y,p_y[2];
 
  p_y[0] = params[0];
  p_y[1] = x;

  f_y.function = (double (*) ()) &g_y;
  f_y.params = &p_y[0]; 

  gsl_deriv_central(&f_y,params[1],p_y[0],&y,&ero);

  return y;  
 }

double g_y(double y, double *params)
 {
  const double pi = 4.0*atan(1.0);    
  double yo,y1;

  yo = sin(pi*params[1]/2.0)*sin(pi*y/2);
  y1 = yo*yo;

  return y1;  
 }

- End program t_dn4.c - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - - - - - - - - 

The above programs are not optimized as they are
merely illustrations on a way to use
the functions in gsl-1.8 for numerical
differentiation. I hope you will find the above
useful and worth distributing to other users of
gsl-1.8.


Regards

Rene Girard 

__________________________________________________
Do You Yahoo!?
Tired of spam?  Yahoo! Mail has the best spam protection around 
http://mail.yahoo.com