GLF-NBH.nb

In[1]:=

(* Preliminary Setup : Define Modules *)Off[General :: "spell1"] Cre ... Series[qb[h], {h, 0, i}]//.R] ... p1[h] = Simplify[Series[pb[h], {h, 0, i}]//.R]

Out[10]=

q0 + H^(0, 1)[q0, p0] h + 1/2 (-H^(0, 2)[q0, p0] H^(1, 0)[q0, p0] + H^(0, 1)[q0, p0] H^(1, 1)[ ... )[q0, p0]^2 - H^(1, 0)[q0, p0] H^(1, 2)[q0, p0] + H^(0, 1)[q0, p0] H^(2, 1)[q0, p0])) h^3 + O[h]^4

Out[11]=

p0 - H^(1, 0)[q0, p0] h + 1/2 (H^(1, 0)[q0, p0] H^(1, 1)[q0, p0] - H^(0, 1)[q0, p0] H^(2, 0)[q ... , 1)[q0, p0] (H^(1, 1)[q0, p0] H^(2, 0)[q0, p0] + H^(0, 1)[q0, p0] H^(3, 0)[q0, p0])) h^3 + O[h]^4

In[12]:=

$(* The perturbed Hamiltonian is defined here *)HP[q_, p_] = H[q, p] + h^(i - 1) * G[q, ... [q[s + h], {h, 0, i}]//.R] ... pe[h] = Simplify[Series[p[s + h], {h, 0, i}]//.R]$

Out[16]=

q0 + H^(0, 1)[q0, p0] h + 1/2 (-H^(0, 2)[q0, p0] H^(1, 0)[q0, p0] + H^(0, 1)[q0, p0] H^(1, 1)[ ... 2)[q0, p0] - H^(0, 2)[q0, p0] H^(2, 0)[q0, p0] + H^(0, 1)[q0, p0] H^(2, 1)[q0, p0])) h^3 + O[h]^5

Out[17]=

p0 - H^(1, 0)[q0, p0] h + 1/2 (H^(1, 0)[q0, p0] H^(1, 1)[q0, p0] - H^(0, 1)[q0, p0] H^(2, 0)[q ... q0, p0] + 2 H^(0, 1)[q0, p0] H^(2, 1)[q0, p0]) - H^(0, 1)[q0, p0]^2 H^(3, 0)[q0, p0]) h^3 + O[h]^5

In[18]:=

(* Expand the difference between method and exact *)DQ = Expand[Collect[q1[h] - qe[h], ... n q for the first equation, and a constant in p for the second equation *)ARE_EQ = T1 - T2

Out[20]=

${G[q, p] 1/24 (-H^(0, 2)[q, p] H^(1, 0)[q, p]^2 + 2 H^(0, 1)[q, p] H^(1, 0)[q, p] H^(1, 1)[q, p] + 2 H^(0, 1)[q, p]^2 H^(2, 0)[q, p])}$

Out[21]=

${G[q, p] 1/24 (-H^(0, 2)[q, p] H^(1, 0)[q, p]^2 + 2 H^(0, 1)[q, p] H^(1, 0)[q, p] H^(1, 1)[q, p] + 2 H^(0, 1)[q, p]^2 H^(2, 0)[q, p])}$

Out[22]=

In[23]:=

(* ROUND TWO : COMPUTE THE SECOND TERM *)(* Set i to be the top - ... Series[qb[h], {h, 0, i}]//.R] ... p1[h] = Simplify[Series[pb[h], {h, 0, i}]//.R]

Out[25]=

q0 + H^(0, 1)[q0, p0] h + 1/2 (-H^(0, 2)[q0, p0] H^(1, 0)[q0, p0] + H^(0, 1)[q0, p0] H^(1, 1)[ ... q0, p0]^3 H^(1, 0)[q0, p0] H^(3, 2)[q0, p0] + 2 H^(0, 1)[q0, p0]^4 H^(4, 1)[q0, p0])) h^5 + O[h]^6

Out[26]=

p0 - H^(1, 0)[q0, p0] h + 1/2 (H^(1, 0)[q0, p0] H^(1, 1)[q0, p0] - H^(0, 1)[q0, p0] H^(2, 0)[q ... [q0, p0] + 2 H^(2, 0)[q0, p0] H^(3, 1)[q0, p0] + H^(0, 1)[q0, p0] H^(5, 0)[q0, p0]))) h^5 + O[h]^6

