[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Help-gsl] integrating a function with externally defined coefficients

From: Nathan Moore
Subject: [Help-gsl] integrating a function with externally defined coefficients
Date: Sat, 12 Feb 2005 12:16:32 -0600


I'd like to use GSL to numerically integrate a function with a bunch of coefficients that I defined as "global" variables. At present the integration doesn't work, which I imagine might have something to do with the "gsl_function" type - which I don't really understand.

Specifically, I define a set of coefficients as globals,


// and then inside the main function I assign them values,

main() {

    extern double Ax[FREQUENCY_LIMIT];

 for (i = 0; i < FREQUENCY_LIMIT; i++) {
           Ax[i] = gsl_rng_uniform(dice);

// later, I integrate a function which uses these coefficients,

gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
    double *length;
    double *abserror;

    //perform integration with GSL adaptive routine
    //specify function
    gsl_function F;
    F.function = &length_function;
    // integration bounds
    double lower_bound = 0.0;
    double upper_bound = T_parametric_max;
    // desired absolute and relative error
    double epsabs = 0.0;
    double epsrel = 1e-7;
    // maximum number of sub-intervals
    size_t limit = 1000;
    // specify integration key, GSL_INTEG_GAUSS61  (key = 6)
    int key = 6;
    gsl_integration_qag(&F, lower_bound, upper_bound, epsabs,
                        epsrel, limit, key, w, length, abserror);

return 1;

the function  I use looks like,

double length_function(double t)
    extern double Ax[FREQUENCY_LIMIT];
    extern double pie;

    int i = 0;
    double length = 0.0;
    double n, angle, sin_k, cos_k;
    double A_dot_A;

    // sum over frequency modes, i, up to FREQUENCY_LIMIT
    for (i = 0; i < FREQUENCY_LIMIT; i++) {
        n = (double) i;
        angle = 2.0 * pie * n * t / T_parametric_max;
        sin_k = sin(angle);
        cos_k = cos(angle);
        A_dot_A = Ax[i] * Ax[i];
        length += A_dot_A * sin_k * sin_k;    }

    length = sqrt(length);
    return length;

I notice in the documentation that they define the gsl_function function parameters with a statement like,
          F.params = &alpha;
I don't do that, is that why I get bus errors when I try to run the scrpt?

any suggestions would be appreciated.


Nathan Moore

reply via email to

[Prev in Thread] Current Thread [Next in Thread]