Mathematica: check if two functions are equal

[rynlee]

Hello All,

Thanks in advance for any advice.

I'm trying to evaluate whether or not a given wavefunction is a valid solution to the time-dependent schroedinger wave equation, the bottom line being that I want to check that two functions are equal to each other at all points (within a tolerance), and I'm not sure if I'm going about it the right way.

I've defined a function, and then want to evaluate two other functions that act on it, so I have:

(*define wavefunction*)
h = 6.626*10^(-34); w = 10;
psi[x_, t_] := (m*w/Pi/h)^(1/4)*
   Exp[(-m*w/2/h)*(x^2 + (a^2)/2*(1 + Exp[-2*I*w*t]) + I*h*t/m - 
       2*a*x*Exp[-i*w*t])];
(*evaluate hamiltonian operation on psi*)
H = (m/2)*Derivative[2, 0][psi][x, t] + (m/2)*w^2*x^2*psi[x, t];
(*evaluate time-depedence*)
T = I*h*Derivative[0, 1][psi][x, t];

My questions are:
1) am I taking the right approach to determining what I want to?
2) If so, how can I check for the equality of these two functions?

For (2), I tried H==T, and it just spit out the functions back, and I tried plotting H and T and showing that they are the same, but I just flatlined on both of them, so I'm not sure if my code is wrong or if those approaches are ineffective.

Thanks again!

 

[Bill Simpson]

1. Correct the typo Exp[-i*w*t] to Exp[-I*w*t].

2. Consider not using the constants, to avoid underflow, and inspect FullSimplify'd results.

In[1]:=
(*define wavefunction h=6.626*10^(-34);w=10;*)
psi[x_,t_]:=(m*w/Pi/h)^(1/4)*Exp[(-m*w/2/h)*(x^2+(a^2)/2*(1+Exp[-2*I*w*
t])+I*h*t/m-2*a*x*Exp[-I*w*t])];

In[4]:=(*evaluate hamiltonian operation on psi*)
FullSimplify[H=(m/2)*Derivative[2,0][psi][x,t]+(m/2)*w^2*x^2*psi[x,t]]//InputForm

Out[4]//InputForm=
(((m*w)/h)^(5/4)*(a^2*m^2*w - 2*a*E^(I*t*w)*m^2*w*x + E^((2*I)*t*w)*(-(h*m) + (h^2 + m^2)*w*x^2)))/(2*E^((w*((10*I)*h*t + m*(a^2*(1 + E^((-2*I)*t*w)) - (4*a*x)/E^(I*t*w) + 2*x^2)))/(4*h))*h*Pi^(1/4))

In[5]:=(*evaluate time-depedence*)
FullSimplify[T=I*h*Derivative[0,1][psi][x,t]]//InputForm

Out[5]//InputForm=
(w*((m*w)/h)^(1/4)*(E^((2*I)*t*w)*h - a^2*m*w + 2*a*E^(I*t*w)*m*w*x))/(2*E^((w*((10*I)*h*t + m*(a^2*(1 + E^((-2*I)*t*w)) - (4*a*x)/E^(I*t*w) + 2*x^2)))/(4*h))*Pi^(1/4))

Notice that a identical factor is shared between those two. Check the magnitude of that factor and then look at the magnitudes of the remaining expressions to see if the results are equal within your tolerance.

 

[rynlee]

Thanks! Wow, FullSimplify is a really powerful tool. I'm not sure I entirely understand how it works though, it generated a new variable E, how did it do that?

Also, in terms of how to use that, can I take the new Input that it spit out, and assign that to a new function of all the variables (m, w, h, a, and E)? I could then assign arbitrary values for those variables and see if the two functions are equal to within some error. Is that reasonable?

Thanks!

 

[rynlee]

I actually also got this function to work: SameQ[H[...]==T[...], I got that function to return True, which is super heartening. Is that result reliable? If it returned False I would be worried about a rounding error in there, but its return of True makes me think that it analytically compared the two, although I can't figure out how to make it show the comparison that it thinks is true.

That or perhaps it uses numerical methods but has a default tolerance, but I can't figure out how to set the tolerance to make sure it is acceptable.

 

[Bill Simpson]

For the 'E', when I want to stuff Mathematica code into a reply here so that someone else can scrape it off the screen and paste it right back into their notebook I use //InputForm. That doesn't provide the typesetting and fonts that people seem so fond of, but it will let you get the text you need. And it displays as E. In Mathematica 'E' is Euler's constant and the same thing as the other font representations of that. Try N[E] and see what it shows.

I am concerned with your use of SameQ. I believe you either use that as H===T, note the 3 equals there, or SameQ[H,T]. And since SameQ doesn't usually do any manipulation of expressions SameQ[Sin[x]^2 + Cos[x]^2, 1] or Sin[x]^2 + Cos[x]^2===1 return False. So since your two expressions appear different to me I am hesitant to think that SameQ actually gave you confidence the two are exactly the same. The documentation for SameQ claims it checks "structural equivalence" or something like that.

Certainly one way of gaining a little confidence that two complicated expressions are within some tolerance would be to assign "reasonable" values to all the variables in each, Numerically evaluate the two and see how close the difference is.

You could do
H=bigcomplicatemess;
T=biggercomplicatedmess;
H-T//.{w->somenumber,m->someothernumber,h->somenumber,a->somenumber}

And perhaps do that a few dozen or a few million times.

 

[rynlee]

You're absolutely right, I had misused SameQ, it doesn't return I want it to.

I'm pressing forward trying to inspect the coefficients of the two functions having simplified them with fullsimplify. I'm disturbed that the two functions aren't yielding identical FullSimplify results, they should be analytically identical, so I would imagine that FullSimplify would bring them to the same place ideally - this is both with and without values for the variables in. I'll continue to play with it and let you know what I manage.

 

[rynlee]

great news, it worked perfectly! I had a stupid error in my hamiltonian, but once I fixed that I was able to divide the H and T FullSimplify'd functions by each other and pull out identity! Thanks for the help, I'm really excited about using FullSimplify in the future.

 

[Bill Simpson]

Wonderful. Glad it worked for you.

Tip for the future: If you are trying to show two expressions are equal you will usually have more success showing that the difference of the two equals zero. People writing computer algebra tools often put a lot of work into being able to simplify things to zero. Dividing means denominators and those can cause all kinds of grief in calculations.