NUOVO METODO DI INTEGRAZIONE NUMERICA DELL'ING. LAMBERTO BERTOLI - LA NUOVA INTEGRAZIONE NUMERICA PER GLI ELEMENTI FINITI

Ingegneria Strutturale

· LA NUOVA INTEGRAZIONE NUMERICA PER GLI ELEMENTI FINITI

· Gli elementi finiti e la loro integrazione

· LA NUOVA QUADRATURA DIRETTA

· L’integrazione di funzioni di 2 variabili reali

· L’integrazione delle funzioni di 3 variabili reali

· Conclusioni

I fondamenti scientifici del metodo proposto dall’Ing. Lamberto Bertoli, i cui contenuti fondamentali sono discussi in questo articolo, sono stati presentati dai docenti Carmelo Majorana, Stefano Odorizzi e Renato Vitaliani dell’Università di Padova al “Fourth seminar about method and variational methods” tenutosi a Plzen (Praga) tra il 12 e il 15 maggio 1981.

Successivamente la trattazione fu pubblicata in versione ridotta sulla rivista “Advanced in Engineering Software” (C. Majorana, S. Odorizzi, R. Vitaliani, “Shortened quadrature rules for finite elements”, Advanced in Engineering Software, 1982, Vol. 4, N.2) e, come selected paper, nel volume “Software in Engineering Problems”.

Le soluzioni del problema dell’integrazione numerica, a partire dai casi più semplici, oggetto della tesi di laurea, fino ai più sofisticati, che sono stati affrontati nelle ricerche degli ultimi anni, sono state pubblicate dall’Ing. Bertoli nel testo:

“Quadratura diretta degli integrali multipli”, 2006, Ed. Libreria Internazionale Cortina, Padova.

Gli elementi finiti e la loro integrazione

Il metodo degli elementi finiti pone il suo fondamento nel principio di minima energia che Jacques Bernoulli (1654-1705) intuì alla fine del ‘600 riflettendo sulla forma assunta da una fune inestensibile tesa alle sue estremità. Egli comprese che il baricentro del sistema in equilibrio doveva necessariamente situarsi alla minima altezza rispetto a qualunque altra configurazione congruente coi vincoli ai suoi estremi.

Da questo assunto, il grande studioso concluse pertanto che l’equilibrio della catena impone la sua minima energia potenziale rispetto a qualunque forma assumibile dalla medesima. Nel caso di una fune estensibile, oltre all’energia potenziale è presente l’energia elastica di deformazione, per cui la complessiva è ottenuta aggiungendo quest’ultima all’energia potenziale gravitazionale. Joseph Louis Lagrange, ritenuto il più grande matematico della seconda metà del XVIII secolo, comprese le potenzialità di questo principio, e sviluppò il calcolo variazionale già anticipato dalle ricerche pionieristiche di Bernoulli. L’energia di una fune sospesa dipende infatti dalla sua configurazione e quando raggiunge un minimo, ne consegue che per ogni cambiamento infinitesimo della sua curvatura, la variazione di energia deve essere nulla. La forma assunta dalla fune in equilibrio deve quindi necessariamente soddisfare la condizione che per qualunque spostamento virtuale dei suoi punti, l’energia deve rimanere costante. Il principio dei lavori virtuali che viene così frequentemente utilizzato nella Scienza delle Costruzioni, si fonda proprio sul principio di minima energia, e ne rappresenta una delle più importanti applicazioni.

Molti fisici matematici dopo le fruttuose ricerche di Lagrange proseguirono gli studi sul calcolo variazionale, la cui diffusione si è estesa dopo la nascita degli elaboratori elettronici. Il metodo degli elementi finiti si fonda su questo tipo di calcolo e deve il suo crescente interesse alla sua generalità, in quanto è in grado di affrontare anche i problemi in cui si rende improponibile il calcolo analitico, che al contrario deve limitarsi allo studio di un numero limitato di situazioni idealizzate.

