(*

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

The code calculates imaginary-time single-particle one-dimensional effective potential W[x, xbar; eps, t] to level eps^(p-1) for time dependent potential V[x, t]. The 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[p, x, xbar, eps, t] = Sum[c[m, k] eps^(m - k) xbar^(2 k) + c[m - 1/2, k] eps^(m - 1 - k) xbar^(2 k + 1), {m, 0, p - 1}, {k, 0, m}]] ,

where the coefficients c[m, k] and c[m-1/2, k] are calculated by iterating the recursions given in Eqs. (49) and (50) from the paper "Fast Converging Path Integrals for Time-Dependent Potentials" by A. Balaz, I. Vidanovic, A. Bogojevic, and A. Pelster, available from http://arxiv.org/abs/0912.2743

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

*)

p = 10;
CurrentTime = FromDate[Date[]];
StartTime = CurrentTime;
c[_, _]=0;
For[m = 0, m < p, m++, c[m, m]=D[V[x, t], {x, 2 m}] / (2 m+1)!; For[k = m-1, k ≥ 0, k--, c[m, k]=Expand[((2 k+2) (2 k+1) c[m, k+1] + D[c[m-1, k],{x, 2}] - 
Sum[ D[c[l, r], {x, 1}] D[ c[m-l-2, k-r], {x, 1}], {l, 0, m-2}, {r, Max[0, k-m+l+2], Min[k, l]}] - Sum[ D[c[l+1/2, r], {x, 1}] D[ c[m-l-5/2, k-r-1], {x, 1}], {l, 1, m-4}, {r, Max[0, k-m+l+1], Min[k, l-1]}] -
Sum[2 r (2 k - 2 r + 2)  c[l, r] c[m-l-1, k-r+1], {l, 1, m-2}, {r, Max[0, k-m+l+2], Min[k, l]}] - Sum[(2 r +1) (2 k - 2 r + 1)  c[l+1/2, r] c[m-l-3/2, k-r], {l, 1, m-3}, {r, Max[0, k-m+l+1], Min[k, l]}] + 8 If[EvenQ[m-k], 1, 0] D[D[V[x, t], {x, 2 k}], {t, m - k}] / (2k)! /(m-k)! / 2^(m-k)) / (8 (m+k+1))]]; For[k = m, k ≥ 0, k--, c[m+1/2, k]=Expand[((2 k+3) (2 k+2) c[m+1/2, k+1] + D[c[m-1/2, k],{x, 2}] - 
Sum[ D[c[l+1/2, r], {x, 1}] D[ c[m-l-2, k-r], {x, 1}], {l, 0, m-2}, {r, Max[0, k-m+l+2], Min[k, l]}] - Sum[ D[c[l, r], {x, 1}] D[ c[m-l-3/2, k-r], {x, 1}], {l, 1, m-4}, {r, Max[0, k-m+l+1], Min[k, l-1]}] -
Sum[(2 r +1) (2 k - 2 r + 2)  c[l+1/2, r] c[m-l-1, k-r+1], {l, 1, m-2}, {r, Max[0, k-m+l+2], Min[k, l]}] - Sum[2 r (2 k - 2 r + 1)  c[l, r] c[m-l-3/2, k-r], {l, 1, m-3}, {r, Max[0, k-m+l+1], Min[k, l]}] + 8 If[OddQ[m-k], 1, 0] D[D[V[x, t], {x, 2 k+1}], {t, m - k}] / (2k+1)! /(m-k)! / 2^(m-k)) / (8 (m+k+2))]];
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=10 calculation to a file, the following command was used: *)
Definition[c] >>> "EffectiveAction-1d-Vt-p10.m";

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