In[27]:=

$(* The perturbed Hamiltonian is defined here *)HP[q_, p_] = (HP[q, p] /. T1) ... [q[s + h], {h, 0, i}]//.R] ... pe[h] = Simplify[Series[p[s + h], {h, 0, i}]//.R]$

Out[31]=

q0 + H^(0, 1)[q0, p0] h + 1/2 (-H^(0, 2)[q0, p0] H^(1, 0)[q0, p0] + H^(0, 1)[q0, p0] H^(1, 1)[ ... q0, p0]^3 H^(0, 2)[q0, p0] H^(4, 0)[q0, p0] + 16 H^(0, 1)[q0, p0]^4 H^(4, 1)[q0, p0]) h^5 + O[h]^7

Out[32]=

p0 - H^(1, 0)[q0, p0] h + 1/2 (H^(1, 0)[q0, p0] H^(1, 1)[q0, p0] - H^(0, 1)[q0, p0] H^(2, 0)[q ... p0] + 10 H^(2, 0)[q0, p0] H^(3, 1)[q0, p0] + 8 H^(0, 1)[q0, p0] H^(5, 0)[q0, p0])))) h^5 + O[h]^7

In[33]:=

(* Expand the difference between method and exact *)DQ = Expand[Collect[q1[h] - qe[h], ... n q for the first equation, and a constant in p for the second equation *)ARE_EQ = T1 - T2

Out[35]=

${G[q, p] 1/5760 (7 H^(0, 4)[q, p] H^(1, 0)[q, p]^4 + 12 H^(0, 3)[q, p] H^(1, 0)[q, p]^ ... 0)[q, p] + 32 H^(0, 1)[q, p]^3 H^(1, 0)[q, p] H^(3, 1)[q, p] - 8 H^(0, 1)[q, p]^4 H^(4, 0)[q, p])}$

Out[36]=

${G[q, p] 1/5760 (7 H^(0, 4)[q, p] H^(1, 0)[q, p]^4 + 12 H^(0, 3)[q, p] H^(1, 0)[q, p]^ ... 0)[q, p] + 32 H^(0, 1)[q, p]^3 H^(1, 0)[q, p] H^(3, 1)[q, p] - 8 H^(0, 1)[q, p]^4 H^(4, 0)[q, p])}$

Out[37]=

In[38]:=

(* ROUND THREE : VERIFY THE EXPANSION *)(* Set i to be the top - order retaine ... [q[s + h], {h, 0, i}]//.R] ... pe[h] = Simplify[Series[p[s + h], {h, 0, i}]//.R]

Out[43]=

q0 + H^(0, 1)[q0, p0] h + 1/2 (-H^(0, 2)[q0, p0] H^(1, 0)[q0, p0] + H^(0, 1)[q0, p0] H^(1, 1)[ ... q0, p0]^3 H^(1, 0)[q0, p0] H^(3, 2)[q0, p0] + 2 H^(0, 1)[q0, p0]^4 H^(4, 1)[q0, p0])) h^5 + O[h]^7

Out[44]=

p0 - H^(1, 0)[q0, p0] h + 1/2 (H^(1, 0)[q0, p0] H^(1, 1)[q0, p0] - H^(0, 1)[q0, p0] H^(2, 0)[q ... [q0, p0] + 2 H^(2, 0)[q0, p0] H^(3, 1)[q0, p0] + H^(0, 1)[q0, p0] H^(5, 0)[q0, p0]))) h^5 + O[h]^7

In[45]:=

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

Out[45]=

Out[46]=

In[47]:=

(* Final perturbed Hamiltonian *)HP[q, p] + O[h]^i

Out[47]=

H[q, p] + (-1/24 (H^(0, 2)[q, p] H^(1, 0)[q, p]^2) + 1/12 H^(0, 1)[q, p] H^(1, 0)[q, p] H^(1, ... H^(0, 1)[q, p]^3 H^(1, 0)[q, p] H^(3, 1)[q, p] - 8 H^(0, 1)[q, p]^4 H^(4, 0)[q, p]) h^4) + O[h]^5

Created by Mathematica (June 11, 2004)