Here's a short comment: using Idx(n+1) seems to work a little better for me,
since I've noticed that if I use Idx(n) with a very complicated expression,
there is a possibility of terms with (n+1) indices ending up on the same side
as the (n) indices.
That sounds like a bug. Can you give a specific example of an expression that
does that?
Aaron Meurer
Sure, I can attempt to give an example. I've been unable to reproduce
this behavior with a simple expression, so I am pasting below the full
expression that I've been using. Of course, it may be something that
I'm doing wrong, so I'm reluctant to cry wolf. But in a similar fashion
to working with a building, if the pipes are a bit leaky, perhaps a
small patch job would be beneficial if warranted.
The output of the example program below using Idx(n) shows that there
are terms with p[i,j,n+1] on both sides of the equality. The equality
== occurs three lines from the bottom of the output. See the first line
of the output for a p[1, 1, 1 + n] term. Here is the output:
WRONG:
p[-1, 1, n]/deltax**4 + p[1, 1, 1 + n]*(3/deltax**4 +
0.5/(deltax**2*deltay**2) - A/(deltat**2*deltax**2) + A*p[1, 1, -1 +
n]/(deltat**2*deltax**2*p[1, 1, n]) - A*p[0, 1,
n]/(deltat**2*deltax**2*p[1, 1, n]) - A*p[2, 1, 1 +
n]/(deltat**2*deltax**2*p[1, 1, n])) + p[2, 1, 1 + n]*(-4/deltax**4 -
1.0/(deltax**2*deltay**2) + 2*A/(deltat**2*deltax**2) - A*p[1, 1, -1 +
n]/(deltat**2*deltax**2*p[1, 1, n])) - 4*p[0, 1, n]/deltax**4 + 3*p[1,
1, n]/deltax**4 + 0.5*p[1, 1, n]/(deltax**2*deltay**2) + 0.5*p[1, -1,
n]/(deltax**2*deltay**2) + 0.5*p[-1, 1, n]/(deltax**2*deltay**2) +
0.5*p[-1, -1, n]/(deltax**2*deltay**2) + 2.0*p[0, 0,
n]/(deltax**2*deltay**2) - 1.0*p[0, 1, n]/(deltax**2*deltay**2) -
1.0*p[0, -1, n]/(deltax**2*deltay**2) - 1.0*p[1, 0,
n]/(deltax**2*deltay**2) - 1.0*p[-1, 0, n]/(deltax**2*deltay**2) -
2*A*p[1, 1, n]/(deltat**2*deltax**2) + A*p[1, 1, 1 +
n]**2/(deltat**2*deltax**2*p[1, 1, n]) - A*p[1, 1, -1 + n]*p[0, 1,
n]/(deltat**2*deltax**2*p[1, 1, n]) + 2*A*p[1, 1, n]*p[0, 1,
n]/(deltat**2*deltax**2*p[1, 1, n]) == (-1/deltax**4 -
0.5/(deltax**2*deltay**2))*p[3, 1, 1 + n] + p[1, 2, 1 +
n]/(deltax**2*deltay**2) + p[2, 3, 1 + n]/(deltax**2*deltay**2) + p[3,
2, 1 + n]/(deltax**2*deltay**2) - 0.5*p[1, 3, 1 +
n]/(deltax**2*deltay**2) - 0.5*p[3, 3, 1 + n]/(deltax**2*deltay**2) -
2.0*p[2, 2, 1 + n]/(deltax**2*deltay**2) - A*p[1, 1, -1 +
n]/(deltat**2*deltax**2)
If we use Idx(n+1) instead, the output appears to be right, and the
p[i,j,n+1] are only on one side of the expression. This seems to be good:
RIGHT:
(deltax**(-4) + 0.5/(deltax**2*deltay**2))*p[3, 1, 1 + n] + p[1, 1, 1 +
n]*(3/deltax**4 + 0.5/(deltax**2*deltay**2) - A/(deltat**2*deltax**2) +
A*p[1, 1, -1 + n]/(deltat**2*deltax**2*p[1, 1, n]) - A*p[0, 1,
n]/(deltat**2*deltax**2*p[1, 1, n]) - A*p[2, 1, 1 +
n]/(deltat**2*deltax**2*p[1, 1, n])) + p[2, 1, 1 + n]*(-4/deltax**4 -
1.0/(deltax**2*deltay**2) + 2*A/(deltat**2*deltax**2) - A*p[1, 1, -1 +
n]/(deltat**2*deltax**2*p[1, 1, n])) + 0.5*p[1, 3, 1 +
n]/(deltax**2*deltay**2) + 0.5*p[3, 3, 1 + n]/(deltax**2*deltay**2) +
2.0*p[2, 2, 1 + n]/(deltax**2*deltay**2) - 1.0*p[1, 2, 1 +
n]/(deltax**2*deltay**2) - 1.0*p[2, 3, 1 + n]/(deltax**2*deltay**2) -
1.0*p[3, 2, 1 + n]/(deltax**2*deltay**2) + A*p[1, 1, 1 +
n]**2/(deltat**2*deltax**2*p[1, 1, n]) == -p[-1, 1, n]/deltax**4 -
3*p[1, 1, n]/deltax**4 + 4*p[0, 1, n]/deltax**4 + p[0, 1,
n]/(deltax**2*deltay**2) + p[0, -1, n]/(deltax**2*deltay**2) + p[1, 0,
n]/(deltax**2*deltay**2) + p[-1, 0, n]/(deltax**2*deltay**2) - 0.5*p[1,
1, n]/(deltax**2*deltay**2) - 0.5*p[1, -1, n]/(deltax**2*deltay**2) -
0.5*p[-1, 1, n]/(deltax**2*deltay**2) - 0.5*p[-1, -1,
n]/(deltax**2*deltay**2) - 2.0*p[0, 0, n]/(deltax**2*deltay**2) - A*p[1,
1, -1 + n]/(deltat**2*deltax**2) + 2*A*p[1, 1, n]/(deltat**2*deltax**2)
+ A*p[1, 1, -1 + n]*p[0, 1, n]/(deltat**2*deltax**2*p[1, 1, n]) -
2*A*p[1, 1, n]*p[0, 1, n]/(deltat**2*deltax**2*p[1, 1, n])
Nicholas
#begin sample code
from sympy import *
from sympy.tensor import Indexed, Idx, IndexedBase
from sympy.matrices import *
p = IndexedBase('p')
var('deltax deltay delta deltat A B')
i,j,n = symbols('i j n', integer = True)
#
# function to generate the FDTD scheme
#
# All p[i,j,n+1] indices should be collected and on the LHS of the
expression
def generate_fdtd(expr, nsolve):
expr0 = expr.lhs - expr.rhs
expr1 = expr0.expand()
w1 = Wild('w1',integer=True)
w2 = Wild('w2',integer=True)
expr2 = collect(expr1, p[w1, w2, nsolve])
expr3 = expr2.as_independent(*filter(lambda t: t.args[-1] ==
Idx(n), expr2.atoms(Indexed) ) )
#Uncomment line for proper behavior
#expr3 = expr2.as_independent(*filter(lambda t: t.args[-1] ==
Idx(n+1), expr2.atoms(Indexed) ) )
left = expr3[0]
right = expr3[1]
left = -1 * expr3[0]
left = left.expand()
left = collect(left, p[w1,w2,nsolve])
right = expr3[1]
return Eq(right, left)
def main():
expr = Eq( 0.5*((-2*p[1, 2, 1 + n] - 2*p[2, 1, 1 + n] - 2*p[2, 3, 1
+ n] - 2*p[3, 2, 1 + n]
+ 4*p[2, 2, 1 + n] + p[1, 1, 1 + n] + p[1, 3, 1 + n] +
p[3, 1, 1 + n]
+ p[3, 3, 1 + n])/(deltax**2*deltay**2) + (-2*p[0, 1,
n] - 2*p[0, -1, n]
- 2*p[1, 0, n] - 2*p[-1, 0, n] + 4*p[0, 0, n] + p[1,
1, n] + p[1, -1, n]
+ p[-1, 1, n] + p[-1, -1, n])/(deltax**2*deltay**2)) +
(-4*p[0, 1, n]
- 4*p[2, 1, 1 + n] + 3*p[1, 1, n] + 3*p[1, 1, 1 + n] +
p[-1, 1, n] + p[3, 1, 1 + n])/deltax**4,
A*(-2*p[1, 1, n] + p[1, 1, 1 + n] + p[1, 1, -1 +
n])*(-p[1, 1, n] - p[1, 1, 1 + n] + p[0, 1, n]
+ p[2, 1, 1 + n])/(deltat**2*deltax**2*p[1, 1, n]) )
print generate_fdtd(expr,n+1)
if __name__ == "__main__":
main()
# end sample code
--
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.
