Przyklad: dwie masy polaczone ,,sprężyną''

Energia potencjalna układu: suma energii grawitacyjnej dla każdej z mas + energia potencjalna sprezyny:

Ep=m1 g z1[t]+m 2g z2[t]  +k/2 ((x1[t]-x2[t])^2+(z1[t]-z2[t])^2)

Energia kinetyczna każdej z mas:

Ek = 1/2 m( Vx1[t]^2 +Vz1[t]^2) + 1/2 m( Vx2[t]^2+Vz2[t]^2)

Funkcja Lagrange'a:

L=Ek-Ep

Równania Lagrange-Eulera:

eq11 = D[D[L,Vx1[t]],t]==D[L,x1[t]]
eq12 = D[D[L,Vz1[t]],t]==D[L,z1[t]]
eq21 = D[D[L,Vx2[t]],t]==D[L,x2[t]]
eq22 = D[D[L,Vz2[t]],t]==D[L,z2[t]]

Kompletny układ równań wraz z wyrażeniami na prędkości:

sys = Union[{eq11,eq12,eq21,eq22}, {Vx1[t]==x1'[t],Vz1[t]==z1'[t],Vx2[t]==x2'[t],Vz2[t]==z2'[t]}]

m1=1

m2=3

k=6

g=10

Warunki początkowe:

ic = {x1[0]==0,z1[0]==0,x2[0]==1,z2[0]==1,Vx1[0]==10,Vz1[0]==20,Vx2[0]==4,Vz2[0]==30}

sol=NDSolve[Union[sys,ic],{x1[t],z1[t],x2[t],z2[t],Vx1[t],Vz1[t],Vx2[t],Vz2[t]},{t,0,10}]

tr1={x1[t],z1[t]} /. sol[[1]]

tr2= {x2[t],z2[t]} /. sol[[1]]

tor=ParametricPlot[{tr1,tr2},{t,0,5}]

Animate[Evaluate[Show[tor,Graphics[           Line[{tr1,tr2}]  ,PlotRange->{{0,50},{0,50}}]]] ,{t,0,5,0.01}]