2
$\begingroup$

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).

$p(x,t)$

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?

$\endgroup$
6
  • 5
    $\begingroup$ 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. $\endgroup$
    – m_goldberg
    Nov 3 '15 at 4:08
  • 1
    $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Nov 3 '15 at 5:17
  • $\begingroup$ @bbgodfrey: what do you mean by curves? $\endgroup$
    – Alex
    Nov 3 '15 at 5:21
  • $\begingroup$ ListPlot of {i, c[i]}, etc. $\endgroup$
    – bbgodfrey
    Nov 3 '15 at 5:23
  • $\begingroup$ 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. $\endgroup$
    – J. M.'s torpor
    Nov 3 '15 at 7:22
3
$\begingroup$

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

enter image description here

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 = Get["~/Downloads/zdata.m"];
w = Get["~/Downloads/wdata.m"];
c = Get["~/Downloads/cdata.m"];
{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], 
   Directive[Gray, Dashed]}, AxesStyle -> Arrowheads[0.02], 
 BaseStyle -> {FontFamily -> "Courier", FontSize -> 15}, 
 AxesLabel -> {"t", "x", "p"}]

enter image description here

$\endgroup$
3
  • $\begingroup$ That's exactly the kind of picture I wanted to get. I wonder why plotting the data as a matrix makes such a big difference? $\endgroup$
    – Alex
    Nov 3 '15 at 16:53
  • $\begingroup$ For the speedup, it was using N on your data that helped, for the plot, it has something to do with the vastly different scales for x and t, I made a post asking about it. If this answers your question go ahead and mark the check next to the answer here. $\endgroup$
    – Jason B.
    Nov 3 '15 at 17:22
  • $\begingroup$ Yes, this answers my question. Thanks again. $\endgroup$
    – Alex
    Nov 3 '15 at 17:27

Not the answer you're looking for? Browse other questions tagged or ask your own question.