# Numerical evaluations of the HamiltonianH_num=sp.lambdify(x,H,"numpy",dummify=False)# Numerical evaluations of the Hamiltonian's gradientdxH_num=sp.lambdify(x,dxH,"numpy",dummify=False)# Numerical evaluations of the structure matrixJx_num=sp.matrix2numpy(Jx,dtype=float)

defget_data(x):"return x1, x2 and H(x1, x2) for all x such that x_i is in [-limsi, limisi], i=(1, 2)"x1=[el[0]forelinx]x2=[el[1]forelinx]ifnp.max(np.abs(x1))>np.max(np.abs(lims[0])):indlimitx1=map(lambdax:x<np.max(np.abs(lims[0])),np.abs(x1)).index(False)else:indlimitx1=len(x1)ifnp.max(np.abs(x2))>np.max(np.abs(lims[1])):indlimitx2=map(lambdax:x<np.max(np.abs(lims[1])),np.abs(x2)).index(False)else:indlimitx2=len(x2)indlimitx=np.min([indlimitx1,indlimitx2])x1=x1[:indlimitx]x2=x2[:indlimitx]H=map(lambdaxk1,xk2:H_num(xk1,xk2),x1,x2)returnx1,x2,Hdefplot_H_surface(ax):''' Plot the surface H(x1, x2). axe is a 3D matplotlib.pyplot.axes (x1, x2, H) for -limsx1<x1<limsx1, -limsx2<x2<limsx2."'''frommatplotlibimportcmnpplot=200xp1=np.linspace(-lims[0],lims[0],npplot)xp2=np.linspace(-lims[1],lims[1],npplot)X1,X2=np.meshgrid(xp1,xp2)HH=np.zeros(X1.shape)foriinrange(npplot):forjinrange(npplot):HH[j,i]=H_num(xp1[i],xp2[j])surf=ax.plot_surface(X1,X2,HH,rstride=1,cstride=1,cmap=cm.BuPu,linewidth=0,antialiased=True,alpha=0.6)defplot_traj_3D(x,title=None):" 3D plot of energy and the given state trajectory, with state space as base"frommpl_toolkits.mplot3dimportAxes3Dfig=plt.figure()ax=fig.gca(projection='3d')plot_H_surface(ax)x1,x2,H=get_data(x)ax.plot_wireframe(np.array(x1),np.array(x2),H)ax.plot_wireframe(np.array(x1),np.array(x2),0)ax.scatter(np.array(x1[:1]),np.array(x2[:1]),H[:1],c='r',marker='x')ax.scatter(np.array(x1[:1]),np.array(x2[:1]),0,c='r',marker='x')ax.set_xlim(-lims[0],lims[0])ax.set_ylim(-lims[1],lims[1])# ax.set_zlim(0, 2*H[0])ax.set_xlabel('x1')ax.set_ylabel('x2')ax.set_zlabel('H(x1, x2)')iftitleisnotNone:plt.title(title)plt.savefig(title+'3D_fs'+str(fs)+'.pdf')defplot_contour():''' Plot the surface H(x1, x2). axe is a 3D matplotlib.pyplot.axes (x1, x2, H) for -limsx1<x1<limsx1, -limsx2<x2<limsx2."'''frompylabimportcontourfNticks=100.axe_x1=np.linspace(-lims[0],lims[0],Nticks)axe_x2=np.linspace(-lims[1],lims[1],Nticks)X,Y=np.meshgrid(axe_x1,axe_x2)Z=H_num(X,Y)cmap='BuPu'alpha=1NLevels=30linestyles='dash'frompylabimportcontourfcontourf(X,Y,Z,NLevels,alpha=alpha,linestyles=linestyles,antialiased=False,cmap=cmap)defplot_traj_2D(x,title=None,ax=None):"Plot the trajectory of state x in state space x1=x[0], x2=x[1]."ifax==None:fig=plt.figure()ax=plt.axes()ax.set_xlabel(r'x1')ax.set_ylabel(r'x2')ax.plot(0,0,'ko')plot_contour()ax.set_xlim(-lims[0],lims[0])ax.set_ylim(-lims[1],lims[1])x1,x2,_=get_data(x)ax.plot(x1,x2,'-b')ax.plot(x1[0],x2[0],'ro')ax.set_xlabel('x1')ax.set_ylabel('x2')iftitleisnotNone:plt.title(title)plt.savefig(title+'2D_fs'+str(fs)+'.pdf')defplot_trajs(x,title=None):plot_traj_3D(x,title=title)plot_traj_2D(x,title=title)defrel_delta_energy(x):"get sequence of H(x) and return sequence of errors on energy conservation"_,_,H=get_data(x)returnabs(np.diff(H)/H[0])defplot_delta_energy(ax=None):# time vectornt_here=int(nt/7)t=np.linspace(0,(nt_here-1)*ts,nt_here)ifax==None:ax=plt.axes()deltaH=rel_delta_energy(x_shp[:nt_here])ax.semilogy(t[:len(deltaH)],deltaH,'-',label=r'Disc. grad.',linewidth=2)deltaH=rel_delta_energy(x_midpoint[:nt_here])ax.semilogy(t[:len(deltaH)],deltaH,'--',label=r'Midpoint',linewidth=2)deltaH=rel_delta_energy(x_trapez[:nt_here])ax.semilogy(t[:len(deltaH)],deltaH,'-.',label=r'Trapezoidal',linewidth=2)deltaH=rel_delta_energy(x_euler_imp[:nt_here])ax.semilogy(t[:len(deltaH)],deltaH,':',label=r'Euler impl.',linewidth=2)deltaH=rel_delta_energy(x_euler_exp[:nt_here])ax.semilogy(t[:len(deltaH)],deltaH,'--',label=r'Euler expl.',linewidth=2)# ax.legend(loc=0)ax.set_title(r'$\epsilon_k = ({\mathrm{H}(\mathbf{x}_{k+1})-\mathrm{H}(\mathbf{x}_k)})/{\mathrm{H}(\mathbf{x}_0)}$',fontsize=18)ax.set_ylabel(r'$\epsilon_k$')ax.set_xlabel(r'time $t$ (s)')importmatplotlibmatplotlib.rcParams.update({'font.size':18})ax.legend(loc=1,fontsize=13)

