Improving the plotting time of a very complicated expression

Question

The quatity $R$ is defined by a complicated expression in terms of $a_1,a_2,\ldots,a_n$. Now each $a_i$ is a function of $x$, $a_i(x)$. I would like to know what is the most efficient way to substitute in R the expression for the $a$s and then plot it as a function of $x$. I am only interested in the plot.

For concreteness consider the following code:

With[{vecup = Table[(Abs[a[i]]^2 - 1)*If[ε^2 >= 2 i, 1, 0], {i, NN}], 
      vecdown = Table[(1 - 1/Abs[a[i]]^2)*If[ε^2 >= 2 i, 1, 0], {i,  NN}]}, 
  R = Abs[ψorbR]^2.vecup/(1 + Abs[ψorbT]^2.vecdown)];
v[n_] := Sqrt[1 - 2 n/ε^2];
a[i_] := -I Sqrt[(1 + Sign[i] v[Abs[i]])/(1 - Sign[i] v[Abs[i]])];

I would like to plot R. The vectors ψorbR and ψorbT are very large expressions, and as I increase NN, these vectors increase in size causing the expressions of their elements also to increase in size, eventually making R impossible to plot.

Are there efficient ways to do operations on large expressions such as R that will help in keeping the plotting time reasonable?

Answer 1

I'm going to make some assumptions/guesses based on the code:

• The . represents the Dot product.
• ψorbR, ψorbT are each a list/vector of complex numbers of length NN.
• Since i runs from 1 to NN, Sign[i] may be omitted.
• Since a[i] is evaluated only in Abs[a[i]] the -I factor may be omitted.
• Since the terms containing a[i] are multiplied by 0 if ε^2 < 2i, two things follow for the terms included (not multiplied by 0): 0 < v[i] < 1, and the Sqrt[..] in a[i] is real; hence the Sqrt cancels with the ^2 in Abs[a[i]]^2.

Aside from the above, I make two changes that speed things up: 1) Replace v[i], a[i] by lists v0, a0 and use UnitStep instead of If[.., 1, 0]. I also wrap the function in Compile.

Block[{NN = 1000},
 SeedRandom[1];
 ψorbR = RandomReal[{-1, 1}, NN] + RandomReal[{-1, 1}, NN] I;
 ψorbT = RandomReal[{-1, 1}, NN] + RandomReal[{-1, 1}, NN] I;
 R2 = Compile[{{ε, _Real}, {ψorbR, _Complex, 1}, {ψ]orbT, _Complex, 1}},
   Module[{v2, v0, a0},
     v2 = 1 - 2 Range[Length[ψorbR]] / ε^2;
     v0 = Sqrt[v2 * UnitStep[v2]];
     a0 = (1 + v0) / (1 - v0);
     Abs[ψorbR]^2 . (a0 - 1) / (1 + Abs[ψorbT]^2 . (1 - 1/a0))
    ], RuntimeOptions -> "Speed"]
 ]

---

Comparison of timings on 100 function calls:

Do[R /. ε -> e, {e, 1., 1001, 10}]; // Timing
Do[R2[e, ψorbR, ψorbT], {e, 1., 1001, 10}]; // Timing

{2.949274, Null}  
{0.008294, Null}

Plotting it is reasonable fast:

Plot[R2[\[Epsilon], \[Psi]orbR, \[Psi]orbT], {\[Epsilon], 0, 50}] // Timing

There are insignificant differences (due to round-off error). Usually the relative error in this example is less than 10^-12:

Table[(# - R2[e, \[Psi]orbR, \[Psi]orbT])/# &[R /. \[Epsilon] -> e] // Chop[#, 10^-12] &,
  {e, 2., 1002, 10}] // Tally

{{0, 99}, {1.98163*10^-12, 1}, {1.91073*10^-11, 1}}

So if it's not exactly the function intended, at least it's an accurate representation of the posted R.

---

My interpretation of R:

With[{
  vecup = Table[(Abs[a[i]]^2 - 1) * If[ε^2 >= 2 i, 1, 0], {i, NN}], 
  vecdown = Table[(1 - 1/Abs[a[i]]^2) * If[ε^2 >= 2 i, 1, 0], {i, NN}]},
 R = Abs[ψorbR]^2 . vecup / (1 + Abs[ψorbT]^2 . vecdown)]

The functions v and a are the same.