# Reducing the computation time used by parallel table when finding a fixed point

I am using Intel Core(™) processor, 32 bit operating system, 4GB RAM, Windows 7. I am running the following code and it is taking hours to obtain the result. Can any one help me in improving the codes so that it takes less time.

μ = 0.01214398779977472

μ1 = 328823.8592475095

Subscript[e, 1] = 0.0167

A = 389.1779396462019

ϕ = 0

Subscript[e, 2] = 0.0549

Vx =
x + μ/(y^2 + (-1 + x + μ)^2)^(3/2) - x μ/(y^2 + (-1 + x + μ)^2)^(3/2) -
μ^2/(y^2 + (-1 + x + μ)^2)^(3/2) - x/(y^2 + (x + μ)^2)^(3/2) -
μ/(y^2 + (x + μ)^2)^(3/2) + x μ/(y^2 + (x + μ)^2)^(3/2) +
μ^2/(y^2 + (x + μ)^2)^(3/2) + x μ1/(2 A^3) + (3 x μ1 Cos[2 ϕ])/A^3 +
(3 y μ1 Sin[2 ϕ])/A^3 - (μ1 Cos[ϕ] Subscript[e, 1]^2)/(2 A^2);

Vy =
y - (y μ)/(y^2 + (-1 + x + μ)^2)^(3/2) - y/(y^2 + (x + μ)^2)^(3/2) +
(y*μ)/(y^2 + (x + μ)^2)^(3/2) + (y*μ1)/(2 A^3) - (3 y μ1 Cos[2 ϕ])/A^3 +
(μ1 Sin[ϕ])/A^2 + (3 x μ1 Sin[2 ϕ])/A^3 - (μ1 Sin[ϕ] Subscript[e, 1]^2)/(2 A^2);

Vxx = D[Vx, x];
Vyy = D[Vy, y];
Vxy = D[Vx, y];
Vyx = D[Vy, x];

newton[{x_, y_}] :=
{x, y} +
{Simplify[-(Vx Vyy - Vy Vxy)/(Vyy Vxx - Vxy^2)],
Simplify[(Vx Vyx - Vy Vxx)/(Vyy Vxx - Vxy^2)]};

sol =
{{ 0.188277, -0.9824026}, { 0.188277, 0.9824026}, { 1.152706, 0},
{-0.998550, 0}, { 0.987856, 0}, { 0.835518, 0}, { -0.012144, 0}};

tab =
ParallelTable[
N[FixedPoint[newton, {i, j}, 10]],
{j, -2, 2, 0.003},
{i, -2, 2, 0.003}]

rules = Rule @@@ Transpose[{sol, Range[Length[sol]]}];

newtab = Map[First @ Nearest[rules, #] &, tab, {2}];

P1 =
ArrayPlot[newtab,
ColorFunction -> "Rainbow",  DataReversed -> True]

• Welcome to Mathematica.StackExchange. Please consult mathematica.stackexchange.com/help/how-to-ask on how to improve your question. As a first remark: Please post code that can be copied so that other user can 1.) really read it and 2.) experiment with it without needing to retype everything. – Henrik Schumacher Dec 16 '17 at 11:05
• I'm voting to close this question as off-topic because, as others already have stated, an image of the screen is not an acceptable question. Please edit the question to replace the screen-shot by the code itself in Mathematica format. – bbgodfrey Dec 16 '17 at 14:37
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey Dec 17 '17 at 3:31

The principle issue with the code in the question is that the definition of newton does not replace x and y in Vx, etc. by the values of the arguments of newton. There are a number of ways to resolve this issue, including defining

newton[{x1_, y1_}] := {x1, y1} + {Simplify[-(Vx Vyy - Vy Vxy)/(Vyy Vxx - Vxy^2)],
Simplify[(Vx Vyx - Vy Vxx)/(Vyy Vxx - Vxy^2)]}/.{{x -> x1, y -> y1}}


With this resolved, we can focus on increasing the computational speed. Since almost all the time goes into repeatedly evaluating newton, our effort should go there. A factor of three or more is obtained by redefining newton as

newton[{x1_, y1_}] := Block[{}, {vx, vy, vxx, vyy, vxy} = {Vx, Vy, Vxx, Vyy, Vxy} /.
{x -> x1, y -> y1}; {x1, y1} + {-(vx vyy - vy vxy), (vx vxy - vy vxx)}/(vyy vxx - vxy^2)]


which takes advantage of the fact that Vxy == Vyx and also eliminates duplicated evaluations of intermediate results. Eliminating Simplify also saves a bit of time, because it is ineffective in the present calculation. (N in the definition of tab also can be eliminated, because it has no effect, but doing so saves negligible time.) With these changes, the computation completes in about eleven minutes on my four-processor PC, yielding

Perhaps, restructuring Vx, Vy, and their derivatives could save a bit more time.

