from sage.calculus.desolvers import desolve_system_rk4 
x1,v1, x2, v2,t=var('x1 v1 x2 v2 t')
a, b = var('a b')
a=10.0
b=1.0
P=desolve_system_rk4([v1,-a*x1+b*x2,v2,-a*x2+b*x1],[x1,v1,x2,v2],ics=[0,1,0,0,0],step=0.01,ivar=t,end_points=60)

Q1=[ [t,x1] for t,x1,v1,x2,v2 in P] #Position de la masse 1 comme fonction du temps
Q2=[ [t,x2] for t,x1,v1,x2,v2 in P] #Position de la masse 2 comme fonction du temps

#Crée les affichages correspondants
R1=list_plot(Q1,plotjoined=True,rgbcolor=(0.0,1.0,0))
R2=list_plot(Q2,plotjoined=True,rgbcolor=(0.0,0.1,0))

#graphics_array([R1,R2])#Affiche x1(t) et x2(t)

Mode1=[ [t,x1+x2] for t,x1,v1,x2,v2 in P] #Le mode +  comme fonction du temps
Mode2=[ [t,x1-x2] for t,x1,v1,x2,v2 in P] #Le mode -  comme fonction du temps

#Crée les affichages correspondants
RM1=list_plot(Mode1,plotjoined=True,rgbcolor=(0.0,1.0,0))
RM2=list_plot(Mode2,plotjoined=True,rgbcolor=(0.0,0.1,0))

graphics_array([RM1,RM2])#Affiche x1(t) et x2(t)

Etot=[ [t,v1**2+v2**2+(a-b)*(x1**2+x2**2)+b*(x1-x2)**2] for
t,x1,v1,x2,v2 in P] #energie totale comme fonction du temps

REtot=list_plot(Etot,plotjoined=True,rgbcolor=(1.0,0.0,0.0),scale='semilogy')
#Création du graphe de l'énergie totale comme fonction du temps; l'axe des valeurs de l'énergie est logarithmique

E1=[ [t,v1**2+(a-b)*x1**2] for t,x1,v1,x2,v2 in P]#L'énergie de la masse 1 comme fonction du temps

E2=[ [t,v2**2+(a-b)*x2**2] for t,x1,v1,x2,v2 in P]#l'énergie de la masse 2 comme fonction du temps

RE1=list_plot(E1,plotjoined=True,rgbcolor=(0.0,1.0,0.0),scale='semilogy')
RE2=list_plot(E2,plotjoined=True,rgbcolor=(0.0,0.0,1.0),scale='semilogy')

#graphics_array([REtot,RE1+RE2])#affiche les énergies comme fonctions du temps