### Mouvement conservatif 1D: lien entre énergie potentielle et portrait de phase¶

In [1]:
Ep=-(-x^3+x^2+3*x-1)*exp(-x);Ep

Out[1]:
(x^3 - x^2 - 3*x + 1)*e^(-x)
In [2]:
plot(Ep,x,0,14,axes_labels=['$x$',r'$E_p(x)$'],fontsize=16)

Out[2]:
In [3]:
f=-diff(Ep,x) # force

In [4]:
solve(f==0,x) #positions d'équilibres

Out[4]:
[x == 4, x == -1, x == 1]
In [5]:
var('t','v')

Out[5]:
(t, v)
In [6]:
sol=[]
for i in range(10) :
sol.append(desolve_system_rk4([v,f],[x,v],ics=[0,i*1,0],end_points=[-10,10],step=1/20,ivar=t))#masse unite

In [7]:
graph=Graphics()
for i in range(10):
graph+=line([(sol[i][j][1],sol[i][j][2]) for j in range(len(sol[i]))], color='blue')
show(graph,xmin=0,xmax=11,axes_labels=['$x$',r'$\dot{x}$'],fontsize=16)

In [8]:
graph_v=plot_vector_field((v,f),(x,0,11),(v,-1.9,1.9),color='magenta')

In [9]:
show(graph_v,axes_labels=['$x$',r'$\dot{x}$'],fontsize=16)

In [10]:
show(graph+graph_v,xmin=0,xmax=11,axes_labels=['$x$',r'$\dot{x}$'],fontsize=16)


On ajoute la courbe correspondant à l'équilibre instable.

In [11]:
sol4=desolve_system_rk4([v,f],[x,v],ics=[0,4,0.01],end_points=[-39,16],step=1/20,ivar=t)
graph4=line([(sol4[j][1],sol4[j][2]) for j in range(len(sol4))], color='green',thickness=2)
show(graph4,xmin=0,xmax=11,axes_labels=['$x$',r'$\dot{x}$'],fontsize=16)

In [12]:
plot(Ep,x,0,11,axes_labels=['$x$',r'$E_p(x)$'],fontsize=16)

Out[12]:
In [13]:
show(graph+graph_v+graph4,xmin=0,xmax=11,axes_labels=['$x$',r'$\dot{x}$'],fontsize=16)


Zoom autour de la position d'équilibre stable

In [14]:
sol1=[]
for i in range(9) :
sol1.append(desolve_system_rk4([v,f],[x,v],ics=[0,(0.1+i*0.1),0],end_points=[0,10],step=1/20,ivar=t))
graph1=Graphics()
for i in range(9):
graph1+=line([(sol1[i][j][1],sol1[i][j][2]) for j in range(len(sol1[i]))], color='blue')
graph_v0=plot_vector_field((v,f),(x,0,3.5),(v,-1.9,1.9),color='magenta')
show(graph1+graph_v0,xmin=0,xmax=3.5,axes_labels=['$x$',r'$\dot{x}$'],fontsize=16)


Tracé du signal temporelle pour une condition initiale proche de l'équilibre x=1 (on prend x=0.9):

In [15]:
graph_eq=line([(sol1[8][i][0],sol1[8][i][1]) for i in range(len(sol1[8]))])
show(graph_eq,axes_labels=['$t$',r'$x(t)$'],fontsize=16)


Tracé du signal temporel pour un mouvement de plus grande amplitude autour de l'équilibre x=1:

(départ sans vitess initiale de x=0.2)

In [16]:
graph_heq=line([(sol1[1][i][0],sol1[1][i][1]) for i in range(len(sol1[1]))])
show(graph_heq,axes_labels=['$t$',r'$x(t)$'],fontsize=16)