There are several issues:

1. newton does not replace x and y. (See bbgrodfrey's answer.)

2. Every call to newton calls also a Simplify which is pretty slow. Simplifying it once at the time of definition should suffice.

You can resolve point 1 and 2 by defining:

 newton[{x_, y_}] :=
Evaluate[{x, y} + {Simplify[-(Vx Vyy - Vy Vxy)/(Vyy Vxx - Vxy^2)],
Simplify[(Vx Vyx - Vy Vxx)/(Vyy Vxx - Vxy^2)]}];


Even that is probably still too slow. Compiling the code helps and it makes parallelization easier.

With[{xnew = N@newton[{x, y}][[1]], ynew = N@newton[{x, y}][[2]]},
newtonize = Compile[{{x0, _Real, 1}, {maxiter, _Integer}},
Block[{x, y, xbuffer},
x = CompileGetElement[x0, 1];
y = CompileGetElement[x0, 2];
Do[
xbuffer = xnew;
y = ynew;
x = xbuffer;
,
{maxiter}];
{x, y}],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"];
];


This CompiledFunction performs maxiter iterations on each input vector. Note that we use the option RuntimeAttributes -> {Listable} which creates a CompiledFunction that threads automatically over the matrix x0data of starting vectors below:

sol = DeveloperToPackedArray[N[{{0.188277, -0.9824026}, {0.188277, 0.9824026}, {1.152706, 0.}, {-0.998550, 0.}, {0.987856, 0.}, {0.835518, 0.}, {-0.012144, 0.}}]];
cols = DeveloperToPackedArray[
List @@@ ColorData["Rainbow"] /@ Rescale[Range[Length[sol]]]];
colorize = Compile[{{a, _Integer}, {cols, _Real, 2}},
CompileGetElement[cols, a],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];

δx = 0.003;
maxiter = 10;
tab = newtonize[
Outer[List, Table[N[j], {j, -2, 2, δx}], Table[N[j], {j, -2, 2, δx}]],
maxiter
]; // AbsoluteTiming // First
newtab = ArrayReshape[
Nearest[sol -> Automatic, Flatten[tab, 1]],
Dimensions[tab][[1 ;; 2]]
]; // AbsoluteTiming // First
size = 600;

P1 = ImageRotate[
ImageResize[Image[colorize[newtab, cols]], {size, size}]
]; // AbsoluteTiming // First
P1


0.578838

0.048854

0.183145

0.1231

Performing that with δx = 0.0005 and maxiter = 30 leads after 2 minutes to this picture:

Note also the following tweaks:

• We do no longer use a nested table for the construction of the initial points but Outer (among other things, this ensures that the generated array is a PackedArray).

• We do not Map a NearestFunction onto tab but rather apply Nearest with a syntax that enables internal low level parallelism.

• We ensure that sol is also a PackedArray so that Nearest can profit from this.

• We use Image in conjunction with a hand written method colorize instead of ArrayPlot in order to improve performance and stability (my FrontEnd had the urgent need to shut down several times).

Update

Corrected a bug in the CompiledFunction: The former version changed the value of x before the new value of y was computed. This was somewhat hidden by the way the code for x and y is injected into the compiled function. Buffering the new value for x (in the variable xbuffer) remedies this.

• Nice solution. FYI, the first issue you identified was introduced in edit 4, when someone was attempting to improve the format. By good fortune, I worked from edit 3, which did not contain this issue. – bbgodfrey Dec 17 '17 at 14:35
• @bbgodfrey I see. So I added edit 5 ;o) – Henrik Schumacher Dec 17 '17 at 14:38
• To understand why our final plots are quite different, I evaluated our respective expressions for newton[{-2., -2.}], both of which yield {-0.238379, -0.248411}. However, my FixedPoint[newton, {-2., 2.}, 10] yields {-0.998511, -0.00312046}, and your newtonize[{-2., -2.}, 10] yields {-0.156382, -1.00541}. To what do you attribute the discrepancy? By the way, where can I find the documentation for "CompileGetElement"? Thanks. – bbgodfrey Dec 17 '17 at 16:05
• @bbgodfrey Thanks for the counterexample! That really helped to debug the code (before that I assumed that this was related to some numerical errors). I corrected the bug and now the resulting picture coincides with yours now. – Henrik Schumacher Dec 17 '17 at 16:36
• @bbgodfrey I am afraid there is no documentation for CompileGetElement =o/ It pops up in the list returned by Sort@CompileCompilerFunctions[] (see also this post: mathematica.stackexchange.com/questions/1096/…) and that raised my curiosity. In some cases, Compile uses it as subsitute for Part in order to speed up code execution. When you compare the actual C-code that is generated then you see that CompileGetElement skips checks for the array bounds when looking up indices. So, this is very fast but also a bit risky to use... – Henrik Schumacher Dec 17 '17 at 16:48