Ack! How to compute double integral?

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Ack! How to compute double integral?

Joe Reinhardt

I am trying to compute an M-dimensional integral using quad().   The
documentation says quad only does integration over a single variable,
so I tried nesting subroutine calls over each of my variables.

This method doesn't seem to work.

I attached some sample code for 2-D showing that my technique fails.
The outer quad() call only evaluates one value of "y".

1) Is there a way to compute M-dimensional integrals in Octave?

2) I used global variables to pass parameters from integrand to
   integrand.   Is there a better way to do this?

Thanks!




--
Joe Reinhardt                            [hidden email]
Department of Biomedical Engineering     1402 Engineering Building
University of Iowa                       Iowa City, IA 52242
Phone: 319-335-5634                      Fax: 319-335-5631

test.m (570 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Ack! How to compute double integral?

John W. Eaton-6
On  3-Nov-1997, Joe Reinhardt <[hidden email]> wrote:

| I am trying to compute an M-dimensional integral using quad().   The
| documentation says quad only does integration over a single variable,
| so I tried nesting subroutine calls over each of my variables.
|
| This method doesn't seem to work.

It won't work because Octave's quad function is ultimately implemented
with some Fortran subroutines that can't be called recursively.

| I attached some sample code for 2-D showing that my technique fails.
| The outer quad() call only evaluates one value of "y".
|
| 1) Is there a way to compute M-dimensional integrals in Octave?

Yes, you can use an M-file specifically designed for multidimensional
integration.  Here is some information that was posted to this list
some time ago:

  Date: Wed, 29 Nov 95 07:49:41 CST
  From: [hidden email] (Vinayak Dutt)
  Subject: Re: How to quad("sin(a*sin(x))*sin(x)",0,pi/2)?
  To: [hidden email]

  Francesco Potorti wrote:
  #
  # I found the following in the matlab archives, form the
  # contribution/integration directory or some like that, if I am not
  # wrong.  They should be public domain.  I also have the sources.
  #
  # ===================================================================
  # Index for integration
  # ===================================================================
  #
  # nit: <directory>
  #
  # ===================================================================
  #
  #  Numerical Integration Toolbox
  #  
  #  MATLAB Toolbox for 1-D, 2-D, and n-D Numerical Integration
  #  
  #  The original 1-D routines were obtained from NETLIB and were
  #  written by
  #           Howard Wilson
  #           Department of Engineering Mechanics
  #           University of Alabama
  #           Box 870278
  #           Tuscaloosa, Alabama 35487-0278
  #           Phone 205 348-1617
  #           Email address: HWILSON @ UA1VM.UA.EDU
  #  
  #  The rest of the routines were written by
  #           Bryce Gardner
  #           Ray W. Herrick Laboratories
  #           Purdue University
  #           West Lafayette, IN 47906
  #           Phone: 317-494-0231
  #           Fax:  317-494-0787
  #           Email:  [hidden email]
  #  
  #  These are the general purpose integration routines:
  #  
  #         quadg.m      -- High accuracy replacement for QUAD and QUAD8 (1-D)
  #         quad2dg.m    -- 2-D integration over a rectangular region
  #         quad2dggen.m -- 2-D integration over a general region
  #         quadndg.m    -- n-D integration over a n-D hyper-rectangular region

  I had fixed these routines sometime back (to work well in Octave)
  and have used them myself and found them quite useful. Gaussian
  Quadrature works a lot faster too :-)

  I had checked with Gardner some time back and he did not have any
  problems if Octave ftp sites kept his software.

  If anyone is interested, I can mail the fixed version of that MATLAB
  package.

I don't remember whether Vinayak ever made his versions available.
Perhaps something like these functions will eventually be distributed
as part of Octave.  Would someone like to volunteer to make that
happen?

| 2) I used global variables to pass parameters from integrand to
|    integrand.   Is there a better way to do this?

Nope, not with quad.  I don't know about the set of M-files for
quadrature mentioned in the message above.  They may allow for passing
additional parameters.

Thanks,

jwe


Reply | Threaded
Open this post in threaded view
|

Re: Ack! How to compute double integral?

Joe Reinhardt
"John W. Eaton" <[hidden email]> writes:

> Yes, you can use an M-file specifically designed for multidimensional
> integration.  Here is some information that was posted to this list
> some time ago:
>
>   Date: Wed, 29 Nov 95 07:49:41 CST
>   From: [hidden email] (Vinayak Dutt)
>   Subject: Re: How to quad("sin(a*sin(x))*sin(x)",0,pi/2)?
>   To: [hidden email]

Thanks for the info.

I can't connect to the http://www.tsc.uvigo.es/GTS/Octave archive
site.  Is there a mirror site somewhere?

I also tried directly contacting Vinayak Dutt, but the message bounced
back from Mayo saying his mailbox was closed.

Thanks for any help -

Joe



--
Joe Reinhardt                            [hidden email]
Department of Biomedical Engineering     1402 Engineering Building
University of Iowa                       Iowa City, IA 52242
Phone: 319-335-5634                      Fax: 319-335-5631


Reply | Threaded
Open this post in threaded view
|

Re: Ack! How to compute double integral?

John W. Eaton-6
On  4-Nov-1997, Joe Reinhardt <[hidden email]> wrote:

| "John W. Eaton" <[hidden email]> writes:
|
| > Yes, you can use an M-file specifically designed for multidimensional
| > integration.  Here is some information that was posted to this list
| > some time ago:
| >
| >   Date: Wed, 29 Nov 95 07:49:41 CST
| >   From: [hidden email] (Vinayak Dutt)
| >   Subject: Re: How to quad("sin(a*sin(x))*sin(x)",0,pi/2)?
| >   To: [hidden email]
|
| Thanks for the info.
|
| I can't connect to the http://www.tsc.uvigo.es/GTS/Octave archive
| site.  Is there a mirror site somewhere?

I just connected (so it appears that it still exists) but the
connection was very slow.  I don't know of any mirrors of that site,
and I don't have time to set up and maintain one myself.  Would
someone else be willing to do that?

| I also tried directly contacting Vinayak Dutt, but the message bounced
| back from Mayo saying his mailbox was closed.

Hmm.  If [hidden email] doesn't work, I don't know how to reach
him.

Thanks,

jwe