(*mean velocity for the tilted particle : analytical results from the Mean First Passage Time (MFPT) approach, see correction*)

$Assumptions=p>0&&p∈Integers&&y0>0

V[y_]:=Cos[2π p y]

Integrate[Exp[-β(V[yp-y0]-V[yp])],{yp,0,1}]

TrigFactor[(V[yp-y0]-V[yp])/.p->1]

ΔV=Simplify/@%

ΔV=-2 Sin[π y0] Sin[π (0y0-2 yp)]

Integrate[Exp[-β ΔV ],{yp,0,1}]

vbar[β_,f_]:=1/(β NIntegrate[Exp[-β f y0]BesselI[0,2 β Sin[π y0]],{y0,0,∞},MaxRecursion->10,WorkingPrecision->10])

βnum=5;vbar[βnum,8]

βnum=40;
plot[βnum]=Plot[Quiet@vbar[βnum,f],{f,0,10},PlotStyle->Darker@Blue]

βnum=20;
plot[βnum]=Plot[Quiet@vbar[βnum,f],{f,0,10},PlotStyle->Blue]

βnum=10;
plot[βnum]=Plot[Quiet@vbar[βnum,f],{f,0,10},PlotStyle->Lighter@Blue]

βnum=5;
plot[βnum]=Plot[Quiet@vbar[βnum,f],{f,0,10},PlotStyle->Lighter@Purple]

βnum=5/2;
plot[βnum]=Plot[Quiet@vbar[βnum,f],{f,0,10},PlotStyle->Purple]

βnum=5/4;
plot[βnum]=Plot[Quiet@vbar[βnum,f],{f,0,10},PlotStyle->Darker@Purple]

Show[plot/@{5/4,5/2,5,10,20,40}]

(*numerical simulations of the Langevin equation in discrete time and continuous space*)

V[y_]=Cos[2π  y]

F[x_,f_]=-V'[x]+f

βnum=40;Tnum=1/βnum;

δtnum=1/100;
steps=3 10^4;
tmax=δtnum steps;
σnum=Sqrt[2 Tnum δtnum];

voff[f_]:=Module[{fnum,tabeta0},
tabeta0=RandomVariate[NormalDistribution[0,σnum],steps+1];
fnum=f;
xnum=N[x0num,10];index=1;
For[tnum=0,tnum≤tmax,tnum+=δtnum;index++,
eta0=tabeta0[[index]];
xnum=xnum+δtnum F[xnum,fnum]+eta0;
];
xnum/tmax]

Timing[voff[8]]

Timing[tv[βnum]=Table[{f,voff[f]},{f,0,10,1/2}]]

Show[ListPlot[tv[5],PlotStyle->Darker@Red],plot[5]]

Show[ListPlot[tv[10],PlotStyle->Darker@Red],plot[10]]

Show[ListPlot[tv[40],PlotStyle->Darker@Red],plot[40]]

(*From discrete master equation to Fokker Planck*)

V[x_]=Cos[2 p π x]

pplus[k_]=δt/δx^2 T Exp[-1/(2T)(V[(k+1)/L]-V[k/L]-f/L)]

pmoins[k_]=δt/δx^2 T Exp[-1/(2T)(V[(k-1)/L]-V[k/L]+f/L)]

P0[k_]=P[k/L]

(*right hand side of the master equation*)

RHS=pplus[k-1]P0[k-1]+pmoins[k+1]P0[k+1]-(pplus[k]+pmoins[k])P0[k]

(*preparation of the large L expansion at fixed x = k/L*)

RHSx=RHS/.k-> x L

Simplify[Series[RHSx,{L,∞,2}]]

Simplify[Normal@%/δt/.δx->1/L]

