LinuxQuestions.org (/questions/)
-   Programming (http://www.linuxquestions.org/questions/programming-9/)
-   -   Python -numerical integration, not getting any output. (http://www.linuxquestions.org/questions/programming-9/python-numerical-integration-not-getting-any-output-868802/)

 C.L. 03-15-2011 08:17 PM

Python -numerical integration, not getting any output.

I'm translating a numerical integration program from Fortran to Python, and I thought that I had gotten rid of all the syntax problems. However, my program runs for a long time without any writing anything into the output file. The output file is created, but there just isn't anything in it. The program has been running for over thirty minutes now. Although it may just be slow, I am new to Python and I am worried that I may have made a mistake somewhere, especially with the output code.
There are some variables that are defined, but not really used. Those are for when the code gets expanded.

Here is the Python code.
Code:

#! /usr/bin/env python # This program solves the parallel coupled electronic dimer from numpy import * import math # System Parameters # np    = number of LRC circuit units in the chain # dll    = inductance # dmm    = mutual inductance # dkk    = proportionality coefficient. dmm = dkk*dll # dcc    = capacitance # omega0 = natural frequency of a unit. omega0 ** 2 = 1/(LC) # GG    = gain/loss. Sqrt(L/C)/R # mu    = rescaled mutual induction. dmm / dll # tmax  = maximum time evolution # dt    = integration time step # ntimes = nsnaps of the number of times we are calling the output # subroutine in order to write the output data # init  = initial excitation GG = float(raw_input('Enter the value of gain/loss.\n')) """np = int(raw_input('Enter the number of LRC circuits.\n')) print np """ np = 2 init = 1 ntimes = 1000 dll = 1.0 dcc = 1.0 dmm = 0.8 mu =0 #dmm / dll tmax = 250.0 dt = 0.01 m = 0 dt2 = dt/2.0 ifile = tmax/(ntimes*dt) ydot = zeros([2,2], dtype=float) #-------------------------------------------- # Initial Condition # Temporary workaround; only for np=2. Not yet general. # y1 and y2 are the first and second circuits. # The first and second elements of each array correspond to Y[i,1] and Y[i,2] in the Fortran code. #-------------------------------------------- y = array([[1,0],[0,0]]) savey = zeros([np,2], dtype=float) phi = zeros([np,2], dtype=float) #z = zeros([np,2], dtype=float) #zt = zeros([np,2], dtype=float) # Equations of Motion ----------------------- # Evaluation of the right hand side of the ODE. def eqmot(np, t, y, ydot, dmm, dll, dcc, GG, mu):     mu2= mu * mu     var1 = (( GG / sqrt(dll * dcc)) + (mu2 * GG / (dll*dcc)) ** 4) / (1-( mu2 / (dll*dcc)**3))     var2 = (1/(dll*dcc))/(1-(mu2 / (dll*dcc)**3))     var3 = (2*GG*mu/(dll*dcc))/(1-(mu2 / (dll*dcc))**3)     var4 = (mu * math.sqrt(dll*dcc) / (dll*dcc)**2)/(1- (mu2 / (dll*dcc)**3))     ydot[0,0]=y[0,1]     ydot[1,0]=y[1,1]     ydot[0,1]= (-var1)*y[0,1]-var2*y[0,0]-var3*y[1,1]+var4*y[1,0]     ydot[1,1]= var1*y[1,1]-var2*y[1,0]+var3*y[0,1]+var4*y[0,0]     return ydot[0,0], ydot[1,0], ydot[0,1], ydot[1][1] # Runge Kutta ------------------------------- def runge(t, n, z, zt, dt2, savey, phi, m):     m = m + 1     if m == 1:         return 1, t, n, z, zt, dt2, savey, phi, m     elif m == 2:         for j in range(1,3):             for i in range(1, n + 1):                 savey[i-1][j-1]=z[i-1][j-1]                 phi[i-1][j-1]=zt[i-1][j-1]                 z[i-1][j-1]=savey[i-1][j-1] + dt2 * zt[i-1][j-1]         t = t + dt2         return 1, t, n, z, zt, dt2, savey, phi, m     elif m == 3:         for j in range(1,3):             for i in range(1,n+1):                 phi[i-1][j-1]=phi[i-1][j-1] + 2.0 * zt[i-1][j-1]                 z[i-1][j-1]=savey[i-1][j-1] + dt2 * zt[i-1][j-1]         return 1, t, n, z, zt, dt2, savey, phi, m     elif m == 4:         for j in range(1,3):             for i in range(1,n+1):                 phi[i-1][j-1] = phi[i-1][j-1] + 2.0 * zt[i-1][j-1]                 z[i-1][j-1] = savey[i-1][j-1] + 2.0 * dt2 * zt[i-1][j-1]         t = t + dt2         return 1, t, n, z, zt, dt2, savey, phi, m     elif m == 5:         for j in range(1,3):             for i in range(1,n+1):                 z[i-1][j-1]=savey[i-1][j-1] + (phi[i-1][j-1] + zt[i-1][j-1])*dt2/3.0         m = 0         return 1, t, n, z, zt, dt2, savey, phi, m     else:         print 'M has gone out of its range.' # crfile ---------------------------------- def crfile(y,ydot,np,ntimes,t,tmax,intex,mu,dcc,dll,dmm):     energyL1 =0.5*dll*y[0][0]**2     energyL2 =0.5*dll*y[1][0]**2         energyC1 =0.5*(-dcc*dll*y[0][1]-dcc*dmm*y[1][1]**2/dcc)     energyC2 =0.5*(-dcc*dll*y[1][1]-dcc*dmm*y[0][1]**2/dcc)     energyM =dmm*y[0][0]*y[1][1]     energytot =energyL1+energyL2+energyC1+energyC2+energyM     return t, energyL1, energyL2, energyC1, energyC2, energytot   fout = open('output.txt', 'w') # Integration ------------------------------------ t = 0.0 for i in range(1,ntimes+1):     for j in range(1, int(ifile+1)):         k = 1         while k == 1:             ydot[0,0], ydot[1,0], ydot[0,1], ydot[1][1] = eqmot(np, t, y , ydot, dmm, dll, dcc, GG, mu)             k, t, np, y, ydot, dt2, savey, phi, m = runge(t, np, y , ydot, dt2, savey, phi, m)         crfile(y,tdot,np,ntimes,t,tmax,i,mu,dcc,dll,dmm)         fout.write(crfile)    f.out.close()
Here is the original Fortran code. There are some errors with variable/array names in the 'crfile' subroutine, but I don't think there are any syntax mistakes.
Code:

Program electro_dimer_evol     !-------------------------------------------------   ! This program solves the coupled electronic dimer   ! This is for the parallel coupled dimer   !   ! - M.C. Zheng, T. Kottos Feb 2011   !--------------------------------------------------   IMPLICIT NONE   INTEGER :: NP, i, ii, j, k, kk, ntimes, ifile, init   INTEGER :: runge, m, nsnaps   REAL(8) :: omega0, GG, mu, tmax   REAL(8) :: dll, dmm, dcc, dkk   REAL(8) :: dt, dt2, t   REAL(8), DIMENSION(:,:), ALLOCATABLE :: Y, YDOT, phi, savey   !-----------------SYSTEM PARAMETERS-----------------   !   ! NP    -- number of LRC circuits units in the chain   ! dll    -- inductance   ! dmm    -- mutual inductance   ! dkk    -- proportionality coefficient dmm = dkk*dll   ! dcc    -- capacitance   ! omega0 -- natural frequency of a unit: \omega0^2 = 1/LC   ! GG    -- gain/loss: gg=R*SQRT[C/L]  use new parameter   ! mu    -- rescaled mutaul induction \mu = MC = dkk/omega0^2, use new para   !--------------------------------------------------   ! Tmax  -- maximum time evolution   ! dt    -- integration time step   ! ntimes -- nsnaps of the number of times we are calling the output subroutine in order to write the output data   ! init  -- initial excitation   !------------------------------------------------------   Print*, 'Enter: GG'   Read(5,*) GG   Np = 2   init = 1   ntimes = 1000 !nn = Tmax/(dt*ntimes); nn INTEGER   dll = 1.d0   dcc = 1.d0 !  dkk = .5d0   dmm = .8d0  ! omega0 = 1.d0/dsqrt(dll*dcc)   mu = dmm/dll !  GG = 1.8d0   Tmax = 250.d0   dt = 0.01d0     !---------------------------------------------------------   ALLOCATE(Y(Np,2),YDot(Np,2),phi(Np,2),savey(Np,2))   OPEN(2, file = 'IE.dat')   m = 0   dt2 = dt/2.d0   ifile = tmax/(ntimes*dt)   !---------Initial Condition----------------------------------     CALL INCOND(Np, Y, init)   !------------LRC Integration-------------------------------               t=0.d0   Do i=1, ntimes     Do j=1, ifile         k = runge(t,NP,Y,YDOT,dt2,savey,phi,m)         Do while (k .EQ. 1)           call eqmot(NP,T,Y,YDOT,dmm,dll,dcc,GG,mu)           k = runge(t,NP,Y,YDOT,dt2,savey,phi,m)          END Do     END Do     call crfile(Y,YDOT,np,ntimes,t,tmax,i,mu,dcc,dll,dmm)   END Do     Close(2)     DEALLOCATE(Y,YDot,phi,savey)   END Program electro_dimer_evol !----------------------------------------------------------- subroutine incond(np,Y,init)         integer np,init,i   real*8 Y(np,2)       do i=1,np         Y(i,1)=0.0         Y(i,2)=0.0       enddo       Y(init,1)=1.d0       return       end ! ======================================================================       SUBROUTINE crfile(Y,YDOT,np,ntimes,t,tmax,intex,mu,dcc,dll,dmm)       implicit none       integer j,NP,ntimes,intex       real*8  Y(NP,2),YDOT(NP,2),tmax       real*8  t,omega0,GG,mu,dcc,dll,dmm       real*8  energyL1, energyL2, energyC1, energyC2, energyM       real*8  energytot       energyL1 = .5d0*dll*Y(1,1)**2       energyL2 = .5d0*dll*Y(2,1)**2             energyC1 = .5d0*(-dcc*dll*YDOT(1,1)-dcc*dmm*YDOT(2,1))**2/dcc       energyC1 = .5d0*(-dcc*dll*YDOT(2,1)-dcc*dmm*YDOT(1,2))**2/dcc       energyM = dmm*Y(1,1)*Y(2,2)       energytot = energyL1+energyL2+energyC1+energyC2+energyM             !write out energies       write(2,*)t,energyL1,energyL2,energyC1,energyC2,energytot       return       end ! ======================================================================       SUBROUTINE eqmot(NP,T,Y,YDOT,dmm,dll,dcc,GG,mu)  !c*********************************************************************** !c  Routine for evaluating the right hand side of the ODE. !c  The number of variables is 4*N. !c  The correspondences of the varibles are, !c  Y(NP,2)=(Q(np),I(np)); YDOT(NP,2)=(dQ(np)/dt,dI(np)/dt) !c  !c***********************************************************************       implicit none          integer j,NP       real*8  Y(NP,2),YDOT(NP,2)       REAL*8  dll, dmm, dcc       real*8  t,GG,mu,pi,pi2       real*8  mu2       real*8 var1,var2,var3,var4 !c******************* EQUATIONS OF MOTION *******************************       mu2 = mu*mu       var1= ((GG/sqrt(dll*dcc))+(mu2*GG/(dll*dcc)**4))/(1-(mu2/(dll*dcc)**3))       var2= (1/(dll*dcc))/(1-(mu2/(dll*dcc)**3))       var3= (2*GG*mu/(dll*dcc))/(1-(mu2/(dll*dcc)**3))       var4= (mu*sqrt(dll*dcc)/(dll*dcc)**2)/(1-(mu2/(dll*dcc)**3))       YDOT(1,1) = Y(1,2)       YDOT(2,1) = Y(2,2)       YDOT(1,2) = -var1*Y(1,2)-var2*Y(1,1)-var3*Y(2,2)+var4*Y(2,1)       YDOT(2,2) =  var1*Y(2,2)-var2*Y(2,1)+var3*Y(1,2)+var4*Y(1,1)       RETURN       END !==================================================================== function runge(t,N,z,zt,dt2,savey,phi,m)   integer :: runge, m, i, j   integer :: N   real*8 dt2,t,phi(N,2),savey(N,2)   real(8), dimension(N,2) :: z,zt     m=m+1   selectcase(m)   case(1)     runge=1     return   case(2)     do  j=1,2         do  i=1,N           savey(i,j)=z(i,j)           phi(i,j)=zt(i,j)           z(i,j)=savey(i,j)+dt2*zt(i,j)         enddo     enddo     t=t+dt2     runge=1     return   case(3)     do  j=1,2         do  i=1,N           phi(i,j)=phi(i,j)+2D0*zt(i,j)           z(i,j)=savey(i,j)+dt2*zt(i,j)         enddo     enddo     runge=1     return   case(4)     do j=1,2         do  i=1,N           phi(i,j)=phi(i,j)+2D0*zt(i,j)           z(i,j)=savey(i,j)+2D0*dt2*zt(i,j)         enddo     enddo     t=t+dt2     runge=1     return   case(5)     do j=1,2         do  i=1,N           z(i,j)=savey(i,j)+(phi(i,j)+zt(i,j))*dt2/3D0         enddo     enddo     m=0     runge=0     return   end select end function runge

 bgeddy 03-16-2011 01:22 AM

The usual way to approach these things is to insert print statements to output values from your loops so you can tell whats happening in your code. Python is very forgiving and just because the program may run doesn't mean you have no syntax errors! Ther errrors may be in a part of your code that is never being reached.For example here:
Code:

        crfile(y,tdot,np,ntimes,t,tmax,i,mu,dcc,dll,dmm) # returned results not being stored anywhere         fout.write(crfile)    # trying to write a function to a file f.out.close() # possibly meant to be fout.close()
If the file write was being reached it would cause a TypeError. These are just initial obervations without going into the code in any detail. Start be putting in print statements in your loops (for and while) to output variables and create counters to see where the code is up to. This will show where it's getting stuck in loops. It looks like you are getting stuck in the inner while loop. Fix these then other errors will become apparent. Also it makes things more readable if you organize your code so you have for example any module imports at the top (as you already have), then function definitions then the main variable definitions then the actual program code getting input and looping when the functions are called. Your code mixes main code, variable definitions and function definitions making it very hard to read. Firstly though put some output so you can start to debug.

 C.L. 03-16-2011 06:54 AM

Thank you. I shall try these things right away.

Edit: At least one of the problems seems to be that k does not change.

 All times are GMT -5. The time now is 09:44 PM.