From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 22426 invoked by alias); 22 Oct 2007 18:59:28 -0000 Received: (qmail 22410 invoked by uid 22791); 22 Oct 2007 18:59:26 -0000 X-Spam-Check-By: sourceware.org Received: from fk-out-0910.google.com (HELO fk-out-0910.google.com) (209.85.128.184) by sourceware.org (qpsmtpd/0.31) with ESMTP; Mon, 22 Oct 2007 18:59:22 +0000 Received: by fk-out-0910.google.com with SMTP id 26so1278242fkx for ; Mon, 22 Oct 2007 11:59:17 -0700 (PDT) Received: by 10.82.156.12 with SMTP id d12mr9064724bue.1193079557163; Mon, 22 Oct 2007 11:59:17 -0700 (PDT) Received: by 10.82.183.15 with HTTP; Mon, 22 Oct 2007 11:59:16 -0700 (PDT) Message-ID: <8d19b8a10710221159x23d93074tfea410a9e3bddaeb@mail.gmail.com> Date: Mon, 22 Oct 2007 18:59:00 -0000 From: "Kevin Jackman" To: gsl-discuss@sourceware.org Subject: multidimensional numerical integration MIME-Version: 1.0 Content-Type: text/plain; charset=ISO-8859-1 Content-Transfer-Encoding: 7bit Content-Disposition: inline Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2007-q4/txt/msg00002.txt.bz2 Hi all, A couple of questions followed by some useful information for others (in my opinion): Q: I realize "epsabs" and "epsrel" control the desired accuracy of numerical integration, but how? Mainly, how do they each affect the GSL routines, and how do they differ from each other? How about on a Quadpack level? I think the GSL manual is kind of vague on how these parameters individually control numerical integration convergence. (No hard feelings, ok?) Some Useful Info: Looking through the GSL documentation, archives, and source code, I soon discovered that GSL does not have specific routines for multidimensional numerical integration, except using Monte Carlo methods. I saw some efforts to initially incorporate cubature routines: http://sourceware.org/ml/gsl-discuss/2007-q2/msg00046.html http://www.cygwin.com/ml/gsl-discuss/2005-q2/msg00020.html However, as cubature has not been formally incorporated into GSL, and that topic is far from my specialty, I moved on. Below is a little tidbit of code which demonstrates how to recursively integrate using 1D QAGS routines and thereby yield an effective 3D integral. The code integrates the volume of a sphere in spherical coordinates. I realize this method can be time and computer resource consuming for complex integrals, but I have yet to see a better GSL alternative suited for my purposes. (Monte Carlo is not good for my needs, long explanation). Also, the absolute error estimation may be lost by doing it this way, though I'm not sure. Aside from telling me which hand I should wish in and which hand I should..., I'd like to see something like this incorporated into the GSL manual, so others don't have to reinvent the wheel. You are welcome to use my code, modify it, etc... to make it up to your manual standards. No props necessary even. compile and build with something like: gcc -o simpleQAGS_3D simpleQAGS_3D.c -L/usr/lib -I/usr/include -lgsl -lgslcblas -lm simpleQAGS_3D.c: #include #include #include double Phi (double phi, void * params) { return 1.0; } double Theta (double theta, void * params) { gsl_integration_workspace * w3 = gsl_integration_workspace_alloc (1000); double result3, error3; double expected3 = 2.0*M_PI; double alpha3 = 1.0; double TwoPi = 2.0*M_PI ; gsl_function F3; F3.function = Φ F3.params = &alpha3; gsl_integration_qags (&F3, 0, TwoPi, 1e-5, 1e-7, 1000, w3, &result3, &error3); printf ("W3 result = %.18f\n", result3); printf ("W3 expected result = %.18f\n", expected3); printf ("W3 estimated error = %.18f\n", error3); printf ("W3 actual error = %.18f\n", result3 - expected3); printf ("W3 intervals = %d\n", w3->size); gsl_integration_workspace_free (w3); return sin(theta)*result3; } double R (double r, void * params) { gsl_integration_workspace * w2 = gsl_integration_workspace_alloc (1000); double result2, error2; double expected2 = M_PI; double alpha2 = 1.0; double Pi = M_PI; gsl_function F2; F2.function = Θ F2.params = &alpha2; gsl_integration_qags (&F2, 0, Pi, 1e-5, 1e-7, 1000, w2, &result2, &error2); printf ("W2 result = %.18f\n", result2); printf ("W2 expected result = %.18f\n", expected2); printf ("W2 estimated error = %.18f\n", error2); printf ("W2 actual error = %.18f\n", result2 - expected2); printf ("W2 intervals = %d\n", w2->size); gsl_integration_workspace_free (w2); return pow(r,2.0)*result2; } int main (void) { gsl_integration_workspace * w1 = gsl_integration_workspace_alloc (1000); double result1, error1; double alpha1 = 1.0; double radius = 2.0; double expected1 = (4.0*M_PI*pow(radius,3.0)) / 3.0; gsl_function F1; F1.function = &R; F1.params = &alpha1; gsl_integration_qags (&F1, 0, radius, 1e-5, 1e-7, 1000, w1, &result1, &error1); printf ("\nW1 result = %.18f\n", result1); printf ("W1 expected result = %.18f\n", expected1); printf ("W1 estimated error = %.18f\n", error1); printf ("W1 actual error = %.18f\n", result1 - expected1); printf ("W1 intervals = %d\n", w1->size); gsl_integration_workspace_free (w1); return 0; } Cheers, Kevin