Wednesday, August 1, 2007

Progress

The midterm evaluation is over and noone has failed, which is great. We started to work more closely on the #sympy at freenode IRC channel and SymPy has made a tremendous progress in past couple of weeks.

Our team was joined by Fredrik and Pearu, so now there are 8 regular contributors.

Fredrik wrote fast floating point arithmetic (in pure Python), faster than the python Decimal module. For example 10 000 decimal digits of Pi can be computed in 0.3s (and one million in less than 20 minutes) with SymPy - just run the pidigits.py in the examples directory. He among other things wrote a fast algorithm for calculating the gamma function, that in some cases is faster than in Mathematica.

We moved to the new core written by Pearu, that has made everything 10x faster or more. Basically, Pearu very cleverly designed everything so that no operation is made in vain. For example my little hobby relativity.py for calculating the Schwarzschild solution now calculates it in 1s (compared to more than 10s before)!


ondra@pc232:~/sympy/examples$ time python relativity.py
----------------------------------------
Christoffel symbols:
(1/2)*(FD(1,)(nu))(r)
(1/2)*(FD(1,)(nu))(r)

(1/2)*exp(-lam(r) + nu(r))*(FD(1,)(nu))(r)
(1/2)*(FD(1,)(lam))(r)
-r/exp(lam(r))
-r*sin(\theta)**2/exp(lam(r))

1/r
1/r
-sin(\theta)*cos(\theta)

1/sin(\theta)*cos(\theta)
1/sin(\theta)*cos(\theta)
1/r
1/r
----------------------------------------
Ricci tensor:
(1/2)*exp(-lam(r) + nu(r))*(FD(1, 1)(nu))(r) - 1/4*exp(-lam(r) + nu(r))*(FD(1,)(nu))(r)**2 + 1/r*exp(-lam(r) + nu(r))*(FD(1,)(nu))(r) + (1/2)*exp(-lam(r) + nu(r))*(-(FD(1,)(lam))(r) + (FD(1,)(nu))(r))*(FD(1,)(nu))(r) + (1/4)*exp(-lam(r) + nu(r))*(FD(1,)(lam))(r)*(FD(1,)(nu))(r)
-1/2*(FD(1, 1)(nu))(r) - 1/4*(FD(1,)(nu))(r)**2 + 1/r*(FD(1,)(lam))(r) + (1/4)*(FD(1,)(nu))(r)*(FD(1,)(lam))(r)
1 - 1/exp(lam(r)) + (1/2)*r/exp(lam(r))*(FD(1,)(lam))(r) - 1/2*r/exp(lam(r))*(FD(1,)(nu))(r)
sin(\theta)**2 - sin(\theta)**2/exp(lam(r)) + (1/2)*r*sin(\theta)**2/exp(lam(r))*(FD(1,)(lam))(r) - 1/2*r*sin(\theta)**2/exp(lam(r))*(FD(1,)(nu))(r)
----------------------------------------
solve the Einstein's equations:
lam(r) = -log(C1 + C2/r)
metric:
-C1 - C2/r 0 0 0
0 1/(C1 + C2/r) 0 0
0 0 r**2 0
0 0 0 r**2*sin(\theta)**2

real 0m1.092s
user 0m1.024s
sys 0m0.068s


Now we are in the phase of polishing things so that everything that used to work before works in the new core the same way (or better). Things that remain are printing and some limits are failing, everything else was already fixed.

One conclusion that can be made is, that Python on today computers is already fast enough to do something useful. Concerning the Pi, it's amusing to read some articles about people trying to calculate many digits, for example in one article, it took them a day to calculate one million digits on Pentium 90 (SymPy can do that in 20 minutes).

Also I found some page, that I cannot rediscover at the moment, that stated it took like 15 minutes to calculate the Schwarzschild with Maple on some older computer, so I am glad SymPy can do that in 1s (and actually the code is very slow, I am doing some calculations again and again, but I didn't find time yet to optimize it).

Anyway, the overall progress on SymPy looks very promising, I have no doubts that it will become a useful symbolic python library.

2 comments:

Naresh said...

Hey I'm looking for sympy code to compute christoffel symbols. Facing a lot of problems currently while using sympy. I hope you can help.

If you can kindly contact me at nareshshah139@gmail.com

We can also talk on skype if necessary. Skype id: nareshshah139

Naresh said...

Here's a copy of what I've written till now:


from sympy import *
def GLA(i,j,a):

theta = Symbol('theta')
chy = Symbol('chy')
sy = Symbol('sy')
t = Symbol('t')
temp1 = 0
temp2 = 0
temp3 = 0
sumTialpha = 0
sumTjalpha = 0

if a == 1:
xa = 'theta'
elif a == 2:
xa = 'sy'
elif a == 3:
xa = 'chy'
elif a == 4:
xa = 't'

csc = exp(-log(sin(theta)))

gl = [ [-exp(4*t), 0, 0, 0],[0, -exp(4*t)*sin(theta)**2-exp(6*t)*cos(theta)**2, -exp(6*t)*cos(theta), 0],[0,-exp(6*t)*cos(theta),-exp(6*t),0],[0, 0, 0, 1]]

#gl = [ [-exp(4*t), 0, 0, 0],[0, -exp(4*t) -exp(6*t)*(theta**2), -exp(6*t)*(theta), 0],[0,-exp(6*t)*(theta),-exp(6*t),0],[0, 0, 0, 1]]

gL = [[-1, 0, 0 , 0],[0,-1*(1/(sin(theta))*sin(theta)),cos(theta)/(sin(theta)*sin(theta)),0],[0,cos(theta)/(sin(theta)*sin(theta)),-1*(1/(sin(theta)*sin(theta))),0],[0,0,0,1]]

#gL = [[-1/exp(2*t), 0, 0 , 0],[0,-1/exp(2*t),theta/exp(2*t),0],[0,theta/exp(2*t),-(theta**2)/exp(2*t) - 1/exp(3*t),0],[0,0,0,1]]

T = [[[0,0,0,0],[0,0,sin(theta)*0.5,0],[0,sin(theta)*0.5,0,0],[0,0,0,0]],[[0,0.5*cot(theta),-0.5*csc,0],[0.5*cot(theta),0,0,0],[-0.5*csc,0,0,0],[0,0,0,0]],[[0,-0.5*csc,0.5*cot(theta),0],[-0.5*csc,0,0,0],[0.5*cot(theta),0,0,0]],[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]]

#T = [[[0,0,0,0],[0,-theta,-0.5,0],[0,-0.5,0,0],[0,0,0,0]],[[0,theta*0.5,0.5,0],[theta*0.5,0,0,0],[0.5,0,0,0],[0,0,0,0]],[[0,(1-theta**2)*0.5,-theta*0.5,0],[(1-theta**2)*0.5,0,0,0],[-theta*0.5,0,0,0],[0,0,0,0]],[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]]

k = diff(gl[i-1][j-1],xa)

while temp1 < 4:
sumTialpha = sumTialpha + T[temp1][i-1][a-1]*gl[temp1][i-1]
temp1 = temp1 + 1

while temp2 < 4:
sumTjalpha = sumTjalpha + T[temp2][j-1][a-1]*gl[temp2][j-1]
temp2 = temp2 + 1

return simplify(k - sumTialpha - sumTjalpha,ratio = 1.7)

def glab(i,j,a,b):

theta = Symbol('theta')
chy = Symbol('chy')
sy = Symbol('sy')
t = Symbol('t')
temp1 = 0
temp2 = 0
temp3 = 0
sumTibeta = 0
sumTjbeta = 0

if a == 1:
xa = 'theta'
elif a == 2:
xa = 'sy'
elif a == 3:
xa = 'chy'
elif a == 4:
xa = 't'

#gl = [ [-1, 0, 0, 0],[0, -1, -cos(theta), 0],[0,-cos(theta),-1,0],[0, 0, 0, 1]]

gl = [ [-exp(2*t), 0, 0, 0],[0, -exp(2*t) -exp(3*t)*(theta**2), -exp(3*t)*(theta), 0],[0,-exp(3*t)*(theta),-exp(3*t),0],[0, 0, 0, 1]]

#gL = [[-1, 0, 0 , 0],[0,-1*(1/(sin(theta))*sin(theta)),cos(theta)/(sin(theta)*sin(theta)),0],[0,cos(theta)/(sin(theta)*sin(theta)),-1*(1/(sin(theta)*sin(theta))),0],[0,0,0,1]]

gL = [[-1/exp(2*t), 0, 0 , 0],[0,-1/exp(2*t),theta/exp(2*t),0],[0,theta/exp(2*t),-(theta**2)/exp(2*t) - 1/exp(3*t),0],[0,0,0,1]]

#T = [[[0,0,0],[0,0,sin(theta)/2],[0,0,0]],[[0,cos(theta)/(2*sin(theta)),-1/sin(theta)],[0,0,0],[0,0,0]],[[0,cos(theta)/(2*sin(theta)),-1/sin(theta)],[0,0,0],[0,0,0]]]

T = [[[0,0,0],[0,-theta,-1/2],[0,-1/2,0]],[[0,theta/2,1/2],[theta/2,0,0],[1/2,0,0]],[[0,1/2 - (theta**2)/2,-(theta/2)],[1/2 - (theta**2)/2,0,0],[-theta/2,0,0]],[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]]

gla = GLA(i,j,a)

k = diff(gla,xa)

while temp1 < 4:
sumTibeta = sumTibeta + T[i-1][b-1][temp1-1]
temp1 = temp1+1

while temp2 < 4:
sumTjbeta = sumTjbeta + T[j-1][b-1][temp2-1]
temp2 = temp2+1

return k - sumTibeta*gla - sumTjbeta*gla