# Very slow mathematica finite differences

I have a simple code which solves an equation by an explicit method (FTCS). It takes mathematica several minutes (mathematica 10.0.2) to finish the calculation while the same code in Fortran runs less than a second (on the same machine, no parallelization). I expected Mathematica to be slower, but that slower? I think something weird is going on. Can somebody help?

ClearAll["Global*"];
(* define viscosity function *)
fGetNu := Compile[{r},
1.34*10^14*(r/AU2cm)^(3/2)
];
(* define conversion functions *)
yr2sec = Compile[{x},
3.155*10^7*x
];
AU2cm = Compile[{x},
1.496*10^13*x
];

calc = Compile[{},
(* define global constants *)
ngrid = 102;
rin = AU2cm[0.1];
rout = AU2cm;
tmax = yr2sec[2*10^6];
tout = yr2sec[1*10^5];
tcounter = 0;
c0 = 0.5;
(* arrays *)
x = ConstantArray[0, ngrid];
v = ConstantArray[0, ngrid];
r = ConstantArray[0, ngrid];
Σ = ConstantArray[0, ngrid];
vNew = ConstantArray[0, ngrid];
ν = ConstantArray[0, ngrid];
(* main part *)
(* setup the grid *)
dX = 2.0*(Sqrt[rout] - Sqrt[rin])/(ngrid - 2);
For[i = 1, i <= ngrid, i++,
x[[i]] = 2 Sqrt[rin] + (i - 1.5)*dX;
r[[i]] = 0.25*x[[i]]^2;
];
(* Define zero-torque (Sigma=0) boundary conditions *)
Σ[] = 0; Σ[[ngrid]] = 0;
(* Define initial conditions *)
For[i = 2, i <= ngrid - 1, i++,
If[TrueQ[r[[i]] <= AU2cm],
Σ[[i]] = 1000*(r[[i]]/AU2cm)^(-3/2);,
Σ[[i]] = 1*10^-4;
];
ν[[i]] = fGetNu[r[[i]]];
v[[i]] = 3/2 ν[[i]]*Σ[[i]]*Sqrt[r[[i]]];
];
(* initial time *)
t = 0; tnext = 0;
(* main loop *)
While[TrueQ[t <= tmax],
dtmin = 1*10^20;
For[i = 2, i <= ngrid - 1, i++,
(* this is what is in Armitage code, but NB.
he uses explicit FTCS and this breaks von Neumann condition for \
stability, which has strict inequality,
also the further from the limit, the better *)
dtzone = c0*(dX*x[[i]])^2/(24*ν[[i]]);
dtmin = Min[dtmin, dtzone];
];
dt = dtmin;
(* Evolve one timestep,store result in array Vnew *)
(* This is FTCS time step *)
For[i = 2, i <= ngrid - 1, i++,
vNew[[i]] =
v[[i]] +
12*ν[[i]]*
dt (v[[i - 1]] - 2*v[[i]] + v[[i + 1]])/(x[[i]]*dX)^2;
];
For[i = 2, i <= ngrid - 1, i++,
v[[i]] = vNew[[i]];
Σ[[i]] = 2/3 v[[i]]/(ν[[i]]*Sqrt[r[[i]]]);
];
t = t + dt;
If[TrueQ[t >= tnext],
Print["Writing a slice at t = ", t/yr2sec];
(*For[i=2,i≤ngrid-1,i++,
Print[{t,
r〚i〛,Σ\
〚i〛,ν〚i\
〛}];
];*)
tcounter = tcounter + 1;
tnext = tnext + tout;
];
];
Print["Wrote ", ngrid - 2, " radial results per time step."];
Print["Wrote ", tcounter - 1, " time steps."];
];


gives

In:= calc[] // AbsoluteTiming

During evaluation of In:= Writing a slice at t = 16.6278

During evaluation of In:= Writing a slice at t = 100016.

During evaluation of In:= Writing a slice at t = 200016.

During evaluation of In:= Writing a slice at t = 300015.

During evaluation of In:= Writing a slice at t = 400015.

During evaluation of In:= Writing a slice at t = 500015.

During evaluation of In:= Writing a slice at t = 600014.

During evaluation of In:= Writing a slice at t = 700014.

