# Storing and taking average of random number generated plots [closed]

I am working on a code which generates a number of random numbers which I have to use in a set of matrices. This set is then multiplied with some other matrix set in a particular order. And finally I have to plot one element of the resultant product matrix. But I have to generate random numbers for each run of the code and then after a numbers of runs I have to take the average. At this moment, I am running and plotting the data for each run of the code separately but this method is extremely time consuming as I have to do this 100 times! So I was wondering that is there a way that I can ask MATHEMATICA to save the resultant matrix element for each run on its own and later on when 100 realizations of the code are completed then I recall all final matrices to take the average of the desired element? Below I am pasting a part of code for a single realization of random number generation and required plot:

TA = {{1 - I*β/(ωeg - δ), -I*β/(ωeg - δ)}, {I*β/(ωeg - δ),  1 + I*β/(ωeg - δ)}};

{L10, L9, L8, L7, L6, L5, L4, L3, L2, L1} = Sort[RandomReal[10, 10], Greater]
(* {9.63622, 9.3111, 8.82698, 7.81967, 6.85424, 2.93255, \
2.56615, 2.46661, 1.51094, 0.0989493} *)

TF1 = {{Exp[4*π*L1*I*δ], 0}, {0,  Exp[-4*π*L1*I*δ]}};
TF2 = {{Exp[4*π*(L2 - L1)*I*δ], 0}, {0, Exp[-4*π*(L2 - L1)*I*δ]}};
TF3 = {{Exp[4*π*(L3 - L2)*I*δ], 0}, {0, Exp[-4*π*(L3 - L2)*I*δ]}};
TF4 = {{Exp[4*π*(L4 - L3)*I*δ], 0}, {0, Exp[-4*π*(L4 - L3)*I*δ]}};
TF5 = {{Exp[4*π*(L5 - L4)*I*δ], 0}, {0, Exp[-4*π*(L5 - L4)*I*δ]}};
TF6 = {{Exp[4*π*(L6 - L5)*I*δ], 0}, {0, Exp[-4*π*(L6 - L5)*I*δ]}};
TF7 = {{Exp[4*π*(L7 - L6)*I*δ], 0}, {0, Exp[-4*π*(L7 - L6)*I*δ]}};
TF8 = {{Exp[4*π*(L8 - L7)*I*δ], 0}, {0, Exp[-4*π*(L8 - L7)*I*δ]}};
TF9 = {{Exp[4*π*(L9 - L8)*I*δ], 0}, {0, Exp[-4*π*(L9 - L8)*I*δ]}};
TF10 = {{Exp[4*π*(L10 - L9)*I*δ], 0}, {0, Exp[-4*π*(L10 - L9)*I*δ]}};
{{m11, m12}, {m21, m22}} = TF1.TA.TF2.TA.TF3.TA.TF4.TA.TF5.TA.TF6.TA.TF7.TA.TF8.TA.TF9.TA.TF10.TA;

m22;
t10a = 1/m22 /. {β -> 0.16, ωeg -> 1};(*α=k*L, δ is the plotting frequency*)
T10a = Abs[t10a]^2;
R10a = 1 - T10a;

Plot[{R10a, T10a}, {δ, 0, 2},
PlotStyle -> {{Thickness[0.015], Red, Dashed}, {Thickness[0.015],
Black}}, Frame -> True, FrameLabel -> {"Δ", ""},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 28},
ImageSize -> {600, 500}, PlotRange -> All,
PlotLabel -> Style["", FontSize -> 28]]


I shall be extremely thankful for help and sorry if I am asking something very basic as I am a beginner in MATHEMATICA.

• You can shorten the code a lot by storing the 10 TF# and L# in a list and then compute the matrix m like so: L = Sort[RandomReal[10, 10]]; TF = {{Exp[4 π # I δ], 0}, {0, Exp[-4 π # I δ]}} & /@ ({First@L}~Join~Differences@L); m = Block[{TA}, Dot @@ Riffle[TF, TA, {2, 20, 2}]] – LLlAMnYP Nov 16 '15 at 7:15

As per my comment, this approach works:

