From 0b29ecefb29e04dfdd3fad7c3e862d58dfdb5d6e Mon Sep 17 00:00:00 2001 From: Rainer Wittmaack Date: Thu, 8 Oct 2009 10:05:43 +0200 Subject: [PATCH] * upgraded euler to euler-maruyama * derivws() calculates a step along the probabilit axis * removed this ability from derivs() instead * as a consequence rk4() can't calculate stochastic equations anymore (and really couldn't before) * updated eul to reflect the changes in ruk --- bifandvel.py | 6 +- calculator.cpp | 349 ++++++++++++++++++++++++++++++++++++--------------------- calculator.h | 12 +- detdiff.py | 57 ++++++---- dimdialog.ui | 73 ++++++------ dimer.py | 4 +- pinv.py | 8 +- 7 files changed, 313 insertions(+), 196 deletions(-) diff --git a/bifandvel.py b/bifandvel.py index 98c8f38..10c772f 100644 --- a/bifandvel.py +++ b/bifandvel.py @@ -7,8 +7,8 @@ import matplotlib.mlab as mlab import matplotlib.axes as axe import matplotlib.text as txt -DATAX = '../rk_bifurcate.dat' -DATAV = '../rk_bif_avg_velo.dat' +DATAX = '../eu_bifurcate.dat' +DATAV = '../eu_bif_avg_velo.dat' FOUT = '../test.png' # some macros @@ -31,7 +31,7 @@ ytkl1 = ["0","L/2","L"] ylim1 = [0-ypad,L+ypad] #ylim2 = [-1-ypad,1+ypad] xlbl1 = r'' -xlbl2 = r'$\eta_1$' +xlbl2 = r'$l$' ylbl1 = r'$x(k\tau)$ mod $L$' ylbl2 = r'$v_x$' res = (1440/80,900/80) # default dpi is 80 diff --git a/calculator.cpp b/calculator.cpp index 7466385..22624cf 100644 --- a/calculator.cpp +++ b/calculator.cpp @@ -50,15 +50,19 @@ void Calculator::derive_init() file_rk_variance = "../rk_variance.dat"; file_rk_distance = "../rk_msd.dat"; file_rk_defangle = "../rk_defangle.dat"; + file_eu_variance = "../eu_variance.dat"; + file_eu_distance = "../eu_msd.dat"; + file_eu_defangle = "../eu_defangle.dat"; file_rk_dualbif = "../rk_dualbif.dat"; + file_eu_dualbif = "../eu_dualbif.dat"; v = c->valuesMap.value(1); } void Calculator::derivs(const Doub x, VecDoub_I &y, VecDoub_O &dydx, - const Doub params[], Normaldev &gauss) + const Doub params[]) { //params[0] = a; //params[1] = A; @@ -81,13 +85,12 @@ void Calculator::derivs(const Doub x, VecDoub_I &y, VecDoub_O &dydx, //params[18] = U; //static Doub G, xi, xi1, xi2, c, s, st, ax, ay, e0, e1, e2, e3, eta, l1, l2; static Doub xi, G, s, c, st, ax, ay, cp, sp, m1, m2, p1, p2; - static Doub x1, y1, x2, y2, F1x, F1y, F2x, F2y, stx, sty, xix1, xix2, xiy1, xiy2; + static Doub x1, y1, x2, y2, F1x, F1y, F2x, F2y, stx, sty; //for (Int i = 0;i<10;i++) // qDebug() << "params["< & map, const Doub & eps) void Calculator::IC_setup(VecDoub &y, Doub x0, Doub y0, Doub xdot0, Doub ydot0, Doub phi0, Doub params[]) { - Int n = y.size(); + //Int n = y.size(); const Int sys = c->valuesMap.value(1); /* @@ -600,6 +631,7 @@ void Calculator::calc2(const Int n) Doub phipi[samples][s]; VecDoub y(n,0.0); VecDoub dydx(n,0.0); + VecDoub dW(n, 0.0); VecDoub yout(n,0.0); VecDoub avg(n,0.0); @@ -619,25 +651,28 @@ void Calculator::calc2(const Int n) if (xtest >= k*tau) { k++; } - derivs(x,y,dydx,params,gauss); // sets first dydx[] for rk4() - rk4(y,dydx,x,h,yout,params,gauss); + derivs(x,y,dydx,params); // sets first dydx[] for rk4() + rk4(y,dydx,x,h,yout,params); + for (Int j=0;j= s*fromt) { switch(v) { case 1: - rkValues_x.insert(x,y[0]); + rkValues_x.insert(x,yout[0]); break; case 4: //xpi = ABS(fmod(yout[0],(2*PI))); //ypi = ABS(fmod(yout[1],(2*PI))); - rkValues_x.insert(x,y[0]); - rkValues_y.insert(x,y[1]); + rkValues_x.insert(x,yout[0]); + rkValues_y.insert(x,yout[1]); //rkValues_x.insert(x,xpi); //rkValues_y.insert(x,ypi); - rkValues_phi.insert(x,(180/PI)*ABS(fmod(y[2],(2*PI)))); - x1 = y[0] + params[7]*cos(y[2]); - y1 = y[1] + params[7]*sin(y[2]); - x2 = y[0] - params[8]*cos(y[2]); - y2 = y[1] - params[8]*sin(y[2]); + rkValues_phi.insert(x,(180/PI)*ABS(fmod(yout[2],(2*PI)))); + x1 = y[0] + params[7]*cos(yout[2]); + y1 = y[1] + params[7]*sin(yout[2]); + x2 = y[0] - params[8]*cos(yout[2]); + y2 = y[1] - params[8]*sin(yout[2]); rkValues_p1x.insert(x,x1); rkValues_p1y.insert(x,y1); rkValues_p2x.insert(x,x2); @@ -645,9 +680,6 @@ void Calculator::calc2(const Int n) break; } } - for (Int j=0;j=" << k << "*" << tau << "("<= k*tau) { hnew = k*tau-x; - derivs(x,y,dydx,params,gauss); - euler(y,dydx,x,hnew,yout); + derivs(x,y,dydx,params); + derivw(x,y,dW,params,gauss); + euler(y,dydx,dW,x,hnew,yout); k++; x += hnew; } else { - derivs(x,y,dydx,params,gauss); - euler(y,dydx,x,hh,yout); + derivs(x,y,dydx,params); + derivw(x,y,dW,params,gauss); + euler(y,dydx,dW,x,hh,yout); x += hh; } for (Int j=0;j= k*tau) { - hnew = k*tau-x; - xnew = x+hnew; - derivs(x,y,dydx,params,gauss); // sets first dydx[] for rk4() - euler(y,dydx,x,hnew,yout); - if (k >= s*fromt) { - switch(v) { - case 1: - euValues_x.insert(xnew,yout[0]); - break; - case 4: - //xpi = ABS(fmod(yout[0],(2*PI))); - //ypi = ABS(fmod(yout[1],(2*PI))); - euValues_x.insert(xnew,yout[0]); - euValues_y.insert(xnew,yout[1]); - //rkValues_x.insert(x,xpi); - //rkValues_y.insert(x,ypi); - euValues_phi.insert(xnew,yout[2]); - x1 = yout[0] + params[7]*cos(yout[2]); - y1 = yout[1] + params[7]*sin(yout[2]); - x2 = yout[0] - params[8]*cos(yout[2]); - y2 = yout[1] - params[8]*sin(yout[2]); - euValues_p1x.insert(xnew,x1); - euValues_p1y.insert(xnew,y1); - euValues_p2x.insert(xnew,x2); - euValues_p2y.insert(xnew,y2); - break; - } - } k++; } - derivs(x,y,dydx,params,gauss); - euler(y,dydx,x,eu_h,yout); // xold + derivs(x,y,dydx,params); // sets first dydx[] for rk4() + derivw(x,y,dW,params,gauss); + euler(y,dydx,dW,x,h,yout); for (Int j=0;j= s*fromt) { + switch(v) { + case 1: + euValues_x.insert(x,yout[0]); + break; + case 4: + //xpi = ABS(fmod(yout[0],(2*PI))); + //ypi = ABS(fmod(yout[1],(2*PI))); + euValues_x.insert(x,yout[0]); + euValues_y.insert(x,yout[1]); + //rkValues_x.insert(x,xpi); + //rkValues_y.insert(x,ypi); + euValues_phi.insert(x,yout[2]); + x1 = yout[0] + params[7]*cos(yout[2]); + y1 = yout[1] + params[7]*sin(yout[2]); + x2 = yout[0] - params[8]*cos(yout[2]); + y2 = yout[1] - params[8]*sin(yout[2]); + euValues_p1x.insert(x,x1); + euValues_p1y.insert(x,y1); + euValues_p2x.insert(x,x2); + euValues_p2y.insert(x,y2); + break; + } + } } // end while qDebug() << "sizeof(euValues_x) =" << euValues_x.size(); qDebug() << "sizeof(euValues_y) =" << euValues_y.size(); @@ -1042,10 +1071,10 @@ void Calculator::calc2(const Int n) IC_setup(y,x0,xdot0,y0,ydot0,phi0,params); bifparam(v,param1,r,params); // inner for: iterate over runge-kutta - Int k = 0; - Doub hnew; + k = 0; L = 2*PI; tau = 2*PI/params[4]; + Doub hnew; //qDebug() << "tau: " << tau << "\tL: " << L; while (k= s*fromt) { - switch(v) { - case 1: - xpi[z][k] = ABS(fmod(yout[0],(2*PI))); - euValues_x.insert(r,xpi[z][k]); + derivs(x,y,dydx,params); + derivw(x,y, dW, params,gauss); + euler(y,dydx,dW,x,hnew,yout); + switch(v) { + case 1: + xpi[z][k] = yout[0]; + if (k >= s*fromt) { + euValues_x.insert(r,ABS(fmod(xpi[z][k],L))); euValues_xdot.insert(r,dydx[0]); - break; - case 4: - xpi[z][k] = ABS(fmod(yout[0],(2*PI))); - ypi[z][k] = ABS(fmod(yout[1],(2*PI))); - euValues_x.insert(r,xpi[z][k]); - euValues_y.insert(r,ypi[z][k]); - euValues_phi.insert(r,yout[2]); + } + break; + case 4: + xpi[z][k] = yout[0]; + //qDebug() << "xpi["<= s*fromt) { + euValues_x.insert(r,ABS(fmod(xpi[z][k],L))); + euValues_y.insert(r,ABS(fmod(ypi[z][k],L))); + euValues_phi.insert(r,ABS(fmod(phipi[z][k],L))); euValues_xdot.insert(r,dydx[0]); euValues_ydot.insert(r,dydx[1]); euValues_phidot.insert(r,dydx[2]); - break; } - } + break; + } // end switch k++; } - derivs(x,y,dydx,params,gauss); - euler(y,dydx,x,eu_h,yout); // xold + derivs(x,y,dydx,params); + derivw(x,y,dW, params,gauss); + euler(y,dydx,dW,x,eu_h,yout); // xold for (Int j=0;j trajectory_length) { + velocity_angle = PI; + } else { + velocity_angle = (acos((yout[0]-x0)/trajectory_length)); + } + bias_angle = params[10]; + deflection = (180/PI)*ABS(velocity_angle - bias_angle); euAvgvelo_x.insert(r,avg[0]); euAvgvelo_y.insert(r,avg[1]); euAvgvelo_phi.insert(r,avg[2]); + euDefAngle.insert(r,deflection); break; } // end switch } // samples loop + VecDoub msd_x(samples,0.0); + VecDoub msd_y(samples,0.0); + VecDoub msd_phi(samples,0.0); + Doub sigma_x[s-1]; + Doub sigma_y[s-1]; + Doub sigma_phi[s-1]; + //qDebug() << "samples:" << samples << "steps:" << s; + for(Int i=0;i=0;k--) { + euVariance_x.insert(r,sigma_x[k]); + euVariance_y.insert(r,sigma_y[k]); + euVariance_phi.insert(r,sigma_phi[k]); + } + //qDebug() << euVariance_x; } // end outer for /* * we need to do something about this @@ -1112,18 +1190,23 @@ void Calculator::calc2(const Int n) write_to_file(file_eu_bif, euValues_x, euValues_y, euValues_phi, params); write_to_file(file_eu_biv, euValues_xdot, euValues_ydot, euValues_phidot, params); write_to_file(file_eu_avg, euAvgvelo_x, euAvgvelo_y, euAvgvelo_phi, params); + write_to_file(file_eu_variance, euVariance_x, euVariance_y, euVariance_phi,params); + write_to_file(file_eu_distance, euDist_x, euDist_y, euDist_phi, params); + write_to_file(file_eu_defangle, euDefAngle, NullMap, NullMap, params); break; case 3: qDebug() << "dual bifurcation mode"; qDebug() << "bifurcation parameters: " << param1<<"\t"< NullMap; // dummy var for use with write_to_file() @@ -86,7 +91,8 @@ private: file_rk_biv, file_rk_conv, file_rk_variance, file_rk_distance, file_rk_defangle, file_rk_dualbif; QString file_eu_bif, file_eu_avg, file_eu_out, file_eu_pt1, file_eu_pt2, - file_eu_biv, file_eu_conv, file_eu_defangle, file_eu_dualbif; + file_eu_biv, file_eu_conv, file_eu_variance, file_eu_distance, + file_eu_defangle, file_eu_dualbif; Int v; Doub T; diff --git a/detdiff.py b/detdiff.py index 9f6da20..161833d 100644 --- a/detdiff.py +++ b/detdiff.py @@ -8,18 +8,20 @@ import random from pylab import * from matplotlib.pyplot import * -FIN = '../rk_variance.dat' -FOUT = '../diffusion.png' +FVAR = '../data/rk_variance.dat' +#FMEAN = '../rk_mean.dat' +FOUT = '../latex/images/diffusion.png' L = 2*np.pi; tau = L/0.1; samples = 3 -xlabel1 = r'$k\,[\tau]$' -ylabel1 = r'$\langle {x(k\tau)}^2 - \langle x(k\tau) \rangle^2\rangle$' +ylabel1 = r'$k\,[\tau]$' +#ylabel1 = r'$\langle {x(k\tau)}^2 - \langle x(k\tau) \rangle^2\rangle$' +xlabel1 = r'$\langle x \rangle$' xlabel2 = r'$k\,[\tau]$' ylabel2 = r'$D$' -res = (1440/80,900/80) # default dpi is 80 +res = (1024/80,768/80) # default dpi is 80 def get_trajectory(r,x): print "---------------------------------------------------------" @@ -40,24 +42,31 @@ def get_trajectory(r,x): print "reading data..." -data = mlab.csv2rec(FIN, comments='#', delimiter='\t') -r = data.r -x = data.x -y = data.y +var = mlab.csv2rec(FVAR, comments='#', delimiter='\t') +#mean = mlab.csv2rec(FMEAN, comments='#', delimiter='\t') + +r = var.r +var_x = var.x +var_y = var.y +#mean_x = mean.x +#mean_y = mean.y print "done." print "extracting trajectories..." -tra_x = get_trajectory(r,x) -tra_y = get_trajectory(r,y) +vartra_x = get_trajectory(r,var_x) +vartra_y = get_trajectory(r,var_y) +#meantra_x = get_trajectory(r,mean_x) +#meantra_y = get_trajectory(r,mean_y) print "done." -t = range(len(tra_x[0])) +t = range(len(vartra_x[0])) print t +unique_r = unique(r) fig = plt.figure(figsize=res) -ax1 = fig.add_subplot(121) -ax2 = fig.add_subplot(122) +#ax1 = fig.add_subplot(121) +ax2 = fig.add_subplot(111) D_x =[ [] for DUMMYVAR in range(len(unique(r))) ] @@ -73,31 +82,37 @@ D_y =[ [] for DUMMYVAR in range(len(unique(r))) ] for i in range(len(unique(r))): for j in t: - tmp = tra_x[i][j]/(2*(j+1)*tau) + tmp = vartra_x[i][j]/(2*(j+1)*tau) #print tmp,"=",tra_x[i][j],"/",2*(j+1)*tau D_x[i].append(tmp) for i in range(len(unique(r))): for j in t: - tmp = tra_y[i][j]/(2*(j+1)*tau) + tmp = vartra_y[i][j]/(2*(j+1)*tau) D_y[i].append(tmp) #print tra_x #print D_x -for i in range(len(unique(r))): - ax1.plot(t,tra_x[i]) +#for i in range(len(unique(r))): +# ax1.plot(meantra_x[i],t) for i in range(len(unique(r))): ax2.plot(t,D_x[i]) # ax2.plot(t,D_y[i]) +for i in range(len(unique(r))): + #tmp_tra = meantra_x[i] + tmp_D = D_x[i] + #ax1.text(tmp_tra[max(t)],max(t)+3,unique_r[i],) + ax2.text(max(t)+3,tmp_D[max(t)],unique_r[i],) + #ax1.set_xlim(0,len(tra_x[0])) -ax1.set_ylim(0,14001) -ax1.set_xlabel(xlabel1) -ax1.set_ylabel(ylabel1) +#ax1.set_ylim(0,14001) +#ax1.set_xlabel(xlabel1) +#ax1.set_ylabel(ylabel1) #ax2.set_xlim(0,len(tra_x[0])) ax2.set_ylim(0,2) ax2.set_xlabel(xlabel2) diff --git a/dimdialog.ui b/dimdialog.ui index e9cf2b0..c83089c 100644 --- a/dimdialog.ui +++ b/dimdialog.ui @@ -285,9 +285,9 @@ 170 - 170 - 101 - 23 + 173 + 261 + 20 @@ -335,6 +335,9 @@ 22 + + 3 + @@ -378,6 +381,38 @@ USpinBox + + + + 120 + 120 + 91 + 23 + + + + 6 + + + 9999.989999999999782 + + + + + + 100 + 120 + 16 + 18 + + + + T + + + tempSpinBox + + @@ -647,38 +682,6 @@ phiSpinBox - - - - 50 - 150 - 221 - 23 - - - - 10 - - - 9999.989999999999782 - - - - - - 10 - 150 - 31 - 18 - - - - T - - - tempSpinBox - - diff --git a/dimer.py b/dimer.py index 3143a50..eb4119d 100644 --- a/dimer.py +++ b/dimer.py @@ -11,7 +11,7 @@ import matplotlib.pyplot as plt import matplotlib.patches as patch import matplotlib.colors as col -FOUT = "../dimer_intro.png" +FOUT = "../latex/images/dimer_intro.png" L = 2*np.pi w = 0.1 @@ -26,7 +26,7 @@ cofx = cofy = L/2 pad = 0.01 # small padding so the ticks get drawn correctly... msize = 20 hw = 0.15 # arrow head width -res = (1440/80,900/80) # default dpi is 80 +res = (800/80,600/80) # default dpi is 80 cmap = cm.gray_r # colormap shader = False # applies ultra cool shadows to the map! show_plot = False diff --git a/pinv.py b/pinv.py index dd3e757..2db9d8b 100644 --- a/pinv.py +++ b/pinv.py @@ -10,14 +10,14 @@ import matplotlib.patches as patch import matplotlib.colors as col FOUT = '../dimer_in_potential.png' -FCOF = '../rk_out.dat' -FP1 = '../rk_point1.dat' -FP2 = '../rk_point2.dat' +FCOF = '../eu_out.dat' +FP1 = '../eu_point1.dat' +FP2 = '../eu_point2.dat' L = 2*np.pi w = 0.1 tau = L/w -evry = 1 # plot periodicity, don't set < 10 (bug) +evry = 1000 # plot periodicity, don't set < 10 (bug) offx = -98 # offset of (x,y) in L offy = 94 tfrom = 80 # from t -- 2.11.4.GIT