Sostituendo ad un sistema continuo un insieme di elementi di dimensioni finite, il principio di minima energia mantiene la sua validità assicurando che la configurazione assunta da questo modello in condizioni di equilibrio riprodurrà la condizione di energia costante per qualunque spostamento virtuale. Attribuendo a ciascun elemento un numero limitato di gradi di libertà determinato da punti collocati sui suoi bordi, è possibile determinarne l’energia di deformazione sulla base degli spostamenti assunti da ciascuno dei suoi punti campione collimanti con quelli degli elementi contigui. La funzione che esprime lo spostamento di ogni singolo punto dell’elemento una volta noto quello del punto di riferimento, sarà costituita da polinomi in grado di annullarsi in tutti gli altri punti base e dalla derivazione di queste funzioni è possibile calcolare la deformazione assunta dal materiale in tutto l’elemento. Da quest’ultima si risale all’energia attraverso il modulo di elasticità del materiale, ed integrando su tutto l’elemento si può infine determinare la matrice di rigidezza relativamente a ciascun grado di libertà stabilito per ogni punto campione.

Il problema principale del metodo degli elementi finiti è proprio il calcolo di questa matrice ottenuta con l’integrazione numerica estesa a tutto l’intervallo delle tensioni e delle deformazioni indotte dagli spostamenti unitari di ciascuno dei suoi gradi di libertà. Questa computazione è infatti così impegnativa da assorbire normalmente oltre la metà del tempo di calcolo di un programma strutturale.

La procedura comunemente usata si basa sulla quadratura di Gauss, che rappresenta la più rapida integrazione teoricamente possibile per funzioni di una sola variabile reale. Tale metodo fu studiato dal suo autore unicamente per funzioni di un’unica variabile, perché a quel tempo non era ancora di attualità la loro quadratura in un insieme a più dimensioni.

Gli elementi finiti invece richiedono normalmente l’integrazione di funzioni di 2, 3 e 4 variabili, per cui la tecnica di Gauss ha dovuto essere generalizzata ricorrendo a quegli stessi polinomi di Legendre su cui si basa la quadratura del grande matematico.

Questo metodo tuttavia non affronta in modo diretto l’integrazione di funzioni dipendenti da più grandezze, e quindi richiede un numero di punti crescente esponenzialmente col numero di variabili in gioco. Per questo motivo l’integrazione di Gauss-Legendre utilizza un numero di punti antieconomico rispetto a quello rigorosamente necessario, producendo così un aumento del tempo di calcolo, soprattutto negli elementi a 3 dimensioni e in quelli in cui è presente la dimensione temporale.

LA NUOVA QUADRATURA DIRETTA

I notevoli tempi di calcolo di programmi sempre più sofisticati, hanno imposto la necessità di affrontare in modo diretto l’integrazione di funzioni definite negli elementi normalizzati di forma quadrata nel caso di 2 variabili e di forma cubica negli elementi a 3 dimensioni.

Tale quadratura è stata affrontata dallo scrivente per la prima volta nel 1977 nella sua tesi di laurea integrando polinomi fino al 7° grado per funzioni di 2, 3 e 4 variabili, e in seguito ampliata estendendone la ricerca a polinomi di grado superiore. Il nuovo metodo consente di ottenere un risparmio di tempo, a parità di precisione, del 25% per funzioni definite nel piano, del 50% per polinomi in 3 dimensioni, e del 77% per la quadratura di polinomi in 4 variabili.

L’integrazione di funzioni di 2 variabili reali

Essa è stata ottenuta a partire dall’espansione di Taylor delle funzioni di due variabili contenenti, oltre ai termini del tipo xⁿ e y^m, anche le componenti miste del tipo x^py^q. Tale sommatoria comprende la presenza di monomi antimetrici rispetto alla variabile x, del tipo x²ⁿ⁺¹, e alla variabile y della forma y^2m+1, i cui effetti si compensano nell’elemento normalizzato di forma quadrata collocando l’origine degli assi nel suo centro, gli assi coordinati paralleli ai lati del dominio, e disponendo i punti campione in modo simmetrico rispetto ai due assi. Anche i termini antimetrici x²ⁿ⁺¹ y^2m e x^2p y^2q+1 vengono annullati perché assumono valori opposti nei semipiani simmetrici rispetto agli assi coordinati.

Di conseguenza tutte le componenti dispari dell’espansione di Taylor vengono eliminate nei calcoli collocando i punti campione in posizioni simmetriche rispetto agli assi coordinati, in quanto ogni valore positivo calcolato nel punto base di un quadrante viene compensato da un uguale valore negativo dovuto al suo simmetrico situato o nel semipiano opposto, o nel quadrante antimetrico. Per esempio, il termine x³ y² produce valori positivi nel semipiano positivo delle x, ma altrettanti valori opposti nel semipiano negativo. Invece il termine x³ y³ produce valori uguali nel I e III quadrante che vengono eliminati dai risultati opposti ottenuti nel II e IV quadrante.

