I'm trying to do analytic calculations in a quantum mechanic harmonic oscillator basis. Specifically I want to be able to evaluate functions of the many particle density.

I define the following functions: [definitions]

(* Harmonic Oscillator Eigenstate *)
basis[i_, z_] = Pi^(-1/4)/Sqrt[2^i Factorial[i]] Exp[-z^2/2] HermiteH[i, z];
(* Orbital *)
orb[i_, z_] = Sum[c[i, k] basis[k, z], {k, 0, t}];
(* Particle density of one orbital *)
orbitaldens[i_, z_] = Conjugate[orb[i, z]] psi[i, z];
(* Total density *)
dens[z_] = Sum[a[i] orbitaldens[i, z], {i, 1, n}];

To help Mathematica I define some assumptions beforehand: [assumptions]

$Assumptions = Element[z, Reals] && (* Real-space coordinate is real *)
               Element[{n, t}, Integers] && n > 0 && t > 0 && (* Just indices *)
               Element[c[i_, j_], Complexes] && (* Complex state vector *)
               Element[a[i_], Reals]; (* Real filling factor *)

What is the correct way to define assumptions on an unknown function? Here c[i_, j_] should take integers i, and j and return a complex number. a[i_] on the other hand should take on integer i and return a positive real number.

Unfortunately, the evaluation of the definitons takes very long if I include these assumptions. Without them the definitions block takes about a second to evaluate.

However, I need those assumptions, so Mathematica can do useful simplifications. For example Conjugate[orb[i,z]] orb[i,z]//FullSimplify should simplify to Abs[orb[i,z]]^2. But the same should work, if I only apply the complex conjugation to c[i,k] inside orb[i,z], because basis[i,z] is real. So far, this takes too long to evaluate.

Is it possible to speed this whole thing up?

Another problem I found is the codomain of the Factorial The following two expressions do not evaluate to True, neither to False. They just repeat the first parameter of Refine:

Refine[Element[Sqrt[Factorial[i]], Reals], Element[i, Integers] && i > 0]
Refine[Factorial[i] > 0, Element[i, Integers] && i > 0]

In words, the factorial of a positive integer is not necessarily positive.

What is wrong there? I don't see how the factorial of a positive integer could become negative. (In symbolic maths)

  • 1
    $\begingroup$ Element[c[i_, j_], Complexes] isn't needed, really; Mathematica assumes it always deals with complex quantities unless told otherwise. $\endgroup$ May 29, 2012 at 15:31
  • 1
    $\begingroup$ FullSimplify[i! > 0, i \[Element] Integers && Positive[i]] and FullSimplify[Sqrt[i!] \[Element] Reals, i \[Element] Integers && Positive[i]] work nicely, it seems. $\endgroup$ May 29, 2012 at 15:34
  • 1
    $\begingroup$ @J.M. I never realized you could use Positive[...] in assumptions. This might explain why it doesn't have Q at the end (Q is for "programming functions", not "math functions"). $\endgroup$
    – Szabolcs
    May 29, 2012 at 15:42
  • $\begingroup$ @Szabolcs: I know i > 0 and i >= 0 are shorter to write, but somehow it looks more readable to me to use Positive[i] and NonNegative[i] in assumptions... of course, this is a matter of taste. :) You're right, it's convenient that they do not immediately evaluate for non-numeric arguments. $\endgroup$ May 29, 2012 at 15:51
  • $\begingroup$ Thanks for your comment J.M. If I enter FullSimplify[i! > 0, i \[Element] Integers && Positive[i]] it simplifies to True, indeed. However, FullSimplify[Conjugate[Sqrt[i!]], i \[Element] Integers && Positive[i]] gives Conjugate[Sqrt[i!]] as an output. I.e. not real. Yet, FullSimplify[Sqrt[i!] \[Element] Reals, i \[Element] Integers && Positive[i]] simplifies to True. On the other hand, FunctionExpand[Conjugate[Sqrt[i!]], i \[Element] Integers && Positive[i]] results in Sqrt[Gamma[1 + i]] where the Conjugate has vanished, so Mathematica has acknowledged that it's real. $\endgroup$
    – Lemming
    May 30, 2012 at 10:46

1 Answer 1


A somewhat related question may be worth linking:

Attaching persistent assumptions to symbol definition

Maybe the following will help with the problem of how to specify assumptions about coefficients in an expansion. I did something like this in another answer, and it could be applied here to your coefficients a[i]. The idea is that once the dimension of the arrays, n, is fixed, we can make dummy arrays of that specific length for the purposes of further symbolic manipulations, and impose the reality "assumption" by using an UpValue definition.

The naming of the individual array elements works by using Unique which automatically picks an unused name (starting with "a" in this case).

n = 3;

a = Table[Unique["a"], {n}];
Map[(# /: Conjugate[#] = #) &, a];

The last line is the crucial definition, saying simply that each element of the array a is equal to its conjugate (hence: real). Now I define the superposition function with x as the variable:

f = a.Sin[2 Pi Range[n] x]

(* ==> a29 Sin[2 Pi x] + a30 Sin[4 Pi x] + a31 Sin[6 Pi x] *)

For simplicity, I chose trigonometric functions as the basis.

Before doing the absolute square, let's define another almost identical function for comparison:

b = Table[Unique["b"], {n}];

g = b.Sin[2 Pi Range[n] x]

(* ==> b32 Sin[2 Pi x] + b33 Sin[4 Pi x] + b34 Sin[6 Pi x] *)

Here I left out the reality assumption for b.

For the last function, the modulus for real argument x is

Assuming[{x ∈ Reals}, Simplify[Conjugate[g] g]]

==> Conjugate[
  b32 Sin[2 Pi x] + b33 Sin[4 Pi x] + 
   b34 Sin[6 Pi x]] (b32 Sin[2 Pi x] + b33 Sin[4 Pi x] + 
   b34 Sin[6 Pi x])

The Conjugate remains in the result. On the other hand, if we do the same for the first function, the reality of a is used:

Assuming[{x ∈ Reals}, Simplify[Conjugate[f] f]]

==> (a29 Sin[2 Pi x] + a30 Sin[4 Pi x] + 
  a31 Sin[6 Pi x])^2


The reason I associated the realness condition with a using TagSet instead of using assumptions is that some tests don't recognize assumptions: see for example

Assuming[{c ∈ Reals}, Simplify[HermitianMatrixQ[{{c}}]]]

which should yield True but doesn't (because it does a literal comparison). This may not be a problem in your application, but it can't hurt to try and avoid such issues by associating the desired properties more closely with the symbol representing your coefficients or functions.

As an alternative, you could make the above definition work with the coefficients b by defining the function

h = Re[b].Sin[2 Pi Range[n] x]

where Re is put in explicitly to wrap the coefficients.

  • $\begingroup$ Thanks for your answer. I tried it out, however, there is one problem. I do not really want to fix the number of basis functions I'm using. However, Table complains when I try your solution without specifying a value of n. Table::iterb: "Iterator {n} does not have appropriate bounds." $\endgroup$
    – Lemming
    May 30, 2012 at 10:41
  • 1
    $\begingroup$ However, the following seems to work: $Assumptions = v[i_] \[Element] Reals; Conjugate[v[k]] // Simplify (* ==> v[k] *) $\endgroup$
    – Lemming
    May 30, 2012 at 10:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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