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?
Please let me know if you have any ideas or suggestions. I've appended
the script which constructs A to the bottom of this email.
all the best
--Colin
from sympy import *
x = Symbol('x',real=True)
y = Symbol('y',real=True)
#transformation (xp,yp) -> (x,y) takes
#right-angled triangle to square
#(x)= (1 1/2)(x')
#(y)= (0 sqrt(3/4))(y')
#so
#(x')= sqrt(4/3)(sqrt(3/4) -1/2)(x) = (1 -sqrt(1/3))(x)
#(y')= (0 1)(y) = (0 sqrt(4/3))(y)
# (1 1/2)(1 -sqrt(1/3))=(1 0)
# (0 sqrt(3/4))(0 sqrt(4/3)) (0 1)
xp = x - sqrt(Rational(1,3))*y
yp = sqrt(Rational(4,3))*y
#Local coordinates in bottom equilateral triangle (arrow up)
# ordering:
# 2
# 0 1
#constructed from local coordinates in right-angled triangle
# ordering:
# 2
# 0 1
B0 = 1-(xp+yp)
B1 = xp
B2 = yp
B = [B0,B1,B2]
def check_val(fn,xi,yi,val):
vali = fn._eval_subs(x,xi)._eval_subs(y,yi)
assert(vali==val)
print 'Checking local coordinates in lower triangle'
#Vertex coordinates
x0 = 0
y0 = 0
x1 = 1
y1 = 0
x2 = Rational(1,2)
y2 = sqrt(Rational(3,4))
X = [x0,x1,x2]
Y = [y0,y1,y2]
delta = eye(3)
for i in range(3):
for j in range(3):
check_val(B[i],X[j],Y[j],delta[i,j])
#Lagrange polynomials
def L1_0(x):
return 1-x
def L1_1(x):
return x
def L2_0(x):
return (2*x-1)*(x-1)
def L2_1(x):
return 4*x*(1-x)
def L2_2(x):
return (2*x-1)*x
print "Checking 1D Lagrange polynomials"
assert(L1_0(0)==1)
assert(L1_0(1)==0)
assert(L2_0(0)==1)
assert(L2_0(Rational(1,2))==0)
assert(L2_0(1)==0)
assert(L2_1(0)==0)
assert(L2_1(Rational(1,2))==1)
assert(L2_1(1)==0)
assert(L2_2(0)==0)
assert(L2_2(Rational(1,2))==0)
assert(L2_2(1)==1)
Nl = []
Nh = []
for a in range(6):
Nl.append(Symbol('Nl'+str(a)))
Nh.append(Symbol('Nh'+str(a)))
#Nl indicates local basis functions with values in lower triangle
#Ordering:
# 5
# 3 4
# 0 1 2
# 2
# 0 1
Nl[0] = L2_2(B0)
Nl[1] = B1*B0*4
Nl[2] = L2_2(B1)
Nl[3] = B2*B0*4
Nl[4] = B1*B2*4
Nl[5] = L2_2(B2)
print "Checking local basis functions in lower triangle"
x0 = 0
y0 = 0
x1 = Rational(1,2)
y1 = 0
x2 = 1
y2 = 0
x3 = Rational(1,4)
y3 = Rational(1,2)*sqrt(Rational(3,4))
x4 = Rational(1,2) + Rational(1,4)
y4 = y3
x5 = Rational(1,2)
y5 = sqrt(Rational(3,4))
X = [x0,x1,x2,x3,x4,x5]
Y = [y0,y1,y2,y3,y4,y5]
delta = eye(6)
for i in range(6):
for j in range(6):
check_val(Nl[i],X[j],Y[j],delta[i,j])
#Integral operators
def int_l(F):
"""
Function to integrate over lower triangle.
Divides into two right-angled triangles
"""
I1 = integrate(integrate(F,(y,0,sqrt(3)*x)),(x,0,Rational(1,2)))
I2 = integrate(integrate(F,(y,0,sqrt(3)-sqrt(3)*x)),
(x,Rational(1,2),1))
return I1+I2
print "Checking lower triangle integral"
assert(int_l(1)==Rational(1,2)*sqrt(Rational(3,4)))
#Construct mass matrix
M = Matrix(6,6, lambda i,j: 0)
print "Constructing mass matrix"
for a in range(6):
for b in range(6):
M[a,b] = int_l(Nl[a]*Nl[b])
Mval = Matrix(6,6, lambda i,j: 0)
for a in range(6):
for b in range(6):
Mval[a,b] = M[a,b].evalf()
print "Checking mass matrix sums to area."
Area = 0
for a in range(6):
for b in range(6):
Area = Area + M[a,b]
assert(Area==0.5*sqrt(Rational(3,4)))
#Construct Laplacian matrix
L = Matrix(6,6, lambda i,j: 0)
Nl_x = [0,0,0,0,0,0]
Nl_y = [0,0,0,0,0,0]
print "Computing derivatives of basis functions"
for a in range(6):
Nl_x[a] = diff(Nl[a],x)
Nl_y[a] = diff(Nl[a],y)
print "Constructing Laplacian matrix"
for a in range(6):
for b in range(6):
L[a,b] = int_l(Nl_x[a]*Nl_x[b] + Nl_y[a]*Nl_y[b])
print "checking Laplacian matrix row sums to zero"
for a in range(6):
Row = 0
for b in range(6):
Row += L[a,b]
assert(Row==0)
print "Constructing x-derivative matrix"
Dx = Matrix(6,6, lambda i,j: 0)
for a in range(6):
for b in range(6):
Dx[a,b] = int_l(Nl[a]*Nl_x[b])
print "Constructing y-derivative matrix"
Dy = Matrix(6,6, lambda i,j: 0)
for a in range(6):
for b in range(6):
Dy[a,b] = int_l(Nl[a]*Nl_y[b])
print "checking x-derivative matrix row sums to zero"
for a in range(6):
Row = 0
for b in range(6):
Row += Dx[a,b]
assert(Row==0)
print "checking y-derivative matrix row sums to zero"
for a in range(6):
Row = 0
for b in range(6):
Row += Dy[a,b]
assert(Row==0)
print "constructing s matrix"
k1 = Symbol('k',real=True,positive=True)
k2 = Symbol('l',real=True,positive=True)
k = Matrix(2,1,(k1,k2))
w = Symbol('w',real=True)
def e(theta):
return Matrix(2,1,(cos(theta),sin(theta)))
# 3 3
# 1 2 -> 2 0
#3 0 3 3 1 3
def myexpI(t):
#return cos(t) + I*sin(t)
return exp(I*t)
def Sn(n):
h = Rational(1,2)
H = sqrt(3)/4
S = Matrix(6,4, lambda i,j : 0)
S[0,3] = 1
S[2,3] = myexpI(k.dot(e(n*pi/3)))
S[5,3]= myexpI(k.dot(e((n+1)*pi/3)))
S[1,n%3] = myexpI(h*k.dot(e(n*pi/3)))
S[3,(n+1)%3] = myexpI(h*k.dot(e((n+1)*pi/3)))
S[4,(n+2)%3] = myexpI(H*k.dot(e((n+h)*pi/3)))
return S
print "Forming global mass matrix"
Mk = Matrix(4,4,lambda i,j: 0)
for n in range(6):
S = Sn(n)
Mk = Mk + S.H*M*S
print "Forming global laplace matrix"
Lk = Matrix(4,4,lambda i,j: 0)
for n in range(6):
S = Sn(n)
Lk = Lk + S.H*L*S
w = Symbol('w',real=True)
#P = (w*Mk - Lk).det()
#P = P.as_real_imag()
for a in range(4):
for b in range(4):
Mk[a,b] = simplify(15*Mk[a,b]/(3**Rational(1,2)))
Lk[a,b] = simplify(15*Lk[a,b]/(3**Rational(1,2)))
A = w*Mk - Lk
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to sy...@googlegroups.com.
To unsubscribe from this group, send email to
sympy+unsubscr...@googlegroups.com.
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