Solo i termini pari del tipo x^2m y²ⁿ producono valori identici in ogni quadrante, e perciò è sufficiente considerarne gli effetti indotti limitandoci al I quadrante individuato dai valori positivi delle ascisse e delle ordinate.

I punti campione su cui va calcolato il valore della funzione sono in grado di risolvere identicamente sia i termini x^2my²ⁿ che le componenti x²ⁿ y^2m, a condizione che essi siano disposti simmetricamente rispetto alla bisettrice del I e III quadrante, collocazione che riduce il numero di incognite e di equazioni risolutrici nelle quadrature di grado più elevato.

Si deve precisare tuttavia che la quadratura rapida dei polinomi di grado uguale o inferiore al quarto è stata possibile rinunciando alla simmetria rispetto alla bisettrice, perché solo in questo modo è stato limitato al minimo teorico il numero di punti campione. Questa circostanza ha tuttavia comportato un incremento delle difficoltà di calcolo, dovendo considerare un maggior numero di termini dell’espansione di Taylor al fine di una risoluzione corretta del problema.

Dal grado 6° in poi, la quadratura diretta è stata invece risolta unicamente con una disposizione simmetrica dei punti sia rispetto agli assi coordinati, sia rispetto alle loro bisettrici, rendendo così possibile l’eliminazione dei termini x^2m y²ⁿ una volta risolta l’integrazione dei termini x²ⁿ y^2m.

Con queste premesse, la quadratura di un qualunque polinomio di grado uguale o inferiore ad un numero prestabilito, richiede che i punti campione siano scelti in modo da integrare esattamente i termini x²ⁿ, y^2m e le componenti miste x^2py^2q fino al raggiungimento di un grado uguale a quello fissato a priori. Infatti, se i punti campione, di cui non si conosce ancora la posizione e i relativi coefficienti di peso, sono in grado di integrare esattamente i singoli monomi di grado pari dell’espansione di Taylor, allora saranno in grado di integrare qualsiasi combinazione lineare di questi ultimi, risolvendo così la quadratura di un qualunque polinomio di grado uguale, ma anche inferiore, a quello prestabilito. Ogni polinomio è infatti per definizione una somma di monomi, per cui è possibile generarne uno qualsiasi attribuendo ai coefficienti dei singoli termini i valori da noi scelti. La quadratura di un qualsiasi polinomio è quindi possibile se i punti campione vengono collocati in posizioni tali da essere in grado di integrare ciascuno dei suoi componenti di grado pari.

Ognuno dei termini indipendenti del polinomio produce un’equazione non lineare in cui al primo membro si calcola la sommatoria, estesa a tutti i punti, dei prodotti del valori assunti dal monomio integrando in ogni punto, moltiplicati per il proprio coefficiente di peso. Per esempio, il termine x² y⁶ dovrà essere calcolato in ciascuno dei punti base incogniti, e tali valori x_i² y_i⁶ saranno moltiplicati per i coefficienti di peso H_i relativi ai medesimi punti e ancora incogniti.

Al secondo membro si colloca il valore assunto dall’integrale del monomio x² y⁶ nell’intervallo delle x e delle y compreso fra zero e uno. Tale valore nel caso del monomio x² y⁶ è pari a 1/21 ed è sempre un valore frazionario per qualsiasi termine considerato.

Estendendo l’integrazione a tutte le componenti di grado pari, con l’esclusione di quelle simmetriche rispetto alla bisettrice del I quadrante, si ottiene con questa procedura un sistema di equazioni non lineari in cui le incognite sono rappresentate dalle coordinate (x_i,y_i) e dai coefficienti di peso H_i di ciascun punto.

Si deve poi precisare che nelle incognite del problema compaiono pure alcuni punti campione situati nelle bisettrici degli assi, le cui incognite si limitano alla sola coordinata x_i e al coefficiente di peso H_i, in quanto y_i coincide con x_i.