(*one reckognizes*)
D[(-f+V'[x])P[x],x]

(*numerical simulations of the Langevin equation in discrete time and discrete space*)

(*p=1*)
V[y_]=Cos[2π  y]

F[x_,f_]=-V'[x]+f

pplusx[x_]=δt/δx^2 T Exp[-1/(2T)(V[x+1/L]-V[x]-f/L)]

pmoinsx[x_]=δt/δx^2 T Exp[-1/(2T)(V[x-1/L]-V[x]+f/L)]

βnum=20;Tnum=1/βnum;

δtnum=10^-4;
steps=10^5;
tmax=δtnum steps;
Lnum=100;δxnum=1/Lnum;

(*check of the worst rate to jump to the right ; we lineraize |V((k+1)/L)-V(k/L)| ~1/L and takef < 20 ; the following probability has to be  <1/2 *)

N[δtnum/δxnum^2 Tnum Exp[-1/(2 Tnum)(-1/Lnum-10/Lnum)]]

(*check of the worst rate to jump to the right ; we lineraize |V((k+1)/L)-V(k/L)| ~1/L and takef < 20 ; the following probability has to be  <1/2 *)

N[δtnum/δxnum^2 Tnum Exp[-1/(2 Tnum)(-1/Lnum+0/Lnum)]]

voffdiscrete[ff_]:=Module[{fnum,tabeta0,Pplus,Pmoins,jump},
fnum=ff;
xnum=N[0,10];
For[tnum=0,tnum≤tmax,tnum+=δtnum,
Pplus=pplusx[xnum]/.{T->Tnum,δt->δtnum,δx->δxnum,f->fnum,L->Lnum};
Pmoins=pmoinsx[xnum]/.{T->Tnum,δt->δtnum,δx->δxnum,f->fnum,L->Lnum};
(*Print[N/@{Pplus,Pmoins}];*)
jump=RandomChoice[{Pplus,Pmoins,1-Pmoins-Pplus}->{δxnum,-δxnum,0}];
xnum=xnum+jump;
];
xnum/tmax]

N@Timing[voffdiscrete[8]]

vbar[βnum,8]

Timing[tvdiscrete[βnum]=N@Table[{f,voffdiscrete[f]},{f,0,10,1}]]

Timing[tvdiscrete[βnum]=N@Table[{f,voffdiscrete[f]},{f,0,10,1}]]

Timing[tvdiscrete[βnum]=N@Table[{f,voffdiscrete[f]},{f,0,10,1}]]

Timing[tvdiscrete[βnum]=N@Table[{f,voffdiscrete[f]},{f,0,10,1}]]

Timing[tvdiscrete[βnum]=N@Table[{f,voffdiscrete[f]},{f,0,10,1}]]

Timing[tvdiscrete[βnum]=N@Table[{f,voffdiscrete[f]},{f,0,10,1}]]

Timing[tvdiscrete[βnum]=N@Table[{f,voffdiscrete[f]},{f,0,10,1}]]

Timing[tvdiscrete[βnum]=N@Table[{f,voffdiscrete[f]},{f,0,10,1}]]

Timing[tvdiscrete[βnum]=N@Table[{f,voffdiscrete[f]},{f,0,10,1}]]

run[1]={{0.`,-0.114`},{1.`,0.1`},{2.`,0.102`},{3.`,0.094`},{4.`,0.124`},{5.`,0.514`},{6.`,0.748`},{7.`,3.346`},{8.`,5.144`},{9.`,7.51`},{10.`,8.51`}};

Show[ListPlot[tvdiscrete[20],PlotStyle->Darker@Red,PlotRange->All],plot[20]]

Show[ListPlot[tvdiscrete[20],PlotStyle->Darker@Red,PlotRange->All],plot[20]]

(*which is not too bad*)

(*averaging now over 10 runs*)

run[1]={{0.`,-0.114`},{1.`,0.1`},{2.`,0.102`},{3.`,0.094`},{4.`,0.124`},{5.`,0.514`},{6.`,0.748`},{7.`,3.346`},{8.`,5.144`},{9.`,7.51`},{10.`,8.51`}};

run[2]={200.237,{{0.,-0.048},{1.,-0.043},{2.,0.061},{3.,0.054},{4.,0.06},{5.,0.163},{6.,1.168},{7.,3.587},{8.,5.456},{9.,6.886},{10.,8.805}}}[[2]];

run[3]={203.329,{{0.,-0.054},{1.,0.047},{2.,0.056},{3.,0.059},{4.,0.063},{5.,0.065},{6.,1.471},{7.,3.775},{8.,5.069},{9.,7.084},{10.,8.557}}}[[2]];

run[4]={202.549,{{0.,-0.057},{1.,-0.046},{2.,0.056},{3.,0.059},{4.,0.068},{5.,0.167},{6.,1.063},{7.,3.382},{8.,5.248},{9.,7.178},{10.,8.861}}}[[2]]

run[5]={199.468,{{0.,-0.048},{1.,0.057},{2.,0.056},{3.,0.055},{4.,0.064},{5.,0.079},{6.,1.868},{7.,3.69},{8.,4.97},{9.,7.112},{10.,8.73}}}[[2]]

run[6]={202.381,{{0.,0.05},{1.,0.05},{2.,0.052},{3.,0.059},{4.,0.065},{5.,0.17},{6.,1.567},{7.,3.477},{8.,5.02},{9.,6.689},{10.,8.502}}}[[2]]

run[7]={203.849,{{0.,-0.052},{1.,0.053},{2.,0.054},{3.,0.061},{4.,0.057},{5.,0.165},{6.,1.675},{7.,3.672},{8.,5.685},{9.,7.19},{10.,8.863}}}[[2]];

run[8]={203.349,{{0.,-0.054},{1.,0.058},{2.,-0.049},{3.,0.057},{4.,0.06},{5.,0.064},{6.,0.967},{7.,3.483},{8.,5.449},{9.,7.476},{10.,8.86}}}[[2]];

run[9]={199.984,{{0.,0.049},{1.,0.045},{2.,0.052},{3.,0.055},{4.,0.064},{5.,0.169},{6.,0.875},{7.,3.473},{8.,5.372},{9.,6.986},{10.,8.665}}}[[2]];

run[10]={199.516,{{0.,-0.048},{1.,0.049},{2.,0.056},{3.,0.065},{4.,0.066},{5.,0.271},{6.,1.065},{7.,3.064},{8.,5.48},{9.,7.27},{10.,8.941}}}[[2]];

imax=10;
Show[ListPlot[Mean[run/@Array[(#1)&,imax]],PlotStyle->Darker@Green],]