This is a continuation of my previous two questions: this one and this one. I would like to plot the following function $$ p(x,t) = \frac{e^{-1/A}}{A}\sum_{i=1}^{500}e^{c_i\, t}m_i(x)\frac{z_i}{w_i}, $$ where $$ m_i(x) = \frac{x}{2}\,e^{1/x}\,W_{1,z_i/2}(2/x), $$ and $c_i$, $z_i$ and $w_i$ are all precomputed constants. Here $W_{a,b}(z)$ denotes the Whittaker W function, and in my case it is real-valued. I'm interested in the case when $x$ is between $0$ and $A=10000$, and $t$ is between $0.001$ and $1000$ (the sum is unstable for small $t$). To evaluate $p$ my strategy is to form a matrix $tt$ whose $(i,j)$-th element is $e^{c_i\,t_j}z_i/w_i$, then form a matrix $xx$ whose $(i,j)$-th element is $m_j(x_i)$, and multiply the second matrix by the first one (in that order) to get $p(x_i,t_j)$. Here's my script:
ClearAll[Evaluate[Context[] <> "*"]];
SetDirectory[NotebookDirectory[]];
$MinPrecision = 5;
A = 10000;
accuracy = $MinPrecision;
myN[x_] := N[x, accuracy];
xmin = 0;
xmax = A;
xCount = 50; (* number of points along the x axis *)
dx = FullSimplify[(xmax - xmin)/(xCount - 1)];
x = Range[xmin, xmax, dx];
tmin = 1 / 1000;
tmax = 10000;
tCount = 50; (* number of points along the t axis *)
dt = FullSimplify[(tmax - tmin)/(tCount - 1)];
t = Range[tmin, tmax, dt];
z = Get["zdata.m"];
w = Get["wdata.m"];
c = Get["cdata.m"];
zCount = Length[z];
ttfun[i_, j_] := myN[(z[[i]] * Exp[t[[j]] * c[[i]]]) / w[[i]]];
xxfun[i_, j_] := Module[{tmp = If [x[[i]] > 0, 2/x[[i]], -1]},
result =
If [tmp > 0,
myN[Power[tmp, -1] * Exp[tmp / 2] *
Re[WhittakerW[1 , z[[j]] / 2, tmp]]], 1];
result
];
Timing[tt = Table[ttfun[i, j], {i, zCount}, {j, tCount}];]
Clear[c];
(* this part is very slow because the Whittaker function is \
apparently too difficult to evaluate *)
Timing[xx = Table[xxfun[i, j], {i, xCount}, {j, zCount}];]
Clear[z];
Clear[w];
(* this matrix matrix multiplication is also very slow if xCount and \
tCount are larger than 100 *)
Timing[pp = Dot[xx, tt];]
pp = Exp[-1/A] * pp / A;
ppList = Flatten[
Table[{x[[i]], t[[j]], pp[[i, j]]}, {i, xCount}, {j, tCount}], 1];
Clear[pp];
Pause[1];
Print[ListPlot3D[ppList,
InterpolationOrder -> 1,
ColorFunction -> Function[{x, y, z}, GrayLevel[0.7 + 0.7 z]],
MeshFunctions -> {#3 &, #1 &, #2 &}, Mesh -> {20, 10, 10},
MeshStyle -> {Black, Directive[Gray, Dashed],
Directive[Gray, Dashed]},
AxesStyle -> Arrowheads[0.02],
BaseStyle -> {FontFamily -> "Courier", FontSize -> 15},
AxesLabel -> {"x", "t", "p"},
MaxPlotPoints -> 200]]
Clear[ppList];
The script requires three datafiles to run: this one , this one, and this one
The script works but suffers from the following three problems: (1) The Whittaker W function is apparently too difficult to evaluate; (2) The matrix product ($xx$ multiplied by $tt$) also takes a lot of time; and (3) Since the grid is only 50 by 50, the 3D plot the script outputs isn't smooth (see the picture below).
One way to improve the smoothness of the plot is to use a grid that's more dense. But even if the grid is 100 by 100, the script takes a lot of time to execute due to problems (1) and (2). Can some sort of interpolation be used to make the plot smoother?
c, w, z
. $\endgroup$ListPlot
of{i, c[i]}
, etc. $\endgroup$HypergeometricU[]
)? Those exponential factors aren't doing you any favors, I think. $\endgroup$