# How do I make a definition/function/method as I would in Java/MATLAB/Python

I've been trying to figure this out for a few hours now and it's driving me nuts. I am trying to confirm that, after transforming a Schrodinger equation, that the new wavefunctions satisfy the new potential.

In pseudocode, I'm trying to do something like this (Here I use a Hermite Polynomial / Harmonic Oscillator Example):

(*For the mth transformation, here let m=0*)
m = 0

(*Define the wavefunction*)
wavefunction[x_,n_]:= (-1)^n*Exp[1/4 x^2]*D[Exp[-1/2 x^2], {x, n}]

(*Define the potential*)
potential[x_] := 1/4 x^2

(*Define the energy eigenvalues*)
energy[n_] := n + 1/2

(*Call a function to check if the first m+10 eigenvalues work*)
ConfirmSolution[wavefunction[x, n], potential[x], energy[n], m]
(*When called, the function will return TRUE or FALSE when evaluated*)

(*Somewhere in my code, the following function exists: *)
Confirm[wavefunction[x, n], potential[x], energy[n], m] :=
For[i = 0, i <= m + 10, i++,
If[FullSimplify[-D[wavefunction[x, i], {x, 2}] + potential[x]*wavefunction[x, i]] != FullSimplify[energy[i]*wavefunction[x, i]],
(*If the Schrodinger equation is not true for ANY eigenstate 'i'*)
(*Return FALSE*)
(*Exit the function call immediately to save computation time*)
];
]:
(*If for all states 0...m+10 the Schrodinger equation is satisfied*)
(*Return TRUE*)
]


I know that I should be able to figure this out by googling it, but I'm not having any luck. Haha

So in Tex, to clarify this... it would look like:

$\psi_{n}(x) = (-1)^{n}e^{\frac{1}{4} x^{2}} \dfrac{d^{n}}{dx^{n}} ( e^{- \frac{1}{2} x^{2}} )$

$V(x) = \dfrac{1}{4} x^{2}$

$E_{n}(x) = n + \dfrac{1}{2}$

Then we would check to see if it satisfies the Schrodinger Equation for state $0$

$LHS = - \dfrac{d^{2}}{dx^{2}} \left( (-1)^{n}e^{\frac{1}{4} x^{2}} \dfrac{d^{n}}{dx^{n}} ( e^{- \frac{1}{2} x^{2}} ) \right) + \left( \dfrac{1}{4} x^{2} \right) \left( (-1)^{n}e^{\frac{1}{4} x^{2}} \dfrac{d^{n}}{dx^{n}} ( e^{- \frac{1}{2} x^{2}} ) \right) \bigg|_{n=0}$

$= \left( \dfrac{1}{2} \right) \left( e^{- \frac{1}{4} x^{2}} \right)$

and

$RHS = \left( n + \dfrac{1}{2} \right) \left( (-1)^{n}e^{\frac{1}{4} x^{2}} \dfrac{d^{n}}{dx^{n}} ( e^{- \frac{1}{2} x^{2}} ) \right) \bigg|_{n=0}$

$= \left( \dfrac{1}{2} \right) \left( e^{- \frac{1}{4} x^{2}} \right)$

So for $n=0$, $LHS = RHS$.

Then we would check to see if it satisfies the Schrodinger Equation for state $1$... then $2$... etc... and if it is satisfied for all states from $0 \rightarrow m+10$, then it would return $TRUE$. However, if for any of the states it $LHS \neq RHS$, then it would immediately abort the method call and return $FALSE$.

EDIT: Thanks for the help! I think I've almost got it. Would the following function work?

checksol[psi_, u_, e_, m_] := Module[{LHS, RHS, x},
res = True;
For[i = 0, i <= m + 10, i++,
LHS = -D[psi[x, i], {x, 2}] + v[x]*(psi[x, i]);
RHS = e[i]*(psi[x, i]);
res = TrueQ[FullSimplify[LHS - RHS] == 0];
If[Not[res], Print["Failed on eigenstate: " <> ToString[i]]; Break[]]
];
res
]


I believe that this would return True or False for all intended cases (or at least I can't seem to get it to return anything unexpected for everything I've tried so far).

I'm wondering if there is a problem defining res=True right away on the off chance that the for loop doesn't make it through the first iteration... But it would crash if that happened would it not? Would it be wise to put the for loop in a try/catch statement?

I was reading and saw that TrueQ may be a poor choice here, but I'm only concerned with the equation to be true within the entire domain of x. Is there a way to specify: check to see if the solution is true for all values of x within the domain x in (a,b)?

I did confirm that this does break out as soon as it fails, and now displays which eigenvalue it fails on.

• I'm wondering how you define something like Confirm[wavefunction[x_, n_], potential[x_], energy[n_]] in any of the languages you mentioned? It's hard to answer your question with the information you gave, could you provide a minimal example which will show your problem? My guess is that you did make definitions for wavefunction,potential and energy and because these are evaluated the definition of Confirm can never match as the results of the former will not have the Heads you are trying to match anymore... – Albert Retey Jun 12 '15 at 20:51
• @AlbertRetey I have edited my question, does this make a bit more sense now? Essentially I'm just trying to make a method that will check each eigenstate for the definitions I provide. It would be really annoying to try to set these up in the other languages I provided which is why I'm using Mathematica. But I'm a bit inexperienced with it and I can't seem to figure out how to make a method or function with appropriate return values. – Kvothealar Jun 15 '15 at 19:36

