[SciPy-User] Scipy dblquad Calculation

[Warren Weckesser]

Scott Barlow wrote:
> 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.
>

Scott,

As you suspected, the problem is your use of arctan, and an appropriate 
use of arctan2 can correct it.  If I change this:

                           t*scipy.sin(theta + 
scipy.arctan((a**2.0/b**2.0)*(1.0/scipy.tan(theta))))


to this:

                           t*scipy.sin(theta + scipy.arctan2( 
a**2.0*scipy.cos(theta), b**2.0*scipy.sin(theta)) )

then the results agree.


Warren

> 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
> ------------------------------------------------------------------------
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>