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