Simulation of an one dimensional driven diffusive system I'm currently writing a simulation in python with scipy and matplotlib to reproduce an one dimensional driven diffusive system described in this paper from M.R. Evans et al.
The system consists of positive, negative and hole particles. On the left side of the system the positive (negative) particles are produced on the left (right) side with a possibility a and destroyed at the right (left) side with the possibility b. 
In some cases the system should show a flip between positive and negative high density states and I'm trying to reproduce this behaviour with my simulation. But all I can see is an increasing current within my simulation and I can't find any problems in my code explaining such behaviour.
Does anyone here have some experience simulation such or similar systems?
Cheers,
Florian
 A: I think your bug comes from your violating my iron-clad rule:
Never modify an array in-place (unless you already have correctly-working code and you are optimizing it for speed).
The reason is, as you're modifying it, the array is a mix of new stuff and old stuff, and it's easy to lose track and make mistakes. (Incidentally, this is one reason that it's easier to write bug-free code when you adopt a functional programming language style -- this programming style more-or-less forbids you from modifying arrays in-place and similar bad ideas.)
If you look at the for i in range(0,N) part of the code, you're checking l[i] vs l[i+1] when the former has already been modified during this iteration while the latter has not. So for example, ...100000... would get swapped to ...000001... in a single timestep - the 1 would march rightward every step of the for loop. And if you iterate i downward instead of upward you'll be doing a different algorithm. I haven't read the paper but I really doubt that you intended this.
In the functional style, you would write a function whose input is the lattice in timestep t and whose output is the lattice in timestep t+1, and where the function does not work by modifying the input but rather starts with an empty output array and enters numbers into the array according to the appropriate function of the input entries.
The problem of modifying an array in place interacts with the more basic problem that you don't have assertions or unit-tests. For example, write out a few small lattices with 5 or 10 entries, figure out with a pencil and paper what they should become at the next time-step, and then write code that runs your algorithm and compares the answers. I think you'll find that the tests fail with your code as written right now.