Il coefficiente di peso H_i ha un significato non solo algebrico, ma anche geometrico. Esso rappresenta infatti il volume sotteso dalla superficie generata da un polinomio in (x, y), avente la proprietà di assumere un valore unitario nel punto P_i e di annullarsi in tutti gli altri. Tale polinomio è di grado uguale o anche inferiore a quello del polinomio di grado massimo integrabile esattamente, ed è costituito da un insieme di termini di grado sia pari che dispari, che riportano a zero il valore assunto dal polinomio in tutti i punti base, tranne che in quello in cui deve essere calcolato il valore della funzione integranda.

Se ora vogliamo quadrare una generica funzione in due variabili di cui sono noti i valori assunti dalla medesima nei punti prestabiliti, il computo prevede che dobbiamo sommare i prodotti H_if(x_i,y_i) valutati in ogni punto P_i.

Questa operazione equivale alla sostituzione della funzione integranda col polinomio che assume nei punti base gli stessi valori della medesima, attribuendo così all’integrale ignoto della funzione data, il corrispondente integrale del polinomio che intercetta i medesimi valori della funzione nei punti campione.

E’ chiaro che quanto più elevato sarà il numero di punti di riferimento, quanto più il polinomio interpolante si confonderà con la funzione integranda anche nel contorno dei medesimi punti, per cui col crescere del numero di punti base, si incrementa pure la precisione della quadratura. La differenza fra il valore dell’integrale esatto e di quello così calcolato viene chiamata resto, e dipende dalle derivate parziali (n+1)-esime della funzione integranda calcolate in un punto interno al dominio di definizione.

Come già precisato, il polinomio interpolante è a sua volta la somma dei polinomi ottenuti moltiplicando per f(x_i,y_i) ciascuno dei polinomi che assumono un valore nullo in tutti i punti, tranne che in P_i, dove valgono 1. Questi nuovi polinomi ottenuti moltiplicando per f(x_i,y_i) quelli generanti con la loro cubatura il coefficiente di peso H_i, descrivono una superficie che ha la proprietà di intercettare nel punto P_i lo stesso valore della funzione integranda, mentre interseca il piano (x, y) in tutti gli altri punti base, dove assume un valore nullo. In questo modo, con la sommatoria dei termini H_i f(x_i,y_i), non si fa altro che calcolare il volume del polinomio che riproduce esattamente nei punti prestabiliti tutti i valori assunti dalla funzione integranda.

La quadratura diretta delle funzioni di due variabili non è altro che una generalizzazione nel piano (x,y) del metodo di Gauss, in quanto affronta la risoluzione di un sistema di equazioni non lineari in cui sia i coefficienti di peso, sia le coordinate dei punti campione costituiscono le incognite del problema.

La equazioni non lineari descritte dall’integrazione gaussiana vengono riprese nella quadratura diretta, ampliando tuttavia sia la loro quantità, sia il numero di incognite, che ora non si limitano più alle posizioni x_i dei punti base di una retta, ma anche alle coordinate y_i dei punti del piano di cui devono essere stabilite entrambe le coordinate (x_i,y_i ). Questo comporta che a parità di grado, il numero di incognite e quindi di equazioni, è superiore rispetto a quanto si presenti nella trattazione di una sola variabile, anche se l’incremento non è così rapido come con l’utilizzo dei polinomi di Legendre.

Si presenta inoltre una nuova complicazione che ha ritardato di molto la soluzione del problema che altrimenti avrebbe potuto essere già risolto fin dagli esordi del metodo degli elementi finiti. Essa riguarda la condizione tassativa di accettare unicamente quelle soluzioni che riproducono unicamente valori positivi dei coefficienti di peso H_i. Valori negativi sono indice di una scelta infelice delle coordinate dei punti base e rendono instabile la soluzione creando problemi nei programmi di calcolo.

Nel metodo di Gauss, i coefficienti di peso sono per loro natura tutti positivi, e questo è reso implicito dalle proprietà intrinseche dei polinomi di Legendre su cui si fonda la quadratura. Ne consegue che una volta impostate le equazioni costitutive, qualunque sia il loro numero, i sistemi di equazioni che ne derivano hanno sempre soluzioni positive che possono essere calcolate mediante un’equazione algebrica associata di grado dipendente sia dal numero di punti base situati nel semiasse positivo delle ascisse che eventualmente da quello collocato nella loro origine.

I metodi dell’analisi numerica non hanno perciò difficoltà a risolvere il sistema di equazioni ottenuto dalla quadratura di Gauss, perché il problema viene trasformato nella risoluzione di un’equazione algebrica di una sola variabile che restituisce tante soluzioni positive quante il numero di punti campione prestabilito.

Diverso è il caso della risoluzione delle equazioni relative all’espansione di Taylor di una funzione in due variabili. Queste equazioni non generano automaticamente valori positivi dei coefficienti di peso e basta anche la presenza di un solo H_i negativo per invalidare l’intera quadratura. Questo problema si aggrava con l’aumento del grado del polinomio interpolante e diventa di ardua difficoltà a partire dai polinomi di 12° grado per funzioni di due variabili e già di 8° grado per polinomi in 3 variabili.

La causa di questo fatto è dovuta alle condizioni sempre più sottili che devono essere rispettate, ma soprattutto individuate, richieste nella collocazione dei punti base del I quadrante, allo scopo di originare in seguito sequenze di risultati sempre positivi.

Questi punti interferiscono in tutte le equazioni risolutrici e sono presenti non solo nell’integrazione dei termini x^2my²ⁿ, ma anche nella quadratura dei monomi del tipo x^2p e y^2q che si incontrano nello sviluppo in serie di Taylor.

Le funzioni del piano governate da una sola variabile come x^2p o y^2q sono delle superfici rigate, ottenute quindi da un insieme di infinite rette parallele e perpendicolari ai piani coordinati, intersecanti la parabola di grado corrispondente definita nel piano (x, z) o (y, z). Queste superfici generano una cubatura nel piano (x, y), e il nostro interesse riguarda l’integrale definito nel dominio di definizione compreso tra 0 e 1 degli assi x e y. E’ chiaro che esse assumono in ogni punto del piano e degli assi dei valori ben precisi che sono calcolabili e successivamente moltiplicabili per i coefficienti di peso dei punti base. I punti campione collocati sugli assi cartesiani ignorano i termini del tipo x²ⁿy^2p e quindi non sono presenti nelle equazioni relative alla loro integrazione. I punti situati sugli assi sono invece rilevanti unicamente nell’integrazione dei termini del tipo x^2pe y^2q e sono necessari per la quadratura di un polinomio in più variabili.

Se i punti P_i(x_i y_i) non sono collocati nel numero e nella posizione più adeguata, saranno proprio i coefficienti H_0i dei punti situati sugli assi coordinati a decadere più facilmente in quei valori negativi che invalidano l’intera soluzione, evidenziando così la criticità delle scelte di base. Poiché il loro valore viene calcolato in modo sequenziale solo dopo la risoluzione delle relazioni dipendenti unicamente dai punti P_i(x_i,y_i), sono proprio quelle contenenti solamente i termini x^2my²ⁿ a determinare la riuscita delle equazioni in cui compaiono anche i punti appartenenti agli assi cartesiani.

L’integrazione delle funzioni di 3 variabili reali

In questo tipo di integrazione sono presenti nuovi termini del tipo x^2my²ⁿz^2p che non esistevano nei polinomi in due variabili e che ora richiedono l’introduzione di nuovi punti campione nello spazio (x,y,z). Ancora una volta il problema può essere semplificato ricorrendo a simmetrie che consentono di isolare le equazioni nel diedro formato dai semipiani positivi delle variabili x, y e z.

Mentre i problemi più sofisticati nel piano vengono risolti da coppie di punti con coordinate simmetriche rispetto alla bisettrice del primo quadrante del tipo P_iA(x_i,y_i) e P_iB(y_i,x_i), nell’integrazione dei polinomi più complessi in 3 variabili, si deve soddisfare una condizione di triplice simmetria con insiemi di punti del tipo P_iA(x_i,y_i,y_i), P_iB(y_i,x_i,y_i) e P_iC(y_i,y_i,x_i). Queste terne sono in grado di integrare correttamente tutti i termini del tipo x^2my²ⁿz^2p quando sono associate ad altri punti base situati sulla diagonale che congiunge l’origine degli assi col vertice di coordinate x, y e z pari a (1,1,1). In questo modo il problema si riconduce all’integrazione dei polinomi in 2 variabili, anche se ora il numero di punti da considerare è notevolmente superiore, e le equazioni risolutrici dipenderanno da un più elevato numero di incognite. Disporre di una terna di punti nello spazio equivale a collocare 2 punti simmetrici rispetto alla bisettrice del piano (x,y), oltre ad un terzo punto sulla medesima retta.

La quadratura prevede la contemporanea risoluzione delle equazioni generate dai termini del tipo x^2my²ⁿz^2p e delle relazioni generate dalla proiezione nel piano (x,y) dei medesimi punti base nello spazio (x,y,z). Pure nel problema spaziale devono infatti essere risolte le equazioni in 2 variabili generate dai termini x^2my²ⁿ, e questa soluzione può essere ottenuta ricorrendo alle proiezioni sul piano (x,y) sia delle terne P_iA, P_iB e P_iC, che dei punti situati sulla diagonale del cubo che pure vengono proiettati sulla bisettrice del quadrato. La soluzione è resa possibile dall’aggiunta di un certo numero di punti sulla medesima bisettrice del piano (x,y) che consentono così di integrare esattamente anche i termini del tipo x^2my²ⁿ. Infine devono ancora una volta essere disposti altri punti sugli assi coordinati, al fine di poter integrare le funzioni del tipo x²ⁿ, y^2m, e z^2p.

Anche nello spazio la principale difficoltà consiste nel posizionare nel modo più opportuno la terna di punti simmetrici rispetto ai tre piani cartesiani, il cui scopo non è solo quello di consentire l’integrazione dei termini x^2my²ⁿz^2p, ma anche dei monomi x^2my²ⁿ definiti nel piano (x,y). Solo una disposizione corretta sarà in grado di produrre unicamente valori positivi dei coefficienti di peso e quindi di rendere accettabile la soluzione.

La validità della configurazione scelta dipenderà dall’individuazione delle sottili e non sempre evidenti condizioni che consentono di indurre una cascata di valori tutti positivi dei coefficienti di peso H_i, sia nei punti situati nel piano (x,y), che sugli assi coordinati. Per questo motivo l’integrazione di un polinomio in 3 variabili è, a parità di grado, più impegnativa della quadratura delle funzioni in 2 variabili, ma la maggior ostilità dei calcoli viene compensata da un più generoso risparmio di punti rispetto alle funzioni di 2 variabili.

Conclusioni

In questo articolo sono stati presentati alcuni dei principi costitutivi del nuovo metodo di integrazione che consente, a parità di precisione, un risparmio di punti del 25% per funzioni di 2 variabili, del 50% per polinomi di 3 variabili, e del 77% per funzioni di variabili.

Come già anticipato, il procedimento qui presentato non è altro che una generalizzazione del metodo che Gauss studiò per l’integrazione numerica delle funzioni di variabile reale. Le equazioni risolutrici sono quindi simili a quelle determinate dalla quadratura gaussiana, cui si aggiungono nuove relazioni che comprendono termini misti del tipo x^2my²ⁿ per le funzioni di 2 variabili, mentre nei polinomi in 3 dimensioni sono stati pure considerati i monomi della forma x^2my²ⁿz^2p.

Il metodo, le cui origini risalgono al 1977, è stato implementato con successo nel Centro di Calcolo del Dipartimento di Costruzioni dell’Università di Padova dal Prof. Renato Vitaliani, Ordinario di Tecnica delle Costruzioni presso la Facoltà di Ingegneria. La diffusione internazionale della nuova quadratura ha convinto lo scrivente a riprendere le ricerche e ad affrontare l’integrazione di polinomi di grado superiore al 7° che non erano ancora stati trattati nella tesi di laurea.

Gli studi avviati dall’autunno del 2000, hanno prodotto nuove conclusioni e sono state risolte positivamente le quadrature di polinomi in 2 dimensioni fino al 15° grado e di 9° grado per funzioni di 3 variabili. Gli esiti ottenuti a partire dalle prime ricerche sono stati pubblicati dallo scrivente nel testo edito dalla Libreria Internazionale Cortina di Padova: “Quadratura diretta degli integrali multipli”, con ISBN 88-7784-263-6.

L’autore del presente articolo si propone di approfondire lo studio di nuovi casi che consentano di risolvere non solo le problematiche relative al piano e allo spazio, ma anche alla dimensione spazio temporale, incrementandone così la precisione finora raggiunta.