During evaluation of In:= Writing a slice at t = 800013.

During evaluation of In:= Writing a slice at t = 900013.

During evaluation of In:= Writing a slice at t = 1.00001*10^6

During evaluation of In:= Writing a slice at t = 1.10001*10^6

During evaluation of In:= Writing a slice at t = 1.20001*10^6

During evaluation of In:= Writing a slice at t = 1.30001*10^6

During evaluation of In:= Writing a slice at t = 1.40001*10^6

During evaluation of In:= Writing a slice at t = 1.50001*10^6

During evaluation of In:= Writing a slice at t = 1.60001*10^6

During evaluation of In:= Writing a slice at t = 1.70001*10^6

During evaluation of In:= Writing a slice at t = 1.80001*10^6

During evaluation of In:= Writing a slice at t = 1.90001*10^6

During evaluation of In:= Writing a slice at t = 2.00001*10^6

During evaluation of In:= Wrote 100 radial results per time step.

During evaluation of In:= Wrote 20 time steps.

Out= {233.787534, Null}


@Jens has provided awesome answer, thank you so much. Documentation is quite scarce in this regard. I have three questions though:

1. You have the following twice in: vNew = Table[0., {ngrid}]; [Nu] = Table[0., {ngrid}]; But if I delete one occurrence the code slows down (with the duplicity the code runs: 0.769044, without: 0.782045). I understand that the time difference is negligible, yet it exists and it confuses me why.
2. I do not really understand the interpretation of input by Mathematica, is 2 an integer by default? Because if so, I can understand the need to use 2. as you have in the code. On the other hand in some places, you do not use 2.0 but only 2.
3. What is the advantage of using Table instead of ConstantArray at the beginning?
4. Next is a more general question, in one of the comments above, @belisarius mentioned to avoid procedural loops, so I read some text about it where Mathematica authors suggest to use Map etc. But I fail to understand how that helps; in essence they are also loops, even when they can be vectorized or parallelized I doubt they would be able to make 2 orders of magnitude of processing time speed up.
5. When I enable the line with Print["Writing a slice at t = ", t]; I get uppon run CompiledFunction::cflist: Nontensor object generated; proceeding with uncompiled evaluation. So I dug a bit deeper and right after the line where I initialize the t variable I put:

t = 0.0; Print[t];

And now,

In:= AbsoluteTiming[calc[]]

During evaluation of In:= t

During evaluation of In:= CompiledFunction::cflist: Nontensor object generated; proceeding with uncompiled evaluation. >>

During evaluation of In:= 0.

During evaluation of In:= Writing a slice at t = 5.24607*10^8

Out= \$Aborted


My conclusion is that for whatever reason there is some sort of "dry run" when Mathematica refuses to do what it is told - set t to 0.0 (and this is not the first time, I came across this strange behaviour with LogPlot), only after that the code runs as it should. This would also explain the necessity to use TrueQ - see below.

@shrx the TrueQ was necessary because Mathematica refused to compile if I had While[t<=tmax,...], it kept complaining. I do not understand, though, how Mathematica can tell during compilation time whether the expression is or is not evaluated.

• "Avoid procedural loops" is the first commandment. See for example mathematica.stackexchange.com/a/18396/193 – Dr. belisarius Mar 19 '15 at 3:33
• Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) 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. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Dr. belisarius Mar 19 '15 at 3:34
• You should expect a procedural solution for this sort of problem to run 2-3 orders of magnitude slower than "normal" Mathematica code – Mike Honeychurch Mar 19 '15 at 3:48
• Can't look at the details, but one thing jumps out in the first line: fGetNu gets re-compiled every single time it's invoked because it's defined with SetDelayed. Probably if you state the problem you want to solve, someone can come up with a better method using built-in functionality. – Jens Mar 19 '15 at 3:50
• You may find this post interesting. – xzczd Mar 19 '15 at 5:35

The slowness is due to the fact that several steps in your code were not compilable because they invoked MainEvaluate. I localized all the variables by adding Module. Then, several variables were mis-recognized as integer when they should be reals. To fix this I added decimal points to some numbers like 1.

