Improving the plotting time of a very complicated expression

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?

-
No code showing your attempt(s)? –  Ｊ. Ｍ. May 24 at 15:29
This is a small step in the code I am doing. I could put this part of the code, but it depends on some other quantities that are not relevant for the question and it would be a bit cumbersome. I decided just to ask the essence of the problem I am facing. Schematically what i do is R= (something in terms of the a_i); then define a_i[x_]:= some function of x. Plot[R,{x,xmin,xmax}]. –  lagoa May 24 at 15:32
Well, there's the obvious approach: if you have a1[x_] := (* stuff *), ... an[x_] := (* stuff *) and R[a1_, a2_, ... an_] := (* stuff *), then you can do Plot[R[a1[x], a2[x], ... an[x]], {x, xmin, xmax}]. –  Ｊ. Ｍ. May 24 at 15:35
Are there typos in your code? For instance, the period in Abs[ψorbR]^2.vecup is interpreted as a decimal point (2.) but the Dot product makes more sense here, so that R is a scalar function of ε. –  Michael E2 May 25 at 1:52
@rm-rf Wow indeed, copying "As Plain Text" mis-interprets the period. Very insidious. –  Jens May 25 at 3:20

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.

-