MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Do I need MathLink to run finite-difference fast enough for

  • To: mathgroup at smc.vnet.net
  • Subject: [mg115846] Re: Do I need MathLink to run finite-difference fast enough for
  • From: James <icorone at hotmail.com>
  • Date: Sat, 22 Jan 2011 03:24:08 -0500 (EST)

Hello Mike,

This is a Brusselator model of chemical reaction diffussion that is employing the Euler numerical method to solve the PDE system:

u_t=D_u lap(u)+a-(b+1)u+u^2v

v_t=D_v lap(v)+bu-u^2v

For suitable values of the parameters, the characteristic Turing patterns of dots and stripes should emerge. 

Using Oliver's Compiler suggestions above under V7, I can run the simulator on a 64X64 grid 10000 times on my machine in about 5.5 seconds.  However, I feel that's still too slow to create a reasonable interactive Demonstration of the Brusselator where the user will be changing the paramters.  Ideally, I'd like to get it down to about 2 seconds.  I'm not sure how fast it would run in V8.  Here is the complete code using Oliver's suggestions with some additional improvements for stripes that runs in 5.5 seconds on my machine.  Can anyone suggest of a way to get it down to two seconds?

n = 64; 
a = 4.5; 
b = 7.5; 
du = 2; 
dv = 16; 
dt = 0.01; 
totaliter = 10000; 
u = a + 0.3*RandomReal[{-0.5, 0.5}, {n, n}]; 
v = b/a + 0.3*RandomReal[{-0.5, 0.5}, {n, n}]; 
cf = Compile[{{uIn, _Real, 2}, {vIn, _Real, 2}, 
     {aIn, _Real}, {bIn, _Real}, {duIn, _Real}, 
     {dvIn, _Real}, {dtIn, _Real}, {iterationsIn, 
      _Integer}}, Block[{u = uIn, v = vIn, lap, dt = dtIn, 
      k = bIn + 1, kern = N[{{0, 1, 0}, {1, -4, 1}, 
         {0, 1, 0}}], du = duIn, dv = dvIn}, 
     Do[lap = RotateLeft[u, {1, 0}] + RotateLeft[u, 
           {0, 1}] + RotateRight[u, {1, 0}] + 
          RotateRight[u, {0, 1}] - 4*u; 
        u = u + dt*(du*lap + a - u*(k - v*u)); 
        lap = RotateLeft[v, {1, 0}] + RotateLeft[v, 
           {0, 1}] + RotateRight[v, {1, 0}] + 
          RotateRight[v, {0, 1}] - 4*v; 
        v = v + dt*(dv*lap + u*(b - v*u)); , 
       {iterationsIn}]; {u, v}]]; 
Timing[c1 = cf[u, v, a, b, du, dv, dt, totaliter]; ]
ListDensityPlot[c1[[1]]]


  • Prev by Date: Re: Help with non-linear system of equations
  • Next by Date: Re: problem with LessEqual::nord:
  • Previous by thread: Re: Converting XML DATEEVENT to Mathematica AbsoluteTime
  • Next by thread: Re: Do I need MathLink to run finite-difference fast enough for