TA={{1-I*β/(ωeg-δ),-I*β/(ωeg-δ)},{I*β/(ωeg-δ),1+I*β/(ωeg-δ)}};
Do[L = Sort[RandomReal[10, 10]];
TF = {{Exp[4 π # I δ], 0}, {0, Exp[-4 π # I δ]}} & /@ ({First@L}~Join~Differences@L);
m = Block[{TA}, Dot @@ Riffle[TF, TA, {2, 20, 2}]];
T[i] = With[{β = 0.16, ωeg = 1}, Evaluate@(Abs[m[[2, 2]]]^-2)];, {i, 1, 100}]


Plotting is ridiculously long, but a quick look at the results is possible like so:

avg[δ_] = Mean[Array[T, 100]];
{#, avg@#} & /@ Range[0., 2., 0.05];
ListPlot@%


A band structure can clearly be seen.

Update:
After this I started this code:

Do[pt[i] = {.001 i, avg@(.001 i)}; count = i;, {i, 0, 2000}]


and went for a coffee. When I came back, it was still running, so I aborted it and ran

count
(* 461 *)
ListPlot[Array[pt, 461]]


Then I let MMA carry on with calculations like so:

Do[pt[i] = {.001 i, avg@(.001 i)}; count = i;, {i, 461, 2000}]


It'll get there eventually.

Update2:
After some time it's finally done, and, turns out, not much of a band structure there:

ListPlot[Array[pt, 2000]]


Using what you have already coded, first put them inside a module (You can shorten this part as suggested by LLlAMnYP).

  f[] := Module[{TA, L10, L9, L8, L7, L6, L5, L4, L3, L2, L1, TF1, TF2,
TF3, TF4, TF5, TF6, TF7, TF8, TF9, TF10, m11, m12, m21, m22, t10a,
T10a, R10a},
TA = {{1 -
I*\[Beta]/(\[Omega]eg - \[Delta]), -I*\[Beta]/(\[Omega]eg - \
\[Delta])}, {I*\[Beta]/(\[Omega]eg - \[Delta]),
1 + I*\[Beta]/(\[Omega]eg - \[Delta])}};

{L10, L9, L8, L7, L6, L5, L4, L3, L2, L1} =
Sort[RandomReal[10, 10], Greater];

TF1 = {{Exp[4*\[Pi]*L1*I*\[Delta]], 0}, {0,
Exp[-4*\[Pi]*L1*I*\[Delta]]}};
TF2 = {{Exp[4*\[Pi]*(L2 - L1)*I*\[Delta]], 0}, {0,
Exp[-4*\[Pi]*(L2 - L1)*I*\[Delta]]}};
TF3 = {{Exp[4*\[Pi]*(L3 - L2)*I*\[Delta]], 0}, {0,
Exp[-4*\[Pi]*(L3 - L2)*I*\[Delta]]}};
TF4 = {{Exp[4*\[Pi]*(L4 - L3)*I*\[Delta]], 0}, {0,
Exp[-4*\[Pi]*(L4 - L3)*I*\[Delta]]}};
TF5 = {{Exp[4*\[Pi]*(L5 - L4)*I*\[Delta]], 0}, {0,
Exp[-4*\[Pi]*(L5 - L4)*I*\[Delta]]}};
TF6 = {{Exp[4*\[Pi]*(L6 - L5)*I*\[Delta]], 0}, {0,
Exp[-4*\[Pi]*(L6 - L5)*I*\[Delta]]}};
TF7 = {{Exp[4*\[Pi]*(L7 - L6)*I*\[Delta]], 0}, {0,
Exp[-4*\[Pi]*(L7 - L6)*I*\[Delta]]}};
TF8 = {{Exp[4*\[Pi]*(L8 - L7)*I*\[Delta]], 0}, {0,
Exp[-4*\[Pi]*(L8 - L7)*I*\[Delta]]}};
TF9 = {{Exp[4*\[Pi]*(L9 - L8)*I*\[Delta]], 0}, {0,
Exp[-4*\[Pi]*(L9 - L8)*I*\[Delta]]}};
TF10 = {{Exp[4*\[Pi]*(L10 - L9)*I*\[Delta]], 0}, {0,
Exp[-4*\[Pi]*(L10 - L9)*I*\[Delta]]}};
{{m11, m12}, {m21, m22}} =
TF1.TA.TF2.TA.TF3.TA.TF4.TA.TF5.TA.TF6.TA.TF7.TA.TF8.TA.TF9.TA.\
TF10.TA;

m22;
t10a = 1/m22 /. {\[Beta] -> 0.16, \[Omega]eg -> 1};(*\[Alpha]=k*
L,\[Delta] is the plotting frequency*)
T10a = Abs[t10a]^2;
R10a = 1 - T10a;
{R10a, T10a}]


The above function returns {R10a, T10a} each time it is executed for a random selection of all L#. Now we can run it n times and collect all the results. Find their mean and plot the results.

 plot[n_] :=
Module[{allMatrix, allR10a, allT10a, mean1, mean2},
allMatrix = Table[f[], {i, 1, n}];
allR10a = Table[allMatrix[[i, 1]], {i, 1, Length[allMatrix]}];
allT10a = Table[allMatrix[[i, 2]], {i, 1, Length[allMatrix]}];
mean1 = Mean[allR10a];
mean2 = Mean[allT10a];
Plot[{mean1, mean2}, {\[Delta], 0, 2},
PlotStyle -> {{Thickness[0.015], Red, Dashed}, {Thickness[0.015],
Black}}, Frame -> True, FrameLabel -> {"\[CapitalDelta]", ""},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 28},
ImageSize -> {600, 500}, PlotRange -> All,
PlotLabel -> Style["", FontSize -> 28]] ]


For n=100 it takes around 3 min on my system to generate the plot shown below.

 plot[100]


• Why do you use Module, but proceed to define global variables within the module. That is not good Mathematica coding style. – m_goldberg Nov 16 '15 at 16:38
• Are you referring to mean1 and mean2 being defined as local variables in the second module?. You are suggesting that since these two variables are the final result that are plotted, they ought to be global? Now I have removed them from the local variable list. – Hubble07 Nov 16 '15 at 17:13
• No. I am suggesting that all the variables to which assignments are made to in your modules be localized. – m_goldberg Nov 16 '15 at 17:16
• OK understood and edited accordingly. Thanks and henceforth I will follow this practice. – Hubble07 Nov 16 '15 at 17:43
• Many thanks LLlAMnYP, Hubble07 and m_goldberg (especially the first two) for answering my question. That was indeed what I was trying to do...And I have learned new things from your reply. Thanks again! For belisarius has settled, MarcoB, Lou, blochwave, gpap: I don't know your reason for placing a hold as I tried to be as explicit and as clear as possible in asking this problem! – Imran Nov 19 '15 at 5:04