On Tue, Jun 29, 2010 at 12:31 PM, Colin <[email protected]> wrote: > Dear Sympy list, > I'm trying to evaluate a numerical dispersion relation for a finite > element method we are working on, which requires construction of a 4x4 > Hermitian matrix A which is a function of l,k,w, and evaluation of > it's determinant. This gives a quartic polynomial which can be solved > for w for each l,k. I want to do this symbolically because then I can > trace the various root branches - if I evaluate A numerically for > various l,k and solve for w it is not so easy to trace the branches (I > could use a continuation method, though). I've tried evaluating the > determinant using various Sympy methods, including taking the LU > decomposition and looking at the diagonals. All of these methods tend > to blow up leading to very long computations which prevent progress. I > was just wondering if anyone on the list has any tips on possible ways > to simplify the problem?
Very interesting. Me, Andy and Mateusz happen to be at this finite element conference right now: http://hpfem.org/events/esco-2010/ I tried your script, here is the output: $ python a.py Checking local coordinates in lower triangle Checking 1D Lagrange polynomials Checking local basis functions in lower triangle Checking lower triangle integral Constructing mass matrix Checking mass matrix sums to area. Computing derivatives of basis functions Constructing Laplacian matrix checking Laplacian matrix row sums to zero Constructing x-derivative matrix Constructing y-derivative matrix checking x-derivative matrix row sums to zero checking y-derivative matrix row sums to zero constructing s matrix Forming global mass matrix Forming global laplace matrix So it finished in about 11s. How long does it take for you? I then printed the matrix A and got: [ -120 + 4*w, w*(exp(-I*k/4)/3 + exp(I*k/4)/3 + exp(-I*k/4 + I*l*3**(1/2)/4)/3 + exp(I*k/8 + I*l*3**(1/2)/8)/3 + exp(I*k/4 - I*l*3**(1/2)/4)/3 + exp(-I*k/8 - I*l*3**(1/2)/8)/3) + (4860 + 4860*exp(I*k/2) + 4860*exp(I*k/4)*exp(-I*k/4 + I*l*3**(1/2)/4) + 4860*exp(I*k/4)*exp(I*k/8 + I*l*3**(1/2)/8) + 4860*exp(I*k/4)*exp(I*k/4 - I*l*3**(1/2)/4) + 4860*exp(I*k/4)*exp(-I*k/8 - I*l*3**(1/2)/8))*exp(-I*k/4)/729, w*(exp(-I*k/4)/3 + exp(I*k/4)/3 + exp(-I*k/8 + I*l*3**(1/2)/8)/3 + exp(-I*k/4 - I*l*3**(1/2)/4)/3 + exp(I*k/8 - I*l*3**(1/2)/8)/3 + exp(I*k/4 + I*l*3**(1/2)/4)/3) + (4860 + 4860*exp(I*k/2) + 4860*exp(I*k/4)*exp(-I*k/8 + I*l*3**(1/2)/8) + 4860*exp(I*k/4)*exp(-I*k/4 - I*l*3**(1/2)/4) + 4860*exp(I*k/4)*exp(I*k/8 - I*l*3**(1/2)/8) + 4860*exp(I*k/4)*exp(I*k/4 + I*l*3**(1/2)/4))*exp(-I*k/4)/729, (9720 + 9720*exp(I*k) + 2430*exp(I*k/2)*exp(I*k/2 + I*l*3**(1/2)/4) + 2430*exp(I*k/2)*exp(I*k/2 - I*l*3**(1/2)/4) + 2430*exp(I*k/2)*exp(-I*k/2 - I*l*3**(1/2)/4) + 2430*exp(I*k/2)*exp(-I*k/2 + I*l*3**(1/2)/4))*exp(-I*k/2)/729 - w*(432*exp(I*l*3**(1/2)) + 432*exp(I*l*3**(1/2)/2) + 864*exp(I*l*3**(1/2)/4) + 864*exp(5*I*l*3**(1/2)/4))*exp(-3*I*l*3**(1/2)/4)/5184] [w*(exp(-I*k/4)/3 + exp(I*k/4)/3 + exp(-I*k/4 + I*l*3**(1/2)/4)/3 + exp(I*k/8 + I*l*3**(1/2)/8)/3 + exp(I*k/4 - I*l*3**(1/2)/4)/3 + exp(-I*k/8 - I*l*3**(1/2)/8)/3) + (4860 + 4860*exp(I*k/2) + 4860*exp(I*k/4)*exp(-I*k/4 + I*l*3**(1/2)/4) + 4860*exp(I*k/4)*exp(I*k/8 + I*l*3**(1/2)/8) + 4860*exp(I*k/4)*exp(I*k/4 - I*l*3**(1/2)/4) + 4860*exp(I*k/4)*exp(-I*k/8 - I*l*3**(1/2)/8))*exp(-I*k/4)/729, -120 + 4*w, w*(exp(I*k/2)/3 + exp(-I*k/2)/3 + exp(-I*k/8 + I*l*3**(1/2)/8)/3 + exp(I*k/8 + I*l*3**(1/2)/8)/3 + exp(-I*k/8 - I*l*3**(1/2)/8)/3 + exp(I*k/8 - I*l*3**(1/2)/8)/3) + (4860 + 4860*exp(I*k) + 4860*exp(I*k/2)*exp(-I*k/8 + I*l*3**(1/2)/8) + 4860*exp(I*k/2)*exp(I*k/8 + I*l*3**(1/2)/8) + 4860*exp(I*k/2)*exp(-I*k/8 - I*l*3**(1/2)/8) + 4860*exp(I*k/2)*exp(I*k/8 - I*l*3**(1/2)/8))*exp(-I*k/2)/729, w*(-exp(-3*I*k/8 + I*l*3**(1/2)/8)/12 - exp(3*I*k/8 - I*l*3**(1/2)/8)/12 - exp(3*I*k/4 - I*l*3**(1/2)/4)/6 - exp(-3*I*k/4 + I*l*3**(1/2)/4)/6) + 10*exp(-5*I*k/8 - I*l*3**(1/2)/8)/3 + 10*exp(5*I*k/8 + I*l*3**(1/2)/8)/3 + 10*exp(I*k/8 - 3*I*l*3**(1/2)/8)/3 + 40*exp(-I*k/4 - I*l*3**(1/2)/4)/3 + 40*exp(I*k/4 + I*l*3**(1/2)/4)/3 + 10*exp(-I*k/8 + 3*I*l*3**(1/2)/8)/3] [w*(exp(-I*k/4)/3 + exp(I*k/4)/3 + exp(-I*k/8 + I*l*3**(1/2)/8)/3 + exp(-I*k/4 - I*l*3**(1/2)/4)/3 + exp(I*k/8 - I*l*3**(1/2)/8)/3 + exp(I*k/4 + I*l*3**(1/2)/4)/3) + (4860 + 4860*exp(I*k/2) + 4860*exp(I*k/4)*exp(-I*k/8 + I*l*3**(1/2)/8) + 4860*exp(I*k/4)*exp(-I*k/4 - I*l*3**(1/2)/4) + 4860*exp(I*k/4)*exp(I*k/8 - I*l*3**(1/2)/8) + 4860*exp(I*k/4)*exp(I*k/4 + I*l*3**(1/2)/4))*exp(-I*k/4)/729, w*(exp(I*k/2)/3 + exp(-I*k/2)/3 + exp(-I*k/8 + I*l*3**(1/2)/8)/3 + exp(I*k/8 + I*l*3**(1/2)/8)/3 + exp(-I*k/8 - I*l*3**(1/2)/8)/3 + exp(I*k/8 - I*l*3**(1/2)/8)/3) + (4860 + 4860*exp(I*k) + 4860*exp(I*k/2)*exp(-I*k/8 + I*l*3**(1/2)/8) + 4860*exp(I*k/2)*exp(I*k/8 + I*l*3**(1/2)/8) + 4860*exp(I*k/2)*exp(-I*k/8 - I*l*3**(1/2)/8) + 4860*exp(I*k/2)*exp(I*k/8 - I*l*3**(1/2)/8))*exp(-I*k/2)/729, -120 + 4*w, w*(-exp(-3*I*k/8 - I*l*3**(1/2)/8)/12 - exp(3*I*k/4 + I*l*3**(1/2)/4)/6 - exp(3*I*k/8 + I*l*3**(1/2)/8)/12 - exp(-3*I*k/4 - I*l*3**(1/2)/4)/6) + 10*exp(-5*I*k/8 + I*l*3**(1/2)/8)/3 + 40*exp(-I*k/4 + I*l*3**(1/2)/4)/3 + 10*exp(5*I*k/8 - I*l*3**(1/2)/8)/3 + 10*exp(-I*k/8 - 3*I*l*3**(1/2)/8)/3 + 40*exp(I*k/4 - I*l*3**(1/2)/4)/3 + 10*exp(I*k/8 + 3*I*l*3**(1/2)/8)/3] [ (9720 + 9720*exp(I*k) + 2430*exp(I*k/2)*exp(I*k/2 + I*l*3**(1/2)/4) + 2430*exp(I*k/2)*exp(I*k/2 - I*l*3**(1/2)/4) + 2430*exp(I*k/2)*exp(-I*k/2 - I*l*3**(1/2)/4) + 2430*exp(I*k/2)*exp(-I*k/2 + I*l*3**(1/2)/4))*exp(-I*k/2)/729 - w*(432*exp(I*l*3**(1/2)) + 432*exp(I*l*3**(1/2)/2) + 864*exp(I*l*3**(1/2)/4) + 864*exp(5*I*l*3**(1/2)/4))*exp(-3*I*l*3**(1/2)/4)/5184, w*(-exp(-3*I*k/8 + I*l*3**(1/2)/8)/12 - exp(3*I*k/8 - I*l*3**(1/2)/8)/12 - exp(3*I*k/4 - I*l*3**(1/2)/4)/6 - exp(-3*I*k/4 + I*l*3**(1/2)/4)/6) + 10*exp(-5*I*k/8 - I*l*3**(1/2)/8)/3 + 10*exp(5*I*k/8 + I*l*3**(1/2)/8)/3 + 10*exp(I*k/8 - 3*I*l*3**(1/2)/8)/3 + 40*exp(-I*k/4 - I*l*3**(1/2)/4)/3 + 40*exp(I*k/4 + I*l*3**(1/2)/4)/3 + 10*exp(-I*k/8 + 3*I*l*3**(1/2)/8)/3, w*(-exp(-3*I*k/8 - I*l*3**(1/2)/8)/12 - exp(3*I*k/4 + I*l*3**(1/2)/4)/6 - exp(3*I*k/8 + I*l*3**(1/2)/8)/12 - exp(-3*I*k/4 - I*l*3**(1/2)/4)/6) + 10*exp(-5*I*k/8 + I*l*3**(1/2)/8)/3 + 40*exp(-I*k/4 + I*l*3**(1/2)/4)/3 + 10*exp(5*I*k/8 - I*l*3**(1/2)/8)/3 + 10*exp(-I*k/8 - 3*I*l*3**(1/2)/8)/3 + 40*exp(I*k/4 - I*l*3**(1/2)/4)/3 + 10*exp(I*k/8 + 3*I*l*3**(1/2)/8)/3, -90 - 5*exp(I*k) - 5*exp(-I*k) - w*(131072 - 2359296*exp(I*k) + 131072*exp(2*I*k) + 131072*exp(I*k)*exp(I*k/2 + I*l*3**(1/2)/2) + 131072*exp(I*k)*exp(I*k/2 - I*l*3**(1/2)/2) + 131072*exp(I*k)*exp(-I*k/2 + I*l*3**(1/2)/2) + 131072*exp(I*k)*exp(-I*k/2 - I*l*3**(1/2)/2))*exp(-I*k)/1048576 - 5*exp(I*k/2 + I*l*3**(1/2)/2) - 5*exp(I*k/2 - I*l*3**(1/2)/2) - 5*exp(-I*k/2 + I*l*3**(1/2)/2) - 5*exp(-I*k/2 - I*l*3**(1/2)/2)] But I don't know how to simplify this. In my experience, if you need to invert some matrices, different methods give different results, I think reasonably simple results gives the ADJ method, see the docstring of the .inv() method. If you happen to know some better method, try to implement it. Just look into sympy/matrices/matrices.py, it's not difficult. Ondrej -- You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to [email protected]. To unsubscribe from this group, send email to [email protected]. For more options, visit this group at http://groups.google.com/group/sympy?hl=en.
