kill(all);

Considero come punto 1 il punto ove è presente l'appoggio e sostituisco allo stesso
la reazione vincolare R: successivamento imporrò la condizione di compatibilità
con il vincolo sostituito.
Per definire le funzioni lineari, utilizzo un'ascissa curvilinea normalizzata xi, che 
vale 0 nel primo punto del tratto considerato e 1 nel punto finale.

errato: 
Mf_23: -R*a + (F-R)*b*xi;  R*a è una quota torcente

Mf_12: -R*a*xi;
Mf_23: (F-R)*b*xi;

Mt_12: 0;
Mt_23: -R*a;

integro sui vari tratti di trave, i contributi di energia potenziale elastica

U: integrate( (Mf_12^2/2/E/J + Mt_12^2/2/G/Jp) *a, xi, 0, 1)
  +integrate( (Mf_23^2/2/E/J + Mt_23^2/2/G/Jp) *b, xi, 0, 1);

calcolo cedimento al punto di applicazione di R, utilizzando Castigliano

delta_R: diff(U,R);

risolvo la reazione incognita R imponendo delta_R nullo

linsolve(delta_R=0,R),globalsolve=true;

ridefinisco U per valutazione di R secondo la nuova ipotesi di cui sopra
(saltando questo passo U potrebbe rimanere funzione di R e non della sua espressione risolta in funzione di P)

U: ev(U);

riduco a denominatore comune per snellire la forma

U:fullratsimp(U);

delta_F: diff(U,F);

dim:	[a=800,b=1000,c=600,F=1000,C=500000,
	 E=210000,G=210000/2/(1+3/10),
	 J=(40^4-36^4)*%pi/64,Jp=(40^4-36^4)*%pi/32];

ev(delta_F,dim)$
%,numer;

controllo su FEM per avere una verifica indipendente