# making 3d listplot smoother? [closed]

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]},
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?

• I'm voting to close this question as off-topic because it is too localized; i.e, it applies only to the local situation and needs of its poster and answers will not benefit others. Nov 3 '15 at 4:08
• Many readers will hesitate to download these files of uncertain provenance. I recommend that you add curves of the three sets of constants,c, w, z. Nov 3 '15 at 5:17
• @bbgodfrey: what do you mean by curves?
– Alex
Nov 3 '15 at 5:21
• ListPlot of {i, c[i]}, etc. Nov 3 '15 at 5:23
• Have you tried reformulating your Whittaker function in terms of the Tricomi function (HypergeometricU[])? Those exponential factors aren't doing you any favors, I think. Nov 3 '15 at 7:22

Before this gets closed, I want to point out that your code doesn't match your plot, it generates this plot:

If you replace tmax with 2.0, then you get the plot you made. Secondly, why do you use such absurd precision in your data files?

z[[2 ;; 3]]
(* {0.\
8150081379877160507750606872010435263248846948854663205945457687523565\
6162974412463360382357487168042859555429703927672801603968001671410295\
7550041837995288657848867802213546137844029126733926997287091998868501\
8260213289172021181200960560820510368192212289561400406306629951121809\
9533063583324055320155562360315353623126962096369263417512260199383489\
94841631489836974055526173012442969345155097303727 I,
1.5300490921611993167246119281091014044269324860781199195645782741357\
9526175846364154362906738587918951902162012211996875129632854348766208\
2613812493369036331017730699075449092730723674354196936762930249228943\
2515507056734471737271051021481328528876368194077590613879628686575260\
8287982960236508858376090219054758660182374286713646644495906121057261\
1932210949894174180746099993127579478941622702409068 I} *)


when this is more than sufficient?

N@z[[2 ;; 3]]
(* {0. + 0.815008 I, 0. + 1.53005 I} *)


Running the code you posted, but changing $tmax=2$ takes about 5 minutes on my machine to produce the plot you gave. Adding the line {z, w, c,A,tmin,tmax} = N /@ {z, w, c,A,tmin,tmax};, brings that down to less than 20 seconds, and gives an identical plot.

For some reason, the generated data looks like hell when you plot it as a list of tuples, but looks just fine when you plot it as a matrix.

Try this out:

A = 10000.0;

xmin = 200.0;
xmax = A;
xCount = 100;
dx = (xmax - xmin)/(xCount - 1);
x = Range[xmin, xmax, dx];

tmin = 1/20.0;
tmax = 2;
tCount = 100;
(*number of points along the t axis*)

dt = (tmax - tmin)/(tCount - 1);
t = Range[tmin, tmax, dt];

{z, w, c} = N /@ {z, w, c};

zCount = Length[z];

ttfun[i_, j_] := (z[[i]]*Exp[t[[j]]*c[[i]]])/w[[i]];

xxfun[i_, j_] :=

0.5 x[[i]] Exp[1.0/x[[i]]] Re[WhittakerW[1.0, 0.5 z[[j]], 2.0/x[[i]]]]

tt = Table[ttfun[i, j], {i, zCount}, {j, tCount}]; // AbsoluteTiming

xx = Table[xxfun[i, j], {i, xCount}, {j, zCount}]; // AbsoluteTiming
pp = Re[Exp[-1/A]/A (xx.tt)];

ListPlot3D[pp, DataRange -> {MinMax@t, MinMax@x},
ColorFunction -> Function[{x, y, z}, GrayLevel[0.7 + 0.7 z]],
MeshFunctions -> {#3 &, #1 &, #2 &}, Mesh -> {20, 10, 10},
MeshStyle -> {Black, Directive[Gray, Dashed],