Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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?

share|improve this question
No code showing your attempt(s)? – J. M. May 24 '13 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 '13 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}]. – J. M. May 24 '13 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 '13 at 1:52
@rm-rf Wow indeed, copying "As Plain Text" mis-interprets the period. Very insidious. – Jens May 25 '13 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},
 ψ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

Timing and plot

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:

  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.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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