[SciPy-user] Double integration using gaussian quadrature query

[David Unitt]

Dear SciPy Users, I would be very grateful of your help,

 

I am trying to do a double integral using Gaussian quadrature, I have put
the code at the bottom of this message, along with the output. The problem I
have is that if the function 'integrand' returns a value calculated using
'a0' the program throws up the exception "ValueError: frames are not
aligned". 'a0' starts out as the integration variable 'a' in function
'integral2'. I have traced the problem to the successive calls of function
'integral1' from 'integral2'. Each time this happens the variable 'aa' in
'integral1' gains an element. Specifying aa[0] doesn't work.

 

Any help on how to get round this would be appreciated

Thanks

david

 

--------CODE-------

#!/usr/bin/env python

 

####################################################################

#calculating a double integral with  fixed variable

#

####################################################################

 

from Numeric import *

from __future__ import division

from scipy import integrate

 

 

#######################

#define integrand 

def integrand(b0,a0,c0):

    print "integrand called"

    return a0*b0*b0*c0 #If this contains a0 the program crashes

        

#calculate the double integral

def integral1(aa,b1,b2,c):

    print "integral1 called"

    print "aa = ",aa

    x1 = integrate.quadrature(lambda b: integrand(b,aa,c), b1, b2)

    return x1[0]

 

 

def integral2(a1,a2,b1,b2,c):

    #a1, a2 are limits of this integral, b1,b2 are limits for 'integral1'

    #c is some constant

    print "integral2 called"

    x2 = integrate.quadrature(lambda a: integral1(a,b1,b2,c), a1, a2)

    return x2

      

 

##

f = integral2(0.0,1.0,0.0,1.0,1.0)

print f[0]

 

 

-----OUTPUT--------

integral2 called

integral1 called

aa =  [ 0.5+0.j]

integrand called

integrand called

integrand called

Took 4 points.

integral1 called

aa =  [ 0.21132487+0.j  0.78867513+0.j]

integrand called

integrand called

Took 3 points.

integral1 called

aa =  [ 0.11270167+0.j  0.5       +0.j  0.88729833+0.j]

integrand called

integrand called

Traceback (most recent call last):

  File "C:/quad_test2.py", line 37, in ?

    f = integral2(0.0,1.0,0.0,1.0,1.0)

  File "C:/quad_test2.py", line 32, in integral2

    x2 = integrate.quadrature(lambda a: integral1(a,b1,b2,c), a1, a2)

  File
"C:\PROGRA~1\python22\Lib\site-packages\scipy\integrate\quadrature.py", line
65, in quadrature

    newval = fixed_quad(func,a,b,args,n)[0]

  File
"C:\PROGRA~1\python22\Lib\site-packages\scipy\integrate\quadrature.py", line
35, in fixed_quad

    return (b-a)/2.0*sum(w*func(y,*args)), None

  File "C:/quad_test2.py", line 32, in <lambda>

    x2 = integrate.quadrature(lambda a: integral1(a,b1,b2,c), a1, a2)

  File "C:/quad_test2.py", line 24, in integral1

    x1 = integrate.quadrature(lambda b: integrand(b,aa,c), b1, b2)

  File
"C:\PROGRA~1\python22\Lib\site-packages\scipy\integrate\quadrature.py", line
65, in quadrature

    newval = fixed_quad(func,a,b,args,n)[0]

  File
"C:\PROGRA~1\python22\Lib\site-packages\scipy\integrate\quadrature.py", line
35, in fixed_quad

    return (b-a)/2.0*sum(w*func(y,*args)), None

  File "C:/quad_test2.py", line 24, in <lambda>

    x1 = integrate.quadrature(lambda b: integrand(b,aa,c), b1, b2)

  File "C:/quad_test2.py", line 18, in integrand

    return a0*b0*b0*c0 #If this contains a) the program crashes

ValueError: frames are not aligned



_____________________________________________________________________
This e-mail has been scanned for viruses by MCI's Internet Managed Scanning Services - powered by MessageLabs. For further information visit http://www.mci.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20030917/3369d3cd/attachment.html>