Scipy dblquad Calculation
![](https://secure.gravatar.com/avatar/6fe9a44b78807ee86f4ca62a695d2977.jpg?s=120&d=mm&r=g)
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
![](https://secure.gravatar.com/avatar/b0f62d137f9ea1d0b6cc4e7e6f61b119.jpg?s=120&d=mm&r=g)
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
![](https://secure.gravatar.com/avatar/6fe9a44b78807ee86f4ca62a695d2977.jpg?s=120&d=mm&r=g)
Warren Weckesser <warren.weckesser <at> enthought.com> writes:
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
Ah, excellent. Thank you very much for your help Warren! When I tried using arctan2, I didn't think about breaking up tangent into its sine and cosine components and instead kept it as tangent or cotangent in the numerator or denominator, as appropriate. Thanks, Scott
participants (3)
-
Scott Barlow
-
Scott Barlow
-
Warren Weckesser