In[1]:=

(* ROUND ONE : COMPUTE THE FIRST TERM *)(* Set i to be the top - order retai ... 371;(* Clear any variables used in method definition *)ClearAll[pb, qb, phalf] ;

Out[5]=

q[s] + (p[s] h)/m - (V^′[q[s]] h^2)/(2 m) + O[h]^4

Out[6]=

p[s] - V^′[q[s]] h - ((p[s] V^′′[q[s]]) h^2)/(2 m) + ((V^′[q[s]] V^′′[q[s]])/(4 m) - (p[s]^2 V^(3)[q[s]])/(4 m^2)) h^3 + O[h]^4

In[8]:=

(* Define the Hamiltonian here *)H0[q[s], p[s]] = p[s]^2/(2 * m) + V[q[s]] ; &#62371 ... p[s + h] = Underoverscript[∑, k = 0, arg3] ((h^k ∂_ {s, k} p[s])/k ! + O[h]^(i + 1))

Out[12]=

q[s] + (p[s] h)/m - (V^′[q[s]] h^2)/(2 m) + (-(p[s] V^′′[q[s]])/(6 m^2) + G^(0, 1)[q[s], p[s]]) h^3 + O[h]^4

Out[13]=

p[s] - V^′[q[s]] h - ((p[s] V^′′[q[s]]) h^2)/(2 m) + ((V^′[q[s]] V^′′[q[s]])/(6 m) - (p[s]^2 V^(3)[q[s]])/(6 m^2) - G^(1, 0)[q[s], p[s]]) h^3 + O[h]^4

In[14]:=

(* Expand the difference between method and exact *)DQ = Expand[Collect[q1[s] - q[s + h], h]] + O[h]^(i + 1) DP = Expand[Collect[p1[s] - p[s + h], h]] + O[h]^(i + 1)

Out[14]=

((p[s] V^′′[q[s]])/(6 m^2) - G^(0, 1)[q[s], p[s]]) h^3 + O[h]^4

Out[15]=

((V^′[q[s]] V^′′[q[s]])/(12 m) - (p[s]^2 V^(3)[q[s]])/(12 m^2) + G^(1, 0)[q[s], p[s]]) h^3 + O[h]^4

In[16]:=

(* Integrate q - diff with respect to p *)Solve[Integrate[DQ, p[s]]  O[h]^(i ... diff with respect to q *)Solve[Integrate[DP, q[s]] O[h]^(i + 1), G[q[s], p[s]]]

Out[16]=

{{G[q[s], p[s]]  (p[s]^2 V^′′[q[s]])/(12 m^2)}}

Out[17]=

{{G[q[s], p[s]]  (-m V^′[q[s]]^2 + 2 p[s]^2 V^′′[q[s]])/(24 m^2)}}

In[18]:=

(* The two expressions should be the same up to a constant in q for the first equation, and ... first perturbation *)H2[q[s], p[s]] = (2 * (p[s]/m)^2 * V''[q[s]] - (V '[q[s]])^2/m)/24

Out[18]=

1/24 (-V^′[q[s]]^2/m + (2 p[s]^2 V^′′[q[s]])/m^2)

In[19]:=

(* ROUND TWO :   COMPUTE THE SECOND TERM *)(* Set i to be the top ... 371;(* Clear any variables used in method definition *)ClearAll[pb, qb, phalf] ;

Out[23]=

q[s] + (p[s] h)/m - (V^′[q[s]] h^2)/(2 m) + O[h]^6

Out[24]=

p[s] - V^′[q[s]] h - ((p[s] V^′′[q[s]]) h^2)/(2 m) + ((V^′[q[s]] V^& ... ^2) + (p[s]^2 V^′[q[s]] V^(4)[q[s]])/(8 m^3) - (p[s]^4 V^(5)[q[s]])/(48 m^4)) h^5 + O[h]^6

In[26]:=

(* Define the Hamiltonian here *)H0[q[s], p[s]] = p[s]^2/(2 * m) + V[q[s]] ; <br />& ... = Underoverscript[∑, k = 0, arg3] ((h^k ∂_ {s, k} p[s])/k ! + O[h]^(i + 1)) 

Out[30]=

q[s] + (p[s] h)/m - (V^′[q[s]] h^2)/(2 m) + (-(p[s] V^′′[q[s]]^2)/(30 m^3) ... [s]] V^(3)[q[s]])/(60 m^3) + (p[s]^3 V^(4)[q[s]])/(180 m^4) + G^(0, 1)[q[s], p[s]]) h^5 + O[h]^6

Out[31]=

p[s] - V^′[q[s]] h - ((p[s] V^′′[q[s]]) h^2)/(2 m) + ((V^′[q[s]] V^& ... q[s]] V^(4)[q[s]])/(15 m^3) - (p[s]^4 V^(5)[q[s]])/(45 m^4) - G^(1, 0)[q[s], p[s]]) h^5 + O[h]^6

In[32]:=

(* Expand the difference between method and exact *)DQ = Expand[Collect[q1[s] - q[s + h], h]] + O[h]^(i + 1) DP = Expand[Collect[p1[s] - p[s + h], h]] + O[h]^(i + 1)

Out[32]=

((p[s] V^′′[q[s]]^2)/(30 m^3) + (p[s] V^′[q[s]] V^(3)[q[s]])/(60 m^3) - (p[s]^3 V^(4)[q[s]])/(180 m^4) - G^(0, 1)[q[s], p[s]]) h^5 + O[h]^6

Out[33]=

((V^′[q[s]] V^′′[q[s]]^2)/(120 m^2) + (V^′[q[s]]^2 V^(3)[q[s]])/(240 ... s]] V^(4)[q[s]])/(120 m^3) + (p[s]^4 V^(5)[q[s]])/(720 m^4) + G^(1, 0)[q[s], p[s]]) h^5 + O[h]^6

In[34]:=

(* Integrate q - diff with respect to p *)Solve[Integrate[DQ, p[s]]  O[h]^(i ... diff with respect to q *)Solve[Integrate[DP, q[s]] O[h]^(i + 1), G[q[s], p[s]]]

Out[34]=

{{G[q[s], p[s]] 1/(720 m^4) (12 m p[s]^2 V^′′[q[s]]^2 + 6 m p[s]^2 V^′[q[s]] V^(3)[q[s]] - p[s]^4 V^(4)[q[s]])}}

Out[35]=

{{G[q[s], p[s]] 1/(720 m^4) (-3 m^2 V^′[q[s]]^2 V^′′[q[s]] + 12 m p[s]^2 V^′′[q[s]]^2 + 6 m p[s]^2 V^′[q[s]] V^(3)[q[s]] - p[s]^4 V^(4)[q[s]])}}

In[36]:=

(* The two expressions should be the same up to a constant in q for the first equation, and ... )^2 * V '[q[s]] * V'''[q[s]]/m - (p[s]/m)^4 * V''''[q[s]] - 3 * (V '[q[s]]/m)^2 * V''[q[s]])/720

Out[36]=

1/720 (-(3 V^′[q[s]]^2 V^′′[q[s]])/m^2 + (12 p[s]^2 V^′′[q[s]]^2)/m^3 + (6 p[s]^2 V^′[q[s]] V^(3)[q[s]])/m^3 - (p[s]^4 V^(4)[q[s]])/m^4)

In[37]:=

(* ROUND THREE : VERIFY THE EXPANSION *)(* Set i to be the top - order retai ... 371;(* Clear any variables used in method definition *)ClearAll[pb, qb, phalf] ;

Out[41]=

q[s] + (p[s] h)/m - (V^′[q[s]] h^2)/(2 m) + O[h]^6

Out[42]=

p[s] - V^′[q[s]] h - ((p[s] V^′′[q[s]]) h^2)/(2 m) + ((V^′[q[s]] V^& ... ^2) + (p[s]^2 V^′[q[s]] V^(4)[q[s]])/(8 m^3) - (p[s]^4 V^(5)[q[s]])/(48 m^4)) h^5 + O[h]^6

In[44]:=

(* Define the Hamiltonian here *)H0[q[s], p[s]] = p[s]^2/(2 * m) + V[q[s]] ; <br />& ... p[s + h] = Underoverscript[∑, k = 0, arg3] ((h^k ∂_ {s, k} p[s])/k ! + O[h]^(i + 1))

Out[48]=

q[s] + (p[s] h)/m - (V^′[q[s]] h^2)/(2 m) + O[h]^6

Out[49]=

p[s] - V^′[q[s]] h - ((p[s] V^′′[q[s]]) h^2)/(2 m) + ((V^′[q[s]] V^& ... ^2) + (p[s]^2 V^′[q[s]] V^(4)[q[s]])/(8 m^3) - (p[s]^4 V^(5)[q[s]])/(48 m^4)) h^5 + O[h]^6

In[50]:=

(* Expand the difference between method and exact *)DQ = Expand[Collect[q1[s] - q[s + h], h]] + O[h]^(i + 1) DP = Expand[Collect[p1[s] - p[s + h], h]] + O[h]^(i + 1)

Out[50]=

O[h]^6

Out[51]=

O[h]^6

Created by Mathematica  (March 27, 2007) Valid XHTML 1.1!