Databases
StackOverflow2013
Postsdbo
PONumerical Integration over a Matrix of Functions, SymPy and SciPy
 

Note that there are some explanatory texts on larger screens.

plurals
  1. PONumerical Integration over a Matrix of Functions, SymPy and SciPy
    primarykey
    data
    text
    <p>From my SymPy output I have the matrix shown below, which I must integrate in 2D. Currently I am doing it element-wise as shown below. This method works but it gets too slow (for both <code>sympy.mpmath.quad</code> and <code>scipy.integrate.dblquad</code>) for my real case (in which <code>A</code> and its functions are much bigger (see edit below):</p> <pre><code>from sympy import Matrix, sin, cos import sympy import scipy sympy.var( 'x, t' ) A = Matrix([[(sin(2-0.1*x)*sin(t)*x+cos(2-0.1*x)*cos(t)*x)*cos(3-0.1*x)*cos(t)], [(cos(2-0.1*x)*sin(t)*x+sin(2-0.1*x)*cos(t)*x)*sin(3-0.1*x)*cos(t)], [(cos(2-0.1*x)*sin(t)*x+cos(2-0.1*x)*sin(t)*x)*sin(3-0.1*x)*sin(t)]]) # integration intervals x1,x2,t1,t2 = (30, 75, 0, 2*scipy.pi) # element-wise integration from sympy.utilities import lambdify from sympy.mpmath import quad from scipy.integrate import dblquad A_int1 = scipy.zeros( A.shape, dtype=float ) A_int2 = scipy.zeros( A.shape, dtype=float ) for (i,j), expr in scipy.ndenumerate(A): tmp = lambdify( (x,t), expr, 'math' ) A_int1[i,j] = quad( tmp, (x1, x2), (t1, t2) ) # or (in scipy) A_int2[i,j] = dblquad( tmp, t1, t2, lambda x:x1, lambda x:x2 )[0] </code></pre> <p>I was considering doing it in one shot like, but I'm not sure if this is the way to go:</p> <pre><code>A_eval = lambdify( (x,t), A, 'math' ) A_int1 = sympy.quad( A_eval, (x1, x2), (t1, t2) # or (in scipy) A_int2 = scipy.integrate.dblquad( A_eval, t1, t2, lambda x: x1, lambda x: x2 )[0] </code></pre> <hr> <p>EDIT: The real case has been made available in <a href="https://www.dropbox.com/s/eyh4zr8ee9lgq84/numerical-integration-over-a-matrix-of-functions-sympy-and-scipy.zip" rel="nofollow noreferrer">this link</a>. Just unzip and run <code>shadmehri_2012.py</code> (is the author from were this example was taken from: <a href="http://www.sciencedirect.com/science/article/pii/S0263822311003606" rel="nofollow noreferrer">Shadmehri et al. 2012</a>). I've started a bounty of 50 for the one who can do the following:</p> <ul> <li>make it reasonably faster than the proposed question<br/></li> <li>manage to run without giving memory error even with a number of terms <code>m=15</code> and <code>n=15</code> in the code), I managed up to <code>m=7</code> and <code>n=7</code> in 32-bit</li> </ul> <p>The current timing can be summarized below(measured with m=3 and n=3). From there it can be seen that the numerical integration is the bottleneck.</p> <p>build trial functions = 0%<br/> evaluating differential equations = 2%<br/> lambdifying k1 = 22%<br/> integrating k1 = 74%<br/> lambdifying and integrating k2 = 2%<br/> extracting eigenvalues = 0%<br/></p> <hr> <p>Related questions: <a href="https://stackoverflow.com/a/10683911/832621">about lambdify</a></p>
    singulars
    1. PO
      singulars
      1. PTAnswer
    1. This table or related slice is empty.
    1. PTQuestion
    1. USCommunity
    1. USSaullo G. P. Castro
    plurals
    1. PL
      singulars
      1. LTLinked
    2. PL
      singulars
      1. LTLinked
    1. PL
      singulars
      1. LTLinked
    2. PL
      singulars
      1. LTLinked
    1. This table or related slice is empty.
    1. PO
      singulars
      1. PTAnswer
    1. VO
      singulars
      1. PONumerical Integration over a Matrix of Functions, SymPy and SciPy
      1. This table or related slice is empty.
      1. VTUpMod
    2. VO
      singulars
      1. PONumerical Integration over a Matrix of Functions, SymPy and SciPy
      1. USSaullo G. P. Castro
      1. VTBountyStart
    3. VO
      singulars
      1. PONumerical Integration over a Matrix of Functions, SymPy and SciPy
      1. This table or related slice is empty.
      1. VTUpMod
    1. COAnd profiling your code with `line_profiler` http://silas.sewell.org/blog/2009/05/28/python-line-by-line-profiler-line_profiler-and-kernprof/ would have helped you.
      singulars
      1. PONumerical Integration over a Matrix of Functions, SymPy and SciPy
      1. USKrastanov
    2. COThe problem is that A is obtained after applying an automatic differentiation over another matrix, using the procedure described [here](http://stackoverflow.com/a/16237073/832621), explaining why it comes from SymPy. I fully agree with you to use directly the numeric modules like `numpy` when applicable.
      singulars
      1. PONumerical Integration over a Matrix of Functions, SymPy and SciPy
      1. USSaullo G. P. Castro
    3. COOh, my bad... And looking closely at the code it seems my issue nb. 1 is also incorrect, as you are not redoing the `lambdify` call each time. So in the hope of actually being helpful unlike in the previous comment: have you actually tried to do the integration symbolically with numpy. If it is only `sin`, `cos` and polynomials it should be easy. Especially if you can permit yourself to rewrite it as exponentials `expr.rewrite(exp).expand()` (there are complex numbers popping out, but after evaluation it is real).
      singulars
      1. PONumerical Integration over a Matrix of Functions, SymPy and SciPy
      1. USKrastanov
 

Querying!

 
Guidance
A row detail

Detail views are divided into sections. All the information in the data section comes from columns in the selected row. The other sections display data from other, related rows.

Related data can be related in a to-one or a to-many fashion. Captions of data related in a to-many fashion link to a list view showing a filtered view of the table.

Try moving around until you find a non-empty to-many entry and click on the label to get to one. You can move back to the root by clicking on the database name in the header.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload