Maxima 2/...

Da CdM_unimore.

Introduzione alla lezione

In questa lezione abbiamo analizzato, con l'ausilio di diversi esercizi, il calcolo delle tensioni su un anello iperstatico con il metodo di Castigliano, la trattazione esaustiva di differenti funzioni tra cui IF, THEN, ELSE, il ciclo FOR, l'approssimazione di funzioni semplici, con Taylor ed altri costrutti per iniziare ad approfondire l'utilizzo del wxmaxima.

Anello con Castigliano, introduzione teorica

Si esamina il caso di un anello sottoposto a due forze verticali uguali e contrarie tali da creare l'equilibrio alla traslazione. E' un anello chiuso, quindi per definzione è una struttura 3 volte iperstatica.

semplificazione di un anello iperstatico

Sono presenti due simmetrie, una orizzontale e una verticale. Si ricorda che, ai fini di considerare solo una parte della struttura, la simmetria deve sia essere geometrica che di carico. Nell'anello in esame abbiamo infinite simmetrie geometriche, ma solo due simmetrie di carico. Possiamo dunque studiare solo un quarto di anello. Si pone il problema del vincolamento del quarto di anello: alle due estremità non è permessa alcuna traslazione al di là degli assi, altrimenti non sarebbero rispettati i vincoli geometrici; inoltre, non è consentita alcuna rotazione, altrimenti si verrebbero a creare delle zone con sovrapposizione di materiale e altre con possibile distacco. Per far fronte a questi vincoli, è possibile utilizzare una delle due soluzioni mostrate in figura:

possibilità di scelta dei vincoli