frommatplotlib.pyplotimportfigure,tight_layout,savefigfig=figure()ax=fig.add_subplot(2,3,1)plot_traj_2D(x_euler_exp,title='Euler exp.',ax=ax)ax.tick_params(axis='x',# changes apply to the x-axiswhich='both',# both major and minor ticks are affected# bottom='off', # ticks along the bottom edge are off# top='off', # ticks along the top edge are offlabelbottom='off')# labels along the bottom edge are offax.set_xlabel('')ax=fig.add_subplot(2,3,2)plot_traj_2D(x_euler_imp,title='Euler imp.',ax=ax)ax.tick_params(axis='x',# changes apply to the x-axiswhich='both',# both major and minor ticks are affected# bottom='off', # ticks along the bottom edge are off# top='off', # ticks along the top edge are offlabelbottom='off')# labels along the bottom edge are offax.tick_params(axis='x',# changes apply to the x-axiswhich='both',# both major and minor ticks are affected# bottom='off', # ticks along the bottom edge are off# top='off', # ticks along the top edge are offlabelbottom='off')# labels along the bottom edge are offax.set_xlabel('')ax.set_ylabel('')ax.set_yticks(list())ax=fig.add_subplot(2,3,3)ax=fig.add_subplot(2,3,3)plot_traj_2D(x_trapez,title='Trapeze',ax=ax)ax.tick_params(axis='both',# changes apply to the x-axiswhich='both',# both major and minor ticks are affected# bottom='off', # ticks along the bottom edge are off# top='off', # ticks along the top edge are offlabelbottom='off')# labels along the bottom edge are offax.set_xlabel('')ax.set_ylabel('')ax.set_yticks(list())ax=fig.add_subplot(2,2,3)plot_traj_2D(x_midpoint,title='Point milieu',ax=ax)ax=fig.add_subplot(2,2,4)plot_traj_2D(x_shp,title='Grad. discret',ax=ax)ax.set_ylabel('')ax.set_yticks(list())#tight_layout()savefig('fs'+str(fs)+'.pdf')