[SciPy-User] Scipy dblquad Calculation
Warren Weckesser
warren.weckesser@enthought....
Wed Feb 24 05:16:52 CST 2010
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@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
More information about the SciPy-User
mailing list