Ricorrendo alla prima configurazione in cui sono stati utilizzati due doppi pendoli si otterrebbe un sistema iperstatico rispetto alla rotazione. Quindi, è più opportuno ricorrere alla seconda soluzione (carrello + doppio pendolo), in cui il vincolo in eccesso viene sostituito da una coppia incognita C, che andremo a calcolare a posteriori imponendo rotazione nulla in quel punto (in analogia con il bipendolo che c'era prima). E' importante ricordare che quando si vincola, bisogna evitare di rendere la struttura labile: """per esempio, considerando la seconda configurazione, se al posto del bipendolo ci fosse un doppio bipendolo il sistema diventerebbe iperstatico (secondo la formula di Grubler) ma labile in realtà, in quanto sarebbero liberi gli spostamenti verticali e orizzontali."""(punto dubbio)

: il conto dei gradi di libertà funziona solo se sono indipendenti (condizione non banale)

Equilibrio del sistema, implementazione del programma

Concio Anello.jpg

Qual è il modo più semplice per determinare il Momento flettente su tale struttura?

Si definisce una coordinata angolare θ per cui vale:

θ=0  all'apice 
θ=π⁄2 all'attacco del doppio pendolo orizzontale

In direzione verticale le due forze sono già equilibrate. Prendendo una generica sezione della struttura ed effettuando l'equilibrio del momento flettente si ottiene l'espressione di Mf in funzione dell'angolo θ:


dove è il braccio tra le due forze di modulo .

Andiamo a calcolare l’energia potenziale elastica, utilizziamo Castigliano, imponendo che la rotazione nel punto di applicazione di C sia nulla.


FIg1.jpg



Di seguito viene riportato il codice del programma:


kill(all);
Mf:C+((P/2)*r*sin(theta));
Mf (theta) := C+P/2*r*sin(theta);


Notiamo che, se non fosse già stato fatto, che il comando kill(all);, necessario per la pulizia della memoria prima di ogni programma, può essere riscritto come kill(all)$, in cui il comando $ compie l'operazione senza dare l'output (in ogni caso è bene lasciare il ";" dopo il kill(all) per sapere se effettivamente ha svolto il suo compito). Qualora non fosse stato eseguito all'avvio del programma, lo stesso può dare problemi nella soluzione del codice. Basta introdurre anche in un secondo momento il codice all'inizio della compilazione per poi digitare Ctrl+R per rieseguire l'intero codice e porre rimedio all'eventuale errore.

Continuando con l'esercizio proposto, mettiamo in evidenza anche il comando := che permette di definire una variabile sotto forma di funzione ad una o più variabili. Nel caso proposto la variabile della funzione è una sola; mentre nel caso a più variabili, queste possono essere separate da una virgola. . Inoltre da questo esempio si nota la differenza tra ":" e ":=" la prima è un'assegnazione per le espressioni, la seconda per le funzioni. Il termine tra parentesi "theta" ha solo la funzione di segnaposto: indica la presenza di una variabile. In questo modo per calcolare la funzione per determinati valori della variabile che nella realtà corrispondono a punti fisici, basta fare in questo modo: Fig2.jpg

Mf (0);
Mf(%pi/4);                                                                                                              

In alternativa la funzione può anche essere definita col comando:

define ( Mf(theta) , C+P/2*r*sin(theta)); 

Ora, a differenza della definizione attraverso il :=, la valutazione della funzione viene effettuata all'atto della definizione stessa. Invece con := la valutazione è ritardata e avviene nel momento in cui vado ad utilizzarla.

L'energia interna del sistema può essere calcolata sfruttando il teorema di Castigliano:

Anello Iperstatico Modificato.png

Nota: il segno meno che compare nella coppia C, indica che il verso supposto inizialmente è sbagliato. Si va quindi ad invertire il verso della coppia nel disegno.

NOTA: il comando linsolve di Maxima risolve (nel caso, come nell'immagine di cui sopra "linsolve(%,C);", in cui non sia specificato) la funzione\espressione come
"f(x)=0", ecco perchè non è servito esplicitarlo.

Valutando il Mf ai due estremi del quarto di anello, possiamo notare che Mf(0)<0 , Mf(%pi/2)>0. Avremo dunque un punto di flesso. Inoltre, essendo Mf monotona tra 0 e %pi/2, il momento flettente massimo (in modulo) si ha per θ=0.

U : integrate ( Mf (theta) ^2*R/2/J/E, theta, 0 , %pi/2);
rot : diff (U,C); 
linsolve ( % , C );                                                                                                                
                                                                                                   

Dopo aver calcolato l’energia interna si può calcolare il delta, ovvero lo schiacciamento dell’anello nel punto di applicazione del carico. Si deriva l’energia interna U rispetto al carico P. (Teorema Di Castigliano).

- Valutazione dello schiacciamento dell'anello sotto il carico P. Applico il :

delta:diff(2*U,P);
delta:ratsimp(ev(delta));


Tale delta, è lo schiacciamento di mezzo anello, ed è il semischiacciamento dell'anello completo o semivalizzazione dell'anello completo.

Nella formula , U è moltiplicato per 2, dato che rappresenta l’energia del quarto di anello. Così facendo si trova il delta del semi anello. Esempio su Maxima:

Max1.JPG

delta: diff(2*U,P);
delta: ratsimp (ev(delta));
PROGRAMMA COMPLETO
kill(all);
Mf:C+((P/2)*r*sin(theta));
Mf (theta) := C+P/2*r*sin(theta);
Mf (0);
Mf(%pi/4);
U : integrate ( Mf (theta) ^2*R/2/J/E, theta, 0 , %pi/2);
rot : diff (U,C); 
linsolve ( % , C );
delta: diff(2*U,P);
delta: ratsimp (ev(delta));

Approssimazione con Taylor e funzioni varie

Si presentano adesso alcuni comandi che permettono di manipolare formule o funzioni in maniera da agevolare alcune operazioni di calcolo. Si ricomincia con Kill(all); per pulire tutto ciò che è avvenuto in precedenza. -Manipolazione formule:

poly:(a+b+c)^3;


  • Forma espansa a somma di monomi:
expand(poly);

Il comando expand è utile per esempio per trovare il cubo di un trinomio, infatti eseguendo questa istruzione viene restituita come output la scomposizione di tale cubo.


  • Funzioni ratsimp e fullratsimp:
ratsimp(poly);

fullratsimp(poly);

riconduco a denominazione comune, effetuando evantauli semplificazioni


Esempio con risultati:

Max1.5.JPG

kill (all);
poly: (a+b+c)^3
expand (poly);
espr1: a/b+c/d;
espr2: fullratsimp(%);


  • Sostituzione di parti di espressioni algebriche con altre:

subst(espressione da sostituire alla vecchia (nuova) , espressione da trovare(vecchia) , espressione entro cui sostituire(ambito) );

ratsubs( espressione da sostituire alla vecchia (nuova), espressione da trovare , espressione entro cui sostituire ); “equivalente del trova/sostitiusci o find/replace”


Esempio:

Max2.JPG

subst (pippo, c/d, espr1);
subst (pippo, c/d, espr2);
ratsubst (pippo, c/d, espr2);

I due comandi sopra si comportano in maniera differente: possiamo dire che ratsubst ha un algoritmo più raffinato del comando subst. Infatti, se per esempio c/d non è esplicito, subst non da' il risultato che ci si aspetta, mentre ratsubst si (internamente è come se avesse operato in questo modo: subst + fullratsimp). Il rat opera dove ci sono delle manipolazioni di tipo razionale.


  • Trigonometria:

Applicando il seguente comando:

sin(2*x)- 2*cos(x)*sin(x);

Maxima fornisce il risultato , tuttavia il risultato corretto sappiamo essere uguale a 0 dalla trigonometria , per questo motivo utilizzo il comando per espandare "trigexpand":


Esempio:

Max3.JPG

sin (2*x) - 2*cos (x) * sin(x);
trigexpand (sin (2*x));

Oppure si può fare:


Max4.JPG

trigexpand (sin (a+2*b));
trigexpand (%);


  • Semplificazioni trigonometriche:


Max6.JPG

sin (x)^2 + cos (x)^2;
trigsimp (%);
“[oppure]: 2*sin(x)^2+cos(x)^2; =>trigsimp(%) => sin(x)^2+1”


  • Funzioni trigonometriche <==> forme esponenziali complesse:

Utilizzando il comando demoivre si trasformano tutti gli esponenziali complessi nelle corrispondenti funzioni trigonometriche:


Esempio:

Max7.JPG

a*exp (%i*b);
demoivre (%);
  • Programmazione

Si introduce l’utilizzo delle funzioni IF, ELSE, ELSEIF e THEN.

Si ‘pulisce con Kill(all);’.

if  =  se;
then = allora;
elseif = altrimenti se;
else = altrimenti.


  • Funzione segno:
sgn(x)= 1    se x>0
sgn(x)= -1   se x<0
sgn(x)= 0    se x=0

Posso scrivere:

sgn(x):= if x>0 then 1 
elseif x<0 then -1 
else 0.

oppure ancora

sgn(x):= if x>0 then 1 
elseif x<0 then -1 
elseif x=0 then 0
else ciccia (esempio sgn(0.0)=> ciccia dato che 0.0 ≠ 0

E' importante sottolineare che 0.0 è diverso da 0, infatti 0 è lo zero simbolico mentre 0.0 è lo zero numerico.

Oppure potrei scrivere:

sgn(x):= if x>0 then 1 
elseif x<0 then -1 
elseif x=0 or x=0.0 then 0
else ciccia

Di seguito viene riportata l'implementazione in Maxima.

Segno if maxima.png

Domanda: Ma se chiedessi sgn(1.5e-16)?? =>1 ?

Risposta: Conviene che il numerico test >0 e <0 sia sempre gestito con un ipotesi di tolleranza, quindi venga sempre fatto con un errore che rimane al di sotto di una certa soglia.

Segno con tolleranza.png

Assegno zero dove sono entro tolleranza ma comunque da il warning all'utente.esempio 0.001 e la tolleranza è 0.01. terremo questo qui come versione definitiva del ‘nostro se’!


  • Funzione blocco => (espressione1,espressione2.......ultimaespressione).

Valuta le espressioni in ordine restituendo come output solo la valutazione dell'ultima espressione. Esempio:

Funzioneblocco.png


  • Espansione in serie di Taylor:
f(t) :=cos(t);

Se al posto della variabile t si inserisse x il risultato non cambierebbe.

Esempio di sviluppo di f(x) in x nel'intorno di x=0 fino al grado 7:

Max9.JPG

tp7:taylor (f (x), x, 0, 7);
x0: %pi/3;
taylor (f (x), x, x0, 7);

dove f(x)= funzione,

    x= variabile,
    0= intorno,
    7= grado

Si può graficare la funzione di espansione in serie di Taylor fino a un certo ordine


  • Grafico della funzione f(x) ===> (funzione,[variabile, intervallo]).
 wxplot2d( funzione da rappresentare , [variabile , inizio intervallo , fine intervallo] );
 wxplot2d(f(x),[x,-2*%pi,2*%pi]);


N.B. Se non scrivo la sigla "wx" davanti a "plot" il grafico appare su una finestra a parte, mentre con wxplot2d l'output viene inserito nella stessa finestra del programma.

ATTENZIONE: I grafici plottati qui di seguito utilizzano per il polinomio di taylor tp7 un x0=0 , seguendo alla lettera questo wiki invece avremmo x0=%pi/3 come ultima definizione. Ovviamente taylor è sensibile al punto x0 che scegliamo. Come ci confermeranno i grafici mettere x0=0 significa che il tp7 approssimerà bene la nostra f(x) in un intorno di 0.


Max10.JPG


  • Grafico di 2 funzioni f(x) e tp7 a confronto
 wxplot2d( [funzione1 , funzione2] , [varibile , inizio intervallo , fine intervallo] );
 wxplot2d([f(x),tp7],[x,-2*%pi,2*%pi]);

Max11.JPG

Posso anche scrivere :

 wxplot2d([f(x),tp7],[x,-2*%pi,2*%pi],[y,-1.5,1.5]);

Il plot si fa solo di espressioni/funzioni che hanno una valutazione numerica, per la funzione f(ax), per esempio, non si può fare. Nel comando plot2d inserisco due parametri: il primo è la funzione da graficare (se voglio plottare due funzioni o più metto una lista di funzioni) , il secondo è una lista che contiene la variabile da far scorrere (ad esempio x) e l'intervallo entro cui farla scorrere. In pratica limito l’intervallo di y.

Max12.JPG

Possiamo aggiungere una legenda aggiungendo all'interno di wxplot2d:

[ legend, "stringa da associare alla funzione1" , "stringa da associare alla funzione2" ] 

In questo modo possiamo definire manualmente i nomi delle funzioni

wxplot2d([f(x),tp7],[x,-2*%pi,2*%pi],[y,-1.5,1.5],[legend,"fun","tp7"]);

In questo caso si è anche posto un limite all'ordinata y con il comando "[y,-1.5,1.5]".

E' bene ricordare che x può essere sostituita con una qualsiasi variabile (ad. esempio t), y invece è fisso, cioè non è una variabile ma una parola chiave per definire l'asse y. Inoltre per "plottare" una funzione è importante definire una legenda.

Max13.JPG


  • Utilizzo del ciclo for per calcolare le espansioni in serie di taylor di f(x) fino all'ordine n:

Implemento uno stoccaggio in forma di lista, pronta per essere plottata.

Inizializzo liste di stoccaggio sulla sola funzione:

plottami_fun:[f(x)];

deve contenere almeno una funzione da plottare

plottami_legenda:[legend,"fun"];

Qui calcolo e accodo a tali liste i polinomi di taylor (e le relative etichette):

N:7;
for i:0 thru N step 1 do disp(i);

la riga precedente vuol dire: per i che va da 0 a N con passo 1, esegui il comando "disp(i)" (che impone di mostrare i risultati a video)

Per mettere più elementi devo aprire un blocco:


      for i:0 thru N step 1 do ( 
      disp(i),
      tmp:taylor(f(x),x,x0,i),
      plottami_fun :append(plottami_fun,[tmp]),
      plottami_legenda :append(plottami_legenda,[i])
      );

N.B : qui x0 torna ad essere x0=%pi/3 , quindi se non lo si è già fatto va riassegnato

Il comando tmp indica l'operazione da ripetere, in questo caso "taylor(f(x),x,x0,i)".

Il comando append serve per concatenare due funzioni. In questo caso è stato utilizzato per allungare elementi a una lista, ed in particolare per espandere il polinomio fino al grado i.


 plottami_fun;

Max14.JPG

 plottami_legenda;


Max15.JPG

 wxplot2d(plottami_fun,[x,-2*%pi,2*%pi],[y,-1.5,1.5],plottami_legenda);


Max16.JPG


  • Parti reali ed immaginarie di numeri complessi:
      complesso : %i*sin(t)+cos(t);

Per prendere solo la parte reale:

      realpart(complesso);

Per prende la sola parte immaginaria:

      imagpart(complesso);

Per quanto riguarda l'integrazione numerica generalmente si usa il metodo dei trapezi (lineare) o di Simpsons (parabolico). Per applicare il metodo dei trapezi bisogna seguire i passaggi: (1) dividiamo la funzione in n intervalli (2) campioniamo n+1 punti sulla curva corrispondenti agli estremi degli intervalli (3) costruiamo dei trapezi che approssimano l'integrale della funzione In Maxima c'è già un algoritmo implementato che usa il metodo dei trapezi.

Integrazione numerica:

      quad_qag(sin(x)*x,x,0,%pi,1);

il primo termine è la funzione da integrare, il secondo la variabile d'integrazione, il terzo e il quarto sono gli estremi dell'intervallo d'integrazione e l’"1" finale indica il metodo dei trapezi (2 per il metodo di Simpson ecc...)


L'output si presenterà nel modo seguente:

       [risultato integrale, errore compiuto, numero di raffinamenti per arrivare a precisione, 0 se tutto è andato bene]

Max17.JPG

quad_qag (sin (x)*x, x, 0, %pi, 1);

Cicli sugli elementi di una lista

E' possibile creare rapidamente una lista attraverso il comando:

    lista:makelist(i^2,i,1,10);

questo fa variare la i da 1 a 10 creando così una lista del tipo:

    [1,4,9,16,25,36,49,64,81,100]

E' possibile anche stabilire un passo (es. 2) di variazione della i, indicandolo come ultimo termine:

    lista1:makelist(i^2,i,1,10,2);

con un output quindi del seguente tipo:

    [1,9,25,49,81]

Posso a questo punto sfruttare gli elementi della lista creata in un ciclo for, dopo aver inizializzato a 0 una variabile "poly":

    for i:1 thru length (lista) do poly:poly+x^lista[i];

Poly1.png

il comando lenght(lista) estrae la lunghezza della lista in esame, e il ciclo for costruisce un polinomio in cui la x è elevata a ciascuno dei termini della lista creata in precedenza. La lunghezza della lista è il numero degli elementi, per scorrere il numero degli elementi uso la variabile i.

Metodo equivalente a quello appena proposto è il seguente:

Poly.png

Matrice B per gli elementi triangolari a tre nodi

La matrice B determina il legame tra spostamenti e deformazioni. Gli spostamenti nodali possono essere definiti come:

u(x,y):= a1*x+a2*y+a3;  spostamenti orizzontali
v(x,y):=a4*x+a5*y+a6;   spostamenti verticali

Successivamente devo trovare un legame fra i coefficienti a1, a2,......, a6, cioè ricavo le costanti a1,..,a6 in funzione degli spostamenti campionati ai nodi. Creo una lista di 3 equazioni in tre incognite per definire lo spostamento orizzontale:

eqns_u:[u(x_i,y_i)=u_i,      
       u(x_j,y_j)=u_j,    
       u(x_k,y_k)=u_k
      ];
linsolve (eqns_u, [a1,a2,a3]);

Legame coef a.png

Alternativamente se si vuole visualizzare il sistema in forma matriciale si può usare il comando:

augcoefmatrix(eqns_u, [a1,a2,a3]); 

Questo comando mi da la matrice aumentata dei coefficienti cioè la matrice dei coefficienti più la colonna dei termini noti cambiata di segno, se voglio solo la matrice dei coefficienti tolgo aug.


NB: PER CONTINUARE A SCRIVERE IL PROGRAMMA CANCELLARE O COMMENTARE IL COMANDO AUGCOEFMATRIX!!!

Ora si vuole dare alle incognite a1,a2,a3 il valore precedentemente trovato:

linsolve (eqns_u, [a1,a2,a3]), globalsolve=true;

Di conseguenza si crea una lista di 3 equazioni in tre incognite per definire lo spostamento verticale:

eqns_v:[v(x_i,y_i)=v_i,      
       v(x_j,y_j)=v_j,    
       v(x_k,y_k)=v_k
      ];
linsolve (eqns_v, [a4,a5,a6]), globalsolve=true;

Per il calcolo delle deformazioni si procede con la differenziazione degli spostamenti: u rispetto ad x, v rispetto ad y ed infine con la derivata mista di u e v. Si passa ora all'estrazione delle deformazioni:

eps_x: diff(u(x,y),x);
eps_y: diff(v(x,y),y);
gamma_xy: diff(u(x,y),y)+ diff(v(x,y),x);

Spostamenti epsilon.png

Devo verificare inoltre che non ci sia dipendenza tra x e y, definendo una lista:

eps:[eps_x,eps_y,gamma_xy];

derivo rispetto a x e y e verifico che non ci sia dipedenza:

diff(eps,x);
diff(eps,y);


invece per verificare la dipendenza dagli spostamenti nodali procedo come segue

diff (eps,u_i);

lo eseguo per tutti gli spostamenti nodali.

Dipendenza spostamenti.png

Nell'espressione non sono presenti spostamenti nodali (u_i........v_k), per cui la relazione tra eps e spostamenti nodali è lineare, per cui è esprimibile come pre moltiplicazione per la matrice dei coefficienti.

Una volta raccolte le componenti degli spostamenti entro una lista/vettore chiamata delta, posso estrarre la matrice dei coeficienti che premoltiplicata per delta da eps.

delta:[u_i,v_i,u_j,v_j,u_k,v_k];
B: coefmatrix(eps,delta);

eps= B * delta


eps è un vettore di 3 componenti (eps_x, eps_y,gamma_xy),B (detto "operatore differenziale") e' una matrice 3x6, delta è un vettore a 6 componenti (u_i,v_i,u_j,v_j,u_k,v_k).


(Per motivi di spazio non è stato possibile introdurre l'output di maxima)

Ora estraggo il termine 2*A(duearea), non faccio la lista di liste ma uso il comando matrix:

duearea:
determinant(
matrix( 
       [1,1,1],
       [x_i,x_j,x_k],
       [y_i,y_j,y_k]
       ));

Duearea.png

Trovo la forma con il termine due area raccolto:

ratsubst(2*A, duearea,B);

In B trovo due area e sostituisco con 2*A

Ratsubst.png


Il calcolo dell'area del triangolo è estendibile al volume dei tetraedro. L'area del triangolo, infatti, è esprimibile come:

 determinant(
 matrix( 
       [1,1,1],
       [x_i,x_j,x_k],
       [y_i,y_j,y_k],
       ))* 1/2;

Aggiungendo una riga (perché ho anche la variabile z) e una colonna si può trovare il volume del tetraedro:

 determinant(
 matrix( 
       [1,1,1,1],
       [x_i,x_j,x_k,x_1],
       [y_i,y_j,y_k,y_1],
       [z_i,z_j,z_k,z_1],
       ))* 1/6;

Bibliografia

Utili per il teorema di Castigliano, poco noto, e in generale per la teoria della scienza delle costruzioni, base della progettazione assistita teorica, sono i testi:

  • Rif. 1 Michele Capurso, Scienza delle costruzioni, Pitagora editrice, Bologna 1971.

e il libro del prof. Strozzi che propone la teoria per la progettazione meccanica:

  • Rif. 2 Antonio Strozzi, Costruzioni di Macchine, Pitagora editrice.