Re: Compiling numerical iterations
- To: mathgroup at smc.vnet.net
- Subject: [mg129950] Re: Compiling numerical iterations
- From: Dana DeLouis <dana01 at icloud.com>
- Date: Wed, 27 Feb 2013 23:43:34 -0500 (EST)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- Delivered-to: l-mathgroup@wolfram.com
- Delivered-to: mathgroup-newout@smc.vnet.net
- Delivered-to: mathgroup-newsend@smc.vnet.net
Hi. I see you've changed the code a little.
Here are just a few suggestions going forward.
> time = AbsoluteTiming[
> tabwalk = Table[0, {i, 1, maxidx}];
> a = 0;
> For[i = 1, i <= maxidx, i++,
> tabwalk[[i]] = a;
> If[RandomReal[{-1, 1}] > 0, a = a + 1, a = a - 1];
> ];
> ];
As mentioned,
tabwalk = Table[0, {maxidx}];
All the above code might be written as:
SeedRandom[1234567890];
NestWhileList[#-Log[1-RandomReal[]]&,0,#<100&] ;
%== tabzeit
True
That gets rid of 1 of your 'For statements.
This statement also has a loop:
Table[{tabzeit[[i]],tabwalk[[i]]},{i,1,maxidx}]
It might be better written as
Transpose[{tablet,tab walk}]
For debugging, and monitoring your code, one technique is to add a label to your outputs.
Instead of
Print[tabwalk];
Print a label, as well as a 'New Line" for spacing.
Print["tabwalk \n\n", tabwalk];*)
Add a style if you want to be fancy.
Print[Style["Corr: \n\n",Bold,Blue], corr, "\n"]];
Your new Corr array is now different, and really hard to follow, so here is just a programming technique.
One technique is to intercept the data, and change it to a symbolic array.
Then continue with the corr function.
Then, look at the output to see if there are any simpler patterns.
That's one way to see that all your hard work in your first code was really just a simple vector dot product.
Here's what I mean.
tabwalk = { has all your data at this point, but we'll change it}
tabwalk=Array[tabwk,Length[tabwalk]];
(*continue with code *)
For[k=1,k<=zeitmax/10,k++,
For[n=1,n<=zeitmax/10*9,n++,
corr[[k]=corr[[k]]+tabwalk[[SucheIndex[tabzeit,n-1,maxidx]]]*tabwalk[[SucheIndex[tabzeit,n+k-2,maxidx]]]/zeitmax*10/9/numteilchen;
]
];
corr[2]
The first few terms are:
1/90 tabwk[1] tabwk[4]+tabwk[4]^2/45+1/90 tabwk[4] tabwk[5]+1/90 tabwk[5] tabwk[7]+1/90 tabwk[7] tabwk[10] .
Your complex answer does not have all the 2 vector terms, and drops the corresponding first terms.
It has tabwk[4], tabwk[5],tabwk[7],tabwk[10], tabwk[12] etc.
It seems to be only keeping the terms of your complex function:
Table[SucheIndex[tabzeit,j,maxidx],{j,20}]//Union
{4,5,7,10,12,13,14,15,16,17,20,22}
I really don't know how to simplify this complex equation, as I don't know what it's really doing.
Don't forget to checkout the new CorrelationFunction in Mathematica 9.
(I haven't figured out exactly what it does just yet)
There is also a Correlation function as well.
Good luck.
= = = = = = = = = =
HTH :>)
Dana DeLouis
Mac & Mathematica 9
= = = = = = = = = =
On Tuesday, February 26, 2013 1:09:32 AM UTC-5, firlefranz wrote:
> Thanks a lot! To be honest, some of the commands Ray is using, I've never seen before. I stopped using mathematica before version 5 came out.
>
>
>
> Coming back to Peters statement of exporting the code from Mathematica to C. How this can be done starting from my or Ray's code? There is an automated C-code-gernerator implemented in Mathematica 9, am I right?
>
>
>
>
>
> As a step forward to the first little peace of code, here is another try, which is not really optimized. The be concrete, I try to simulate the autocorrelation function of a random walk, which is doing a step at none equally distant time steps. This has to been done for a long random walk for many particles. Instead of doing an average over many random walks and calculate one autocorrelation function I want to simulate many correlation functions and make an average over them. Since the time steps are non equal, I wrote a sub-function, which creates a new time axis and taking the necessary value for the random walk from the first table.
>
>
>
> Here is what I come up with. It's running in a reasonable time for one particle, but for a real statistic ensemble, I have to do it over 1.000.000 particles for a long time. Optimizing this or (probably better) exporting it to C would hopefully help a lot. So how to export it?
>
>
>
> Clear["Global`*"]
>
> SeedRandom[1234567890];
>
> zeitmax = 100;(* muss ein Integer sein *)
>
> numteilchen = 1;
>
> tauj = 1;
>
> corr = Table[0, {i, 1, zeitmax/10}];
>
>
>
> SucheIndex[zeitliste_, zeit_, maxindex_] :=
>
> Module[{i},
>
> For[i = 1, i <= maxindex, i++,
>
> If[zeitliste[[i]] > zeit, Break[]];
>
> ];
>
> i - 1
>
> ];
>
>
>
> For[j = 1, j <= numteilchen, j++,
>
> (* Zeitachse generieren von 0 bis zeitmax *)
>
> t = 0;
>
> i = 1;
>
> tabzeit = {};
>
> time = AbsoluteTiming[While[True,
>
> tabzeit = Append[tabzeit, t];
>
> dt = -tauj*Log[1 - RandomReal[]];
>
> If[t > zeitmax, Break[]];
>
> t = t + dt;
>
> i++;
>
> ];
>
> ];
>
> Print[time];
>
> maxidx = i;
>
>
>
> (* Random Walk *)
>
> time = AbsoluteTiming[
>
> tabwalk = Table[0, {i, 1, maxidx}];
>
> a = 0;
>
> For[i = 1, i <= maxidx, i++,
>
> tabwalk[[i]] = a;
>
> If[RandomReal[{-1, 1}] > 0, a = a + 1, a = a - 1];
>
> ];
>
> ];
>
> Print[time];
>
> (*tabwalk=Table[Subscript[b, i],{i,1,maxidx}];*)
>
>
>
> (* Korrelationsfunktion berechnen *)
>
> time = AbsoluteTiming[
>
> For[k = 1, k <= zeitmax/10, k++,
>
> For[n = 1, n <= zeitmax/10*9, n++,
>
> corr[[k]] =
>
> corr[[k]] +
>
> tabwalk[[SucheIndex[tabzeit, n - 1, maxidx]]]*
>
> tabwalk[[SucheIndex[tabzeit, n + k - 2, maxidx]]]/
>
> zeitmax*10/9/numteilchen;
>
> (*Print[corr//N];*)
>
> ];
>
> ];
>
> ];
>
> Print[time];
>
> Print[corr // N];
>
> ];
>
> Table[{tabzeit[[i]], tabwalk[[i]]}, {i, 1, maxidx}]
>
> corr // N