I'm still not sure what you want but I guess you just want to pass functions as arguments. That is so trivial in Mathematica that you might just not have dared to try it. A very simplistic version of your confirm function which does a check as you have described could be written like:

 confirm[psi_, u_, e_, m_] := Module[{lhs, rhs, x},
lhs = -D[psi[x, m], {x, 2}] + u[x]*(psi[x, n] /. n -> m);
rhs = e[m]*(psi[x, n] /. n -> m);
FullSimplify[lhs == rhs]
]


to make clear that the function arguments are really local I deliberatly used other names for them in the definition. You could with that definition and your definitions for wavefunction, potential and energy check the m=0 case like that:

  confirm[wavefunction, potential, energy, 0]


and e.g. test for m=0 to 20 with early exit like this:

  Do[
res = confirm[wavefunction, potential, energy, m];
If[Not[res], Break[]],
{m, 0, 20}
];
res

• Alright, I think I'm catching on. I've been playing with it a bit. I edited my post above because it is much cleaner than putting it in a comment. Do you see any problem with doing it the way I have it above? – Kvothealar Jun 16 '15 at 0:19
• it certainly is OK if it works for you. I personally would use a Do loop instead of For and probably localize res and (when using For) i as well. Other than that I think there is a typo (v vs. u in the definition of LHS). As Jens has mentioned there are more Mathematica-ish ways to do the same thing but I don't think that these will have strong advantages in this case... – Albert Retey Jun 16 '15 at 15:28

Although Do or For loops with Break[] certainly works, it looks clumsy in Mathematica. Here is what I would suggest instead:

 checksol[psi_, u_, e_, m_] := TakeWhile[
Range[0, m + 10],
Simplify[(-D[psi[x, #], {x, 2}] + u[x]*psi[x, #] ==
e[#] psi[x, #])] &]

checksol[wavefunction, potential, energy, m]

(* ==> {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10} *)


I decided to use TakeWhile instead of While because this way you get the list of indices $n$ up to which the check has succeeded. By the way, you could also try to get more general results for arbitrary n instead of iterating over a finite range. But that's not the aim of the question, I guess.

• I like this but I don't think this is what I'm looking for. I don't expect it to fail often and all I need is a true or false. One question I do have though is do you think it is dangerous to use == or TrueQ[] in a case like this? Because == will sometimes not give a True or False answer, but an equation that could be true for specific values of x. TrueQ[] just assumes that to be false in this case, but I'm a bit afraid I may get false negatives. – Kvothealar Jun 16 '15 at 13:20
• If you're more worried about false negatives than false positives, then using TrueQ would indeed be better. But my initial worry would be the opposite. If you get a false negative, you can still identify possible assumptions that you could add, and retry. With a false positive, you'll never know that there was a problem, right? Anyway, to get just one result instead of a list, you could simply test if my output above is equal to Range[0, m + 10]. – Jens Jun 16 '15 at 17:03

Thanks so much for the help guys. I figured I would come back and post what I settled on. Once I had it running properly I wanted to run it in parallel to speed it up.

(*Checkone Function checks a Schrodinger Equation for eigenstate: num*)

checkone[psi_, v_, e_, num_] := Block[{LHS, RHS, x},
LHS = -D[psi[x, num], {x, 2}] + v[x]*psi[x, num];
RHS = e[num]*psi[x, num];
If[TrueQ[FullSimplify[LHS - RHS] == 0], 1, 0;
Print["Failed on eigenstate: " <> ToString[num]]]
]

(*Checksol Function checks a Schrodinger equation for eigenstates: m=[start,end]*)

checksol[psi_, v_, e_, start_, end_] := Block[{res},
res = Parallelize[
Product[ checkone[psi, v, e, m], {m, start, end}]
];
TrueQ[res == 1]
]

V0[x_] := x^2
En[n_] := 2 n + 1
Psi0[x_, n_] := Exp[-(x)^2/2]*(-1)^n*Exp[(x)^2]*D[Exp[-(x)^2], {x, n}]
checksol[Psi0, V0, En, 0, 3]


checksol[] is called for a group of 4 eigenstates, and would Parallelize the product of checkone[] for the four states.

So if in checkone[] if TrueQ[FullSimplify[LHS - RHS] == 0] returned True, it would return 1 to checksol[], or 0 otherwise along with a print statement to say which one failed.

So if any state failed, the product would be 0 rather than 1. Then if the product was 1, checksol[] would return True, if the product was 0, it would return False.

I figured taking out the Break[] was appropriate because just running 4 states at the same time was faster than one at a time and breaking out of it.

• I think it would not be too hard to combine parallel evaluation with an early exit strategy. If you need that, you could use a combination of ParallelSubmit and WaitNext – Albert Retey Jun 17 '15 at 17:55