The externally defined functions at the beginning are best included in the compilation of the main function calc by inlining, so I removed their individual compilation and instead specified CompilationOptions->{"InlineExternalDefinitions"->True}. I also removed the SetDelayed mentioned in my comment.

When testing this in version 8, I also noticed that ConstantArray isn't compilable, so I replaced it by Table. All the Print statements needed to be removed because they can't be compiled, either.

Finally, I added the option to compile to C. There are several additional cosmetic changes one could make, but with the above changes the whole computation indeed takes less than a second now:

ClearAll["Global*"];
(*define viscosity function*)

yr2sec = Function[x, 3.155*10.^7*x];
AU2cm = Function[x, 1.496*10.^13*x];
fGetNu = Function[r, 1.34*10.^14*(r/AU2cm[1.])^(3/2)];

calc = Compile[{},
Module[{rin, rout, tmax, tout, tcounter, c0, v, t = 0.,
au0, Σ, r, vNew, ngrid, ν, dX, i, tnext, dtmin,
dtzone, x, dt}, rin = AU2cm[0.1];
au0 = AU2cm[1.];
rout = AU2cm[100.];
tmax = yr2sec[2*10.^6];
tout = yr2sec[1*10.^5];
ngrid = 102;
x = Table[0., {ngrid}];
v = Table[0., {ngrid}];
r = Table[0., {ngrid}];
Σ = Table[0., {ngrid}];
vNew = Table[0., {ngrid}];
ν = Table[0., {ngrid}];
tcounter = 0;
c0 = 0.5;
(*main part*)(*setup the grid*)
dX = 2.0*(Sqrt[rout] - Sqrt[rin])/(ngrid - 2);
For[i = 1, i <= ngrid, i++, x[[i]] = 2 Sqrt[rin] + (i - 1.5)*dX;
r[[i]] = 0.25*x[[i]]^2;];
(*Define zero-torque (Sigma=
0) boundary conditions*)Σ[] =
0; Σ[[ngrid]] = 0.;
(*Define initial conditions*)
For[i = 2, i <= ngrid - 1, i++,
If[TrueQ[r[[i]] <= AU2cm], Σ[[i]] =
1000.*(r[[i]]/au0)^(-3./2), Σ[[i]] = 1*10.^-4];
ν[[i]] =
fGetNu[r[[i]]];(*1.34*10.^14*(r[[i]]/(1.496*10.^13))^(3/2);*)
v[[i]] = 3./2 ν[[i]]*Σ[[i]]*Sqrt[r[[i]]];];
(*initial time*)t = 0.;
tnext = 0.;
(*main loop*)While[TrueQ[t <= tmax], dtmin = 1.*10.^20;
For[i = 2, i <= ngrid - 1,
i++,(*this is what is in Armitage code,
but NB.he uses explicit FTCS and this breaks von Neumann \
condition for stability,which has strict inequality,
also the further from the limit,the better*)
dtzone = c0*(dX*x[[i]])^2/(24*ν[[i]]);
dtmin = Min[dtmin, dtzone]];
dt = dtmin;
(*Evolve one timestep,
store result in array Vnew*)(*This is FTCS time step*)
For[i = 2, i <= ngrid - 1, i++,
vNew[[i]] =
v[[i]] +
12*ν[[i]]*
dt (v[[i - 1]] - 2*v[[i]] + v[[i + 1]])/(x[[i]]*dX)^2;];
For[i = 2, i <= ngrid - 1, i++, v[[i]] = vNew[[i]];
Σ[[i]] = 2/3 v[[i]]/(ν[[i]]*Sqrt[r[[i]]]);];
t = t + dt;
If[TrueQ[t >= tnext],(*Print["Writing a slice at t = ",
t/(3.155*10.^7)];*)(*For[i=2,i<=ngrid-1,i++,Print[{t,
r[[i]],Σ[[i]],ν[[i]]}];];*)
tcounter = tcounter + 1;
tnext = tnext + tout]];
Σ], CompilationTarget -> "C",
CompilationOptions -> {"InlineExternalDefinitions" -> True}];

AbsoluteTiming[calc[]]

(*
==> {0.815346, {0., 0.0562153, 0.0435909, 0.0302906, 0.0211078,
0.0150447, 0.0109952, 0.00822576, 0.00628334, 0.00488813, 0.0038641,
0.00309785, 0.00251452, 0.00206359, 0.00171019, 0.00142978,
0.0012048, 0.00102248, 0.000873367, 0.0007504, 0.000648218,
0.000562712, 0.000490701, 0.000429693, 0.000377722, 0.000333224,
0.000294942, 0.000261863, 0.000233161, 0.00020816, 0.000186305,
0.000167134, 0.000150263, 0.000135372, 0.000122192, 0.000110493,
0.000100084, 0.0000907991, 0.0000824988, 0.0000750623, 0.0000683862,
0.0000623809, 0.000056969, 0.0000520832, 0.000047665, 0.0000436631,
0.0000400329, 0.000036735, 0.0000337347, 0.0000310017,
0.0000285089, 0.0000262325, 0.0000241512, 0.0000222463,
0.0000205009, 0.0000189001, 0.0000174304, 0.00001608, 0.0000148379,
0.0000136946, 0.0000126414, 0.0000116704, 0.0000107745,
9.94737*10^-6, 9.18322*10^-6, 8.4768*10^-6, 7.82336*10^-6,
7.21858*10^-6, 6.65855*10^-6, 6.13968*10^-6, 5.65872*10^-6,
5.21271*10^-6, 4.79892*10^-6, 4.41488*10^-6, 4.05832*10^-6,
3.72716*10^-6, 3.4195*10^-6, 3.13359*10^-6, 2.86782*10^-6,
2.62072*10^-6, 2.39093*10^-6, 2.1772*10^-6, 1.97838*10^-6,
1.79341*10^-6, 1.62131*10^-6, 1.46117*10^-6, 1.31216*10^-6,
1.17352*10^-6, 1.04451*10^-6, 9.24482*10^-7, 8.12823*10^-7,
7.08963*10^-7, 6.12374*10^-7, 5.22569*10^-7, 4.39093*10^-7,
3.61523*10^-7, 2.89469*10^-7, 2.22564*10^-7, 1.6047*10^-7,
1.02871*10^-7, 4.94722*10^-8, 0.}}
*)


Here, I assumed that you need to output only the list Σ.

An important tool in analyzing where a compilation is failing is this:

Needs["CompiledFunctionTools"]

CompilePrint[calc]


The output of this diagnostic shows you the low-level commands and will alert you when a call to MainEvaluate has to be inserted. When you see that, it means that the function will not be running in compiled mode.

With the use of this tool and the List of compilable functions, you can therefore reach almost the speed of FORTRAN.

Another option that could be added, as mentioned by xzczd in the comment, is RuntimeOptions -> "Speed". This assumes that you don't expect to encounter any numerical instability due to under/overflow, e.g., so the checking for such events can be turned off.

3. Table compiles in version 8, ConstantArray doesn't.
4. There is typically no disadvantage in using For or While loops when creating a compiled function. But there is often a significant speed disadvantage with such loops when running uncompiled. You can generally make much better use of Mathematica's high-level capabilities if you don't compile your main programs, and instead reserve compilation for auxiliary functions that need to be called many times and can't be made more efficient any other way, e.g. by making them Listable.
• Since you used C compiler, RuntimeOptions -> "Speed" will speed up the code even more. – xzczd Mar 19 '15 at 5:33
• @xzczd I forgot to add that - thanks! Of course, in reality FORTRAN is unbeatable for speed, but at least Compile gets reasonable speed. – Jens Mar 19 '15 at 5:47
• @shrx It's unnecessary here, but sometimes it can be useful because without it, While could remain unevaluated if the condition doesn't evaluate to True or False. This is among the "cosmetic changes" I mentioned. – Jens Mar 19 '15 at 15:43
• You can replace the arrays by larger ones that store the information for each time step. E.g., replace Σ[[_]] by Σ[[tcount, _]] everywhere. This means you also have to start with tcount=1 and initialize Σ = Table[0., {100}, {ngrid}]. Here, 100 is a static size that should be chosen to be large enough for your number of time steps (which is variable). In the output, replace Σ by Σ[[;; tcounter - 1]]`. For more information, you should also read How to compile effectively?, which goes beyond the scope of this question. – Jens Mar 22 '15 at 4:16