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

*To*: mathgroup at smc.vnet.net*Subject*: [mg115862] Re: Do I need MathLink to run finite-difference fast enough for*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Sun, 23 Jan 2011 05:35:37 -0500 (EST)

----- Original Message ----- > From: "James" <icorone at hotmail.com> > To: mathgroup at smc.vnet.net > Sent: Saturday, January 22, 2011 2:24:08 AM > Subject: [mg115846] Re: Do I need MathLink to run finite-difference fast enough for > 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]]] Slightly faster, because you sorta left 'a' and 'b' outside: cf = Compile[{{uIn, _Real, 2}, {vIn, _Real, 2}, {a, _Real}, {b, _Real}, {du, _Real}, {dv, _Real}, {dt, \ _Real}, {iterationsIn, _Integer}}, Module[ {u = uIn, v = vIn, lap, k = b + 1, kern = N[{{0, 1, 0}, {1, -4, 1}, {0, 1, 0}}]}, 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]]] Two questions: Do you know why this balks if I raise dt to, say, 0.02? (There may be a sound reason based on PDE stability of this numeric scheme, I'm not sure.) I raise this question because clearly you might be able to get better speed via larger time steps, if the numerics will support that. (2) Since u is updated and then used in finding the new v, I tried reversing these to see if it would change. It did. Is this expected behavior? Daniel Lichtblau Wolfram Research