(*

This code is written by Antun Balaz (antun@ipb.ac.rs) and Ivana Vidanovic (ivanavi@ipb.ac.rs) from the Scientific Computing Laboratory of the Institute of Physics Belgrade, http://www.scl.rs/

The code calculates single-particle one-dimensional effective potential W[x, xbar; eps] to level eps^(p-1) for time independent potential V[x]. he use of level p effective potential for calculation of short-time amplitudes leads to substantially reduces errors, proportional to eps^(p+1-Pd/2), where P is the number of particles, and d is spatial dimensionality. For the case P=d=1 considered in this file, the errors are proportional to eps^(p+1/2). In Path-Integral Monte Carlo calculations of transition amplitudes, the error is proportional to 1/N^p, where N is the number of time steps used in the discretization. The effective potential is defined as

W[x, xbar, eps] = Sum[c[m,k] eps^(m-k) xbar^(2k), {m, 0, p-1}, {k, 0, m}] ,

and the coefficients c[m, k] are calculated by iterating the recursion given in Eq. (14) from the paper "Recursive Schroedinger Equation Approach to Faster Converging Path Integrals" by A. Balaz, A. Bogojevic, I. Vidanovic, and A. Pelster, Phys. Rev. E 79, 036701 (2009), available from http://arxiv.org/abs/0806.4774

The effective potential here is calculated for the specific case of a Modified Poeschl-Teller potential,

V[x] = - c / Cosh[a x]^2

To further simplify calculations, the coefficients c[m, k] are expressed in terms of numbers d[m, k, i], where

c[m, k] = Sum[d[m, k, i] Sech[a x]^(2 m - 2 i + 2) Tanh[a x]^(2 i), {i, 0, m}]

The effective actions to level p=41 are available in the file EffectiveAction-MPT-1d-p41.m from http://www.scl.rs/speedup/

*)

p=40;
CurrentTime=FromDate[Date[]];
StartTime=CurrentTime;
d[ii_, jj_, kk_]=0;
d[0, 0, 0]=- c;
For[m = 1, m < p, m++,
d[m, m, m] = Expand[4 a^2 d[m-1, m-1, m-1] / ( 2 m (2 m+1))];
d[m, m, 0] = Expand[(- 2 a^2 m d[m-1, m-1, 0] + 2 a^2 d[m-1, m-1, 1]) / ( 2 m (2 m+1))];
If[m > 1, d[m, m, m-1] = Expand[(16 a^2 d[m-1, m-1, m-2] - a^2 (12 m-10) d[m-1, m-1, m-1]) / (2 m (2 m+1))]];
If[m > 2, For[i = 1,i ≤ m-2, i++,d[m, m, i] = Expand[(2 a^2 (2 (m-i+1)^2 d[m-1, m-1, i-1] - (4 i (m-i) + m + i) d[m-1, m-1, i] + (i+1) (2 i+1) d[m-1, m-1, i+1])) / (2 m (2 m+1))]]];
For[k = m-1, k ≥  0, k--,
For[i = 0, i ≤ m, i++,
d[m, k, i] = Expand[( 2 (k+1) (2 k+1) d[m, k+1, i] + 2 a^2 (-((m-i) (2 i+1) + i (2 m-2 i+2)) d[m-1, k, i] + 2 (m+1-i)^2 d[m-1, k, i-1] + (i+1) (2 i+1) d[m-1, k, i+1]) - Sum[4 r (k-r+1) d[l, r, n] d[m-l-1, k-r+1, i-n], {l, 1, m-2}, {r, Max[0, k-m+l+2], Min[k, l]}, {n, 0, i}] - 4 a^2 Sum[( (l+1-n) (m-l-i+n) d[l, r, n]  d[m-l-2, k-r, i-1-n] - ((l+1-n) (i-n) +  n (m-l-1-i+n)) d[l, r, n]  d[m-l-2, k-r, i-n] + n (i+1-n) d[l, r, n] d[m-l-2, k-r, i+1-n] ), {l, 0, m-2}, {r, Max[0, k-m+l+2], Min[k, l]}, {n, 0, i+1}]) / (8 (m+k+1))]]];
Print["p = ",m+1," contribution calculated."]; CalculationEndTime=FromDate[Date[]];  CalculationTime=CalculationEndTime-StartTime;
Print["Calculation time = ",CalculationTime," seconds."]; Print["Memory used: ",N[MaxMemoryUsed[]/1024/1024]," MB"]; Print[]]

(* To save the output of p=40 calculation to a file, the following command was used: *)
Definition[d] >>> "EffectiveAction-MPT-1d-p40.m";

(* To recover the coefficients c[m, k], the following command can be used: *)
For[m = 0, m < p, m++, For[k = m, k ≥ 0, k--, c[m, k] = 0; For[i = 0, i ≤ m, i++, c[m, k] += d[m, k, i] Sech[a  x]^(2 m-2 i+2) Tanh[a  x]^(2 i)]]]

(* Level p effective potential. *)
W[p, x_, xbar_, eps_] = Sum[c[m,k] eps^(m-k) xbar^(2k),{m,0,p-1},{k,0,m}]