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

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

p[s] - V^′[q[s]] h + O[h]^3

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

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

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

(* 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)

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

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

(* 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]]]

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

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

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

-(3 p[s] V^′[q[s]])/m

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

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

p[s] - V^′[q[s]] h + O[h]^4

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

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

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

(* 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)

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

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

(* 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]]]

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

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

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

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

(* ROUND THREE :   COMPUTE THE THIRD TERM *)(* Set i to be the top - ... i}] (* Clear any variables used in method definition *)ClearAll[pb, qb] ;

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

p[s] - V^′[q[s]] h + O[h]^5

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

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

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

(* 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)

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

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

(* 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]]]

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

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

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

-(10 p[s] V^′[q[s]] V^′′[q[s]])/m^2

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

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

p[s] - V^′[q[s]] h + O[h]^5

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

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

p[s] - V^′[q[s]] h + O[h]^5

(* 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)

O[h]^5

O[h]^5

Created by Mathematica  (June 7, 2004)