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

Question

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.

Answer 1

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

Answer 2

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.

Answer 3

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.