[SciPy-User] Scipy dblquad Calculation

[Scott Barlow]

Hello all.  I'm writing some code and need to compute a double integral.
I'm trying to calculate this using dblquad, but I'm having an issue with the
result.  I get a result and no error messages, but the result is not as
expected.  I've check the code over and over and I see no mistakes, so I'm
asking for suggestions here (or reasons for the difference).  I realize that
the difference could be from (quadrature) estimations, but I don't think
that is it.

Python version: 2.6.4
SciPy version: 0.7.0
NumPy version: 1.3.0

Please see the code below.  This describes the equation for the bending
rigidity of an elliptical composite laminate of 16 plies under bending.
Note that "a" and "b" would be the length in inches of the major and minor
axes, and I have them both set to 2.0 here, thus this is actually a circle.
The bottom equation (Final_Result_2) is the reduced equation for a circular
tube (no integration required).  Thus when a=b, the result between the two
equations should be the same.

When run, I get the following:
21560941.3345   (Final_Result)
22852651.614    (Final_Result_2)

As mentioned before, I've checked the equations numerous times for typos,
but I haven't found any.  I'm quite sure the equation is input correctly.
I've also tried using arctan2 instead of arctan, but I get the same result.

Here is the python code:
---------------------------------------------------------------------------
-------------
import scipy
import numpy
from scipy.integrate import dblquad

Q11 = 20009576.0
Q22 = 1397653.8
Qbar11 = [Q11, Q11,  Q11,  Q11,  Q22,  Q22,  Q22,  Q22, Q22, Q22, Q22, Q22,
Q11, Q11, Q11, Q11];

a = 2.0;
b = 2.0;
thickness = 0.005;

Final_Result = 0.0;
for i in range(1, 17):
  Temp_Result = dblquad(lambda theta, t:
                        Qbar11[i-1] *
                        scipy.sin(theta)**2.0 *
                        (
                           1.0*((a**2.0 * b**2.0)/((a**2.0 -
b**2.0)*scipy.sin(theta)**2.0 + b**2.0)) +
                           t**2.0 +
                           2.0*((a**2.0 * b**2.0)/((a**2.0 -
b**2.0)*scipy.sin(theta)**2.0 + b**2.0))**(0.5)*
                           t*scipy.sin(theta +
scipy.arctan((a**2.0/b**2.0)*(1.0/scipy.tan(theta))))
                        )**(3.0/2.0),
                        thickness*(i-1),
                        thickness*(i),
                        lambda x: 0.0,
                        lambda x: 2.0*scipy.pi,
                        epsabs = 1.4899999999999999e-12,
                        epsrel = 1.4899999999999999e-12);
  Final_Result += Temp_Result[0];
print Final_Result;

ri = a;
Final_Result_2 = 0.0;
for i in range(1, 17):
  Final_Result_2 += (scipy.pi/4.0) * Qbar11[i-1] * ((ri +
thickness*(i))**4.0 - (ri + thickness*(i-1))**4.0);
print Final_Result_2;
---------------------------------------------------------------------------
-------------

Thanks for any help,
Scott
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100223/3dd53831/attachment.html>