kill(all);

definisco una funzione lunghezza della molla

lunghezza(X1,Y1,X2,Y2) := sqrt( (X2-X1)^2 + (Y2-Y1)^2 );

definisco allungamento molla

dl:
 lunghezza(x_i+u_i,y_i+v_i,x_j+u_j,y_j+v_j)
-lunghezza(x_i    ,y_i    ,x_j    ,y_j    );

definisco energia potenziale elastica della molla in 
configurazione deformata

U : 1/2 * k * dl^2;

utilizzando l'analogia con l'espansione secondo taylor al secondo
ordine di una funzione scalare a variabile vettoe

vedo che K può essere definita come la matrice Hessiana del potenziale
elastico U

K : zeromatrix(4,4);

definisco un vettore (lista) delle incognite

delta : [ u_i , v_i , u_j , v_j ];

definisco una configurazione di riferimento scelta indeformata
(condizione di equilibrio in assenza di carichi esterni)

indeformata : [ u_i = 0 , v_i = 0 , u_j = 0 , v_j = 0 ];

imposto due cicli for nidificati che scorrono i termini della
matrice K

for irow : 1 thru 4 do (
  for jcol : 1 thru 4 do (
    tmp : diff (U , delta[irow] , 1 , delta[jcol] , 1 ),
    tmp : ev(tmp, indeformata),
    K[irow,jcol] : tmp
  )
);

K;