A have a problem with the speed of the program written in Python. The program is "Simulation of ferromagnetic particles in a magnetic field", and more specifically in the magnetically inert liquid. The program works, but very slowly compared to the same program written in C++, but I have it written in Python, on a project to study.

Overall, the program code is based on loops, with lots of floating point calculations. There is a random number of particles (generated randomly position) which interact with each other under the influence of a magnetic field.

This is the initial position:

http://i.stack.imgur.com/T15Bb.jpg

And the final:

http://i.stack.imgur.com/0nU5D.jpg

The main loop (in SymMain.py, with k variable) iterations are time steps are calculated coordinates of the particles in them at the time and the forces acting on it (the attractive force and a small repulsive force). In order to speed up, I wanted to use parallel processing to calculate the number of iterations at the same time, but this is not possible, because the calculation in one iteration depend on the calculations in the previous one.

I did not know that Python is so slower than C + +. As an example, the simulation of 529 particles in one time step is calculated (on my computer):

C + + ~ 0.5s

Python ~ 50s

So, how to speed up the program? Please help.

This code is only calculate, not drawing particles as in the pictures. The number of simulate partcile is on SymConst.py, this is nrH*nrL.

SymMain.py

#coding:windows-1250

from os import system
from SymCalc import *
from SymParticle import *


if __name__ == '__main__':
    App = SymCalc()
    App.MainLoop()

SymParticle.py

#coding:windows-1250

from random import randint
from math import *
from SymConst import *
from SymParticle import *

class SymCalc(object):
    def __init__(self):
        # declaration lists containing the properties of the particles
        ParticleList = []
        ParticleListTemp = []
        t = 0.0

        # the initial values ​​of particle
        for x in range(0, nParticle):
            ParticleList.append(Particle(x+1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 8e-6))
            ParticleListTemp.append(Particle(x+1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 8e-6))

        # generating random coordinates X, Y 
        for x in range(0, nParticle):
            self.Rand(ParticleList[x], x)

        # time steps
        for k in range(0, k0):
            print "Time step = {0}".format(k)

            # calculation of forces
            for i in range(0, nParticle):
                for j in range(0, i-1):
                    self.ParticleCalculate(ParticleList[i], ParticleList[j], 1)
                for j in range(i+1, nParticle):
                    self.ParticleCalculate(ParticleList[i], ParticleList[j], 1)

            # display data
            if(k%Constant == 0):
                for i in range(0, nParticle):
                    self.Print(k, i, ParticleList[i], t)

            # changing the position of the particle
            for i in range(0, nParticle):
                self.ChangeOfPosition (ParticleList[i], dt)

                # reset forces
                for j in range(0, nParticle):
                    self.ParticleCalculate(ParticleList[i], ParticleList[j], 0)
                    self.ParticleCalculate(ParticleList[i], ParticleListTemp[j], 0)

            # next time step
            t += dt


    # random coordinates of the particles
    def Rand(self, part, lp):
        l = lp%nrL    # współrzędna X pola      
        part.x = (l-nrL/2)*(L0/nrL) + ((randint(0,32767)%100)-50)*(L0/nrL-2*part.r)/99    # X

        h = (lp+1-(lp)%nrL)/nrL    # Y      
        part.y = (h-nrH/2)*(H0/nrH) + ((randint(0,32767)%100)-50)*(H0/nrH-2*part.r)/99    # współrzędna Y cząstki


    # calculating function of the force acting on the particles
    def ParticleCalculate(self, part, part_tmp, p):
        # auxiliary variables
        # r_4, dx, dy, rr, sum, sx, sy, chi_eff, tmp, frepx, frepy, f0
        md = []

        if(p == 0):
            part.frx = 0
            part.fry = 0
            part.fx = 0
            part.fy = 0

        if(p == 1):
            # versor coordinates connecting the geometrical center of the particle
            dx = part.x - part_tmp.x
            dy = part.y - part_tmp.y

            # the distance between the geometric means of the particle
            rr = sqrt(dx*dx + dy*dy)

            if(rr < 0.85*(part.r + part_tmp.r)):
               print "ERROR: Invalid distance between the particles! Simulation aborted..."

            if(rr >= 10*part.r):
                # magnetic dipoles
                chi_eff = (3.*(MI_P - 1.))/((MI_P - 1) + 3.)

                md.append((4.*MI_0*pi*part.r*part.r*part.r*chi_eff*M_H0)/3.0)
                md.append((4.*MI_0*pi*part_tmp.r*part_tmp.r*part_tmp.r*chi_eff*M_H0)/3.0)

                tmp = pow(rr,7)

                # first member
                sum = (5.*(md[0]*part.nx*dx + md[0]*part.ny*dy) * (md[1]*part_tmp.nx*dx + md[1]*part_tmp.ny*dy)) / tmp
                sx = sum*dx
                sy = sum*dy
                tmp = tmp / (rr*rr)

                # second member
                sum = (md[0]*part.nx*md[1]*part_tmp.nx + md[0]*part.ny*md[1]*part_tmp.ny) / tmp
                sx -= sum*dx
                sy -= sum*dy

                # third member
                sx -= (md[0]*(md[1]*part_tmp.nx*dx + md[1]*part_tmp.ny*dy)*part.nx + md[1]*(md[0]*part.nx*dx + md[0]*part.ny*dy)*part_tmp.nx)/tmp
                sy -= (md[0]*(md[1]*part_tmp.nx*dx + md[1]*part_tmp.ny*dy)*part.ny + md[1]*(md[0]*part.nx*dx + md[0]*part.ny*dy)*part_tmp.ny)/tmp

                # finally
                tmp = (-3./(4*pi*MI_0))
                sx *= tmp
                sy *= tmp
                part.fx += sx
                part.fy += sy

                # short-range repulsive force
                tmp = pow(rr,15)
                r_4 = pow((part.r+part.r),4)

                f0 = 3.*fabs(md[0])*fabs(md[1])/(4.*pi*MI_0*r_4)
                frepx = pow((part.r+part.r),15)*dx/(tmp*rr)*f0
                frepy = pow((part.r+part.r),15)*dy/(tmp*rr)*f0

                part.frx += frepx;
                part.fry += frepy;


    # change the position of the particle
    def ChangeOfPosition(self, part, dt):
        part.ddx = 0
        part.ddy = 0

        part.ddx = 1/(6*pi*part.r*eta)*(part.fx+part.frx)*dt
        part.ddy = 1/(6*pi*part.r*eta)*(part.fy+part.fry)*dt

        # particles new position value
        part.x += part.ddx
        part.y += part.ddy

        if(part.x < -L0/2):
            part.x += L0
        elif(part.x > L0/2):
            part.x -= L0
        elif(part.y < -H0/2):
            part.y += H0
        elif(part.y > H0/2):
            part.y -= H0



    # display data
    def Print(self, k, i, part, t):
        print "-"*50
        print "\nParticle {0}".format(i+1)
        print "The resultant magnetostatic force fx = {0}".format(part.fx - part.frx)
        print "The resultant magnetostatic force fy = {0}\n".format(part.fy - part.fry)
        if(i == nParticle-1):
            print "\n\t\t---t={0}[s]---".format(t)
        print "-"*50

SymParticle.py

#coding:windows-1250

class Particle(object):
    # generating a particle properties
    def __init__(self, num, x, y, fx, fy, frx, fry, nx, ny, ddx, ddy, r):
        self.num = num     
        self.x = x     
        self.y = y     
        self.fx = fx     
        self.fy = fy     
        self.frx = frx     
        self.fry = fry     
        self.nx = nx     
        self.ny = ny     
        self.ddx = ddx     
        self.ddy = ddy     
        self.r = r     

SymConst.py

#coding:windows-1250

### Constant
M_H0 = 3e4
MI_0 = 12.56e-7     
MI_P = 2000
eta = 0.1           
k0 = 10001
dt = 1e-6           # time step      
H0 = 5.95e-4
L0 = 5.95e-4
nrH = 4
nrL = 4
Constant = 100
nParticle = nrH*nrL    # number of particle