· Meccanica Computazionale

· Nuovo Metodo di Integrazione Numerica dell'Ing. Lamberto Bertoli

· UN NUOVO METODO DI INTEGRAZIONE NUMERICA PER GLI ELEMENTI FINITI

· Lamberto Bertoli

· INTRODUZIONE

· I principi fondanti la quadratura di Newton-Cotes

· L’integrazione di Gauss-Legendre

· UN NUOVO METODO PER LA quadratura diretta degli integrali multipli

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.

INTRODUZIONE

In questo articolo viene presentato un nuovo metodo di integrazione degli integrali multipli, proposto dallo scrivente, per la risoluzione dei problemi riguardanti gli elementi finiti. Tale quadratura costituisce una generalizzazione della procedura di Gauss che, conservandone la precisione, consente di aumentarne la velocità di calcolo e ridurre del 25% il numero di punti campione per gli integrali doppi, del 50% per gli integrali tripli e del 77% per gli integrali a 4 variabili.

Esso è nato dalla risoluzione dei sistemi di equazioni non lineari generati dall’espansione di Taylor delle funzioni di più variabili, che a causa del loro grado estremamente elevato non possono essere affrontati coi tradizionali metodi dell’analisi numerica. Il malcondizionamento di questi sistemi è stato superato attraverso originali e innovative tecniche di calcolo, che hanno consentito di individuare il minimo numero teorico di punti necessario per la quadratura numerica degli integrali multipli. I coefficienti di peso così calcolati sono tutti positivi, conformemente all’integrazione di Gauss di cui rappresenta una estensione a più variabili.

I principi fondanti la quadratura di Newton-Cotes

L’integrazione numerica di una funzione di variabile reale ha avuto le sue prime origini nel diciassettesimo secolo ad opera di Bonaventura Cavalieri che riuscì a individuare una procedura in grado di calcolare l’area sottesa da una parabola passante per tre punti dati. Tale formula, che ha assunto il nome di Cavalieri-Simpson, è stata in seguito generalizzata alla fine del ‘600 da Newton e Cotes che stabilirono un metodo generale per integrare un polinomio di qualunque grado passante per un prefissato numero di punti fra loro equidistanti.

Tale procedimento richiede tuttavia un numero di punti campione piuttosto elevato per la quadratura delle funzioni più sofisticate. Inoltre, a partire da un certo grado di precisione, comporta l’introduzione di alcuni coefficienti di peso negativi, che dal punto di vista della stabilità della risoluzione, può comportare qualche criticità.

Il metodo di Newton-Cotes, utilizzando punti distribuiti uniformemente all’interno del dominio di definizione, genera funzioni interpolanti predefinite che non costituiscono un’incognita del problema. Tali funzioni sono polinomi di Lagrange aventi la proprietà di annullarsi in tutti i punti dell’intervallo, tranne che in quello in cui va valutato il valore della funzione integranda. Esse sono facilmente calcolabili partendo da una produttoria di binomi del tipo (x – x_j), dove x_j rappresenta la coordinata di uno dei rimanenti punti dell’intervallo in cui il polinomio interpolante si annulla. Il risultato così ottenuto va quindi diviso per il numero ottenuto dal prodotto dei fattori (x_k – x_j)

Il coefficiente di peso H_k rappresenta l’area del polinomio di Lagrange avente la proprietà di annullarsi in tutti i punti base dell’intervallo, con esclusione di x_k in cui deve essere testato il valore della funzione integranda, dove assume un valore unitario

Moltiplicando la grandezza del coefficiente di peso H_k per la quantità assunta dalla funzione integranda nel punto x_k, si ottiene così l’area sottesa dal polinomio di Lagrange che assume gli stessi valori della funzione data nel punto x_k, mentre si annulla in tutti gli altri. Se ora sommiamo i risultati ottenuti in ciascun punto campione, si ottiene l’area del polinomio interpolante avente la proprietà di assumere gli stessi valori della funzione integranda in tutti i punti base in cui è stato suddiviso l’intervallo.

L’integrazione numerica di una qualsiasi funzione viene così semplicemente ottenuta moltiplicando i valori assunti dalla medesima nei punti prescelti per i relativi coefficienti di peso e sommando infine tutti i termini dati da questi prodotti. E’ intuitivo che la quadratura numerica sarà tanto più precisa, quanto più alto sarà il numero di punti considerati. Dati n punti campione, l’integrazione di Newton-Cotes determina infatti l’area sottesa dal polinomio interpolante di grado (n-1) passante per gli stessi punti intercettati dalla funzione integranda.

La differenza fra il valore dell’integrale così calcolato e quello della funzione integranda dipenderà dall’estensione dell’intervallo e dal valore assunto in un particolare punto all’interno del medesimo dalla derivata n-esima della funzione data. Il metodo di Newton-Cotes ha avuto un successo incontrastato per secoli, in quanto consente di operare su numeri piuttosto semplici che non determinano calcoli troppo impegnativi nella risoluzione manuale del problema.

L’integrazione di Gauss-Legendre

In seguito Gauss riprese l’integrazione di Newton-Cotes rendendola più precisa a parità di numero di punti e quindi di tempo di calcolo. L’idea del grande matematico consisteva nell’individuare all’interno dell’intervallo di definizione le particolari posizioni dei punti base in grado produrre un valore nullo dell’integrale di qualunque polinomio di grado compreso fra n e (2n-1) passanti per i punti dati.

In questo modo si otteneva con soli n punti l’integrazione esatta di un polinomio di grado (2n-1), mentre con Newton-Cotes lo stesso numero di punti è in grado di integrare esattamente solo un polinomio di grado (n – 1). Per ottenere questo risultato Gauss doveva innanzitutto imporre che il polinomio intercettante tutti i punti base sottendesse un’area nulla. Questo comportava che la posizione dei punti dati non fosse più prestabilita a priori, ma che soddisfasse un’equazione a più incognite costituite appunto dalle coordinate dei punti stessi. Il metodo di Gauss utilizza polinomi di Legendre che sono derivati dai polinomi di Lagrange e contenenti ancora una produttoria di binomi del tipo (x – x_j), dove x_j rappresenta la coordinata di un generico punto campione. Sviluppando tale prodotto, si ottiene un polinomio di grado n in cui i vari coefficienti sono funzioni delle coordinate dei punti prescelti. La condizione di annullamento dell’integrale comporta un’equazione che non ha una sola soluzione, perché il suo numero di incognite corrisponde al numero di punti considerati.

Di conseguenza l’equazione che soddisfa il requisito di sottendere un’area nulla viene soddisfatta scegliendo arbitrariamente (n-1) punti a piacere mentre rimane univocamente determinata solo la coordinata x_i dell’ultimo punto. I rimanenti (n – 1) gradi di libertà ci consentono di imporre nuove condizioni al problema che migliorano ancora di più la precisione dell’integrazione.

La seconda equazione viene individuata imponendo che qualunque polinomio di grado (n + 1) passante per i punti dati abbia pure un’integrale nullo all’interno dell’intervallo. Tale polinomio è vincolato dalla condizione di intercettare gli n punti campione e quindi è necessariamente generato dalla medesima produttoria dei binomi (x – x_j), moltiplicata questa volta per x. In questo modo il grado del prodotto aumenta da n a (n + 1) e il suo sviluppo determina ancora un polinomio i cui termini dipendono dalla posizione di punti prescelti. Imponendo nuovamente che l’integrale del polinomio abbia un valore nullo all’interno del dominio di definizione, si ottiene una seconda equazione che forma sistema con la precedente.

Con la stessa procedura si ottengono una di seguito all’altra tutte le altre equazioni moltiplicando il prodotto dei binomi per le potenze di x fino al grado (n – 1). Poiché la semplice produttoria dei binomi ha grado n, moltiplicandola per x elevato a (n – 1), si ottiene per l’ultimo polinomio il grado (2n – 1). Il sistema sarà quindi costituito da n equazioni, in cui la prima è ottenuta moltiplicando il prodotto dei binomi per x° e l’ultima per x ^n-1

L’imposizione dell’annullamento degli integrali dei polinomi di Lagrange moltiplicati per una potenza di x, produce dei vincoli sulle coordinate dei punti campione che, una volta soddisfatti, generano i polinomi Legendre, coi quali opera l’integrazione di Gauss.

Le n equazioni sono formate da n incognite che rappresentano le coordinate dei punti campione. Non si tratta di un sistema lineare, perché le coordinate dei punti base dopo lo sviluppo dei prodotti vengono moltiplicate fra loro determinando così fattori via via più complessi man mano che aumenta il numero di punti in cui va testata la funzione. Nell’integrazione di Gauss, quindi, anche la posizione dei punti dell’intervallo non è nota, ma costituisce un’incognita del problema e va ricercata mediante la soluzione di un sistema non lineare di equazioni. Una volta risolto il sistema e calcolata la posizione dei punti, il valore dei coefficienti di peso può essere più facilmente determinato con le stesse procedure dell’integrazione di Newton-Cotes.

Un primo metodo consiste nello sviluppo del polinomio di Lagrange avente la proprietà di annullarsi in tutti i punti così calcolati, tranne che nel punto in cui si deve calcolare il coefficiente di peso. Poiché esso assume un valore unitario nel punto x_k in cui si deve testare la funzione integranda, si determina la funzione polinomiale che, sviluppata, ci consente di calcolare l’integrale definito che costituisce il coefficiente di peso H_k. Ripetendo l’operazione solamente per i punti situati sul semiasse positivo delle ascisse, si ottiene la serie dei coefficienti di peso che in un intervallo di n punti ci consente di integrare esattamente qualunque polinomio di grado (2n – 1) passante per i medesimi punti.

Una seconda procedura in grado di fornirci il valore dei coefficienti di peso nasce dalla considerazione che qualunque polinomio passante per gli n punti è costituito dalla somma di termini di grado uguale o inferiore a quello del polinomio interpolante. Affinché la quadratura sia di validità generale, è necessario che essa sia in grado di integrare qualsiasi componente polinomiale di grado inferiore o uguale a quella richiesta. Di conseguenza essa deve essere in grado di integrare qualsiasi parabola di grado inferiore o uguale a (2n-1).

I termini di grado dispari contenuti nei polinomi integrandi non producono alcuna conseguenza, perché i loro effetti sono antimetrici sul semiasse positivo delle x e sul semiasse negativo. Infatti, disponendo in modo simmetrico i punti campione sul semiasse positivo e negativo, e attribuendo lo stesso coefficiente di peso ai punti equidistanti dall’origine degli assi, tutti i termini dispari producono automaticamente integrali nulli perché ciò che incrementano sull’uno viene eliminato dagli effetti contrari prodotti sul semiasse opposto. Proprio per questo motivo l’integrazione di Gauss viene affrontata considerando unicamente le parabole di grado pari appartenenti al polinomio integrando, in quanto la presenza dei punti disposti simmetricamente sul semiasse negativo delle x ci consente di eliminare qualsiasi influenza dovuta ai termini di grado dispari dell’espansione di Taylor della funzione integranda.

Detto questo, il calcolo dei coefficienti di peso può essere effettuato considerando una serie di equazioni, in cui la prima nasce dall’integrazione di una funzione costante, cioè da un polinomio di grado zero. Questo comporta che la somma di tutti i coefficienti di peso compresi fra l’intervallo (-1) e (1) deve essere uguale a 2. Questa uguaglianza produce la prima equazione di un sistema in cui le incognite sono appunto i coefficienti di peso H_k.

La seconda equazione nasce dall’integrazione esatta di una parabola di secondo grado ed è quindi necessario che la somma dei coefficienti di peso moltiplicati per il valore assunto dalla medesima nei punti campione sia uguale a 2/3, che rappresenta l’area intercettata dalla parabola x² definita in un intervallo compreso tra (–1) e (1). Ripetendo l’operazione per una parabola di quarto grado, si trova che la terza equazione deve fornire come risultato 2/5.

Il metodo di Gauss, oltre ad essere più preciso di quello di Newton-Cotes, presenta il vantaggio che i suoi coefficienti di peso sono tutti positivi per qualsiasi numero di punti, e questo vantaggio conferma che esso genera una quadratura non solo più rapida, ma anche più affidabile. I coefficienti di peso positivi sono infatti preferiti dai tecnici perchè non fanno insorgere instabilità o criticità nello svolgimento dei calcoli.

Prima dell’avvento del calcolo elettronico, la quadratura di Gauss non era tuttavia molto praticata perché opera quasi esclusivamente con numeri irrazionali che rappresentano le coordinate dei punti campione e dei coefficienti di peso. Questo inconveniente non ha attualmente alcuna rilevanza e di conseguenza il metodo di Gauss viene oggi preferito a quello di Newton-Cotes per la sua precisione e sicurezza di calcolo.

Tuttavia, il metodo di Gauss, per quanto oggi vincente, è stato studiato dal suo autore unicamente per l’integrazione di funzioni di una sola variabile reale e quindi non può essere utilizzato senza adattamenti agli elementi finiti che invece ricorrono quasi sempre a funzioni di più variabili. Per estenderlo a questi elementi è stato perciò necessario generalizzarlo mediante i polinomi di Legendre, le cui proprietà hanno reso così vantaggiosa la quadratura di Gauss a una sola dimensione. Come abbiamo visto, questi polinomi non sono altro che particolari polinomi di Lagrange ottenuti con una specifica collocazione dei punti campione che consente di migliorare la precisione della quadratura. Attraverso questi ultimi, è possibile costruire un polinomio di 2 variabili semplicemente moltiplicando le produttorie nella variabile x per altrettanti prodotti nella variabile y.

Il polinomio di Legendre in due variabili così ottenuto si annulla in una rete le cui maglie si congiungono nei punti base del quadrato rendendo possibile l’integrazione numerica di funzioni di due variabili con gli stessi principi dell’integrazione di funzioni di una sola variabile.

Il valore dei coefficienti di peso dell’integrazione a due variabili si ottiene semplicemente dall’integrazione di Gauss di una funzione di variabile reale moltiplicando il coefficiente H_k relativo all’ascissa x per il peso Hj relativo all’ordinata y. Naturalmente la procedura può essere generalizzata per qualunque numero di variabili, al prezzo però di aumentare in modo esponenziale il numero di punti campione.

Questo comporta un notevole incremento del tempo di calcolo con l’aumentare del numero di variabili, al punto che l’integrazione numerica determina da sola oltre la metà del tempo impiegato dal metodo degli elementi finiti.

Un’alternativa all’integrazione di Gauss-Legendre è stata formulata nel ventesimo secolo dal matematico John Von Neumann in collaborazione con Ulam. Tale metodo, chiamato di Monte Carlo è una procedura stocastica che tuttavia viene applicata raramente agli elementi finiti perché richiede abilità specialistiche e solo in pochi casi si dimostra competitivo con la quadratura di Gauss-Legendre.

Quest’ultima integrazione conserva quindi un ruolo fondamentale nel metodo degli elementi finiti, anche se essa rappresenta una quadratura indiretta e proprio per questo motivo richiede un numero superfluo di punti campione. Allo scopo di risolvere il problema ricorrendo al numero minimo di punti teoricamente possibile, è invece necessario generalizzare il metodo di Gauss estendendolo a funzioni di più variabili con l’integrazione diretta, che tuttavia solo recentemente è stata risolta grazie anche ai progressi delle tecnologie informatiche.

UN NUOVO METODO PER LA quadratura diretta degli integrali multipli

Lo scrivente ha iniziato a trattare il problema nel 1976 quando il Prof. Renato Vitaliani gli propose di affrontare la quadratura diretta degli integrali multipli come tema per la tesi di laurea. La relazione fu discussa col compianto Prof. Ubaldo Richard, matematico di chiara fama e allora Direttore dell’Istituto di Matematica Applicata della Facoltà di Ingegneria dell’Università di Padova. Il Prof. Richard accettò di seguire la tesi grazie anche alla richiesta del Prof. Lorenzo Contri, suo cordiale amico, che fin dall’inizio aveva creduto nella validità dello studio suggerito dal Prof. Vitaliani.

La ricerca fu portata a termine in tempi rapidi, ma lo scrivente dovette limitarsi a trattare la quadratura degli integrali multipli fino al 7° grado di precisione del polinomio interpolante, date le ardue difficoltà insorgenti nella risoluzione dei gradi più elevati.

Il metodo generalizzava l’integrazione di Gauss considerando l’espansione di Taylor delle funzioni di 2, 3 e 4 variabili che abitualmente si incontrano negli elementi finiti.

Anche le funzioni di 5, 6 e 7 variabili venivano considerate, ma la risoluzione si limitava al 5° grado di approssimazione, dal momento che tali casi vengono incontrati solo in rari problemi specialistici, come ad esempio l’esplosione delle supernove. La quadratura venne testata nel Centro di Calcolo dell’Istituto di Costruzioni Ponti e Strade e si dimostrò così efficace da venire implementata al posto della tradizionale integrazione di Gauss-Legendre.

Nel 1981 il Prof. Vitaliani, assieme a Stefano Odorizzi e a Carmelo Majorana, presentò al Congresso Internazionale di Calcolo Variazionale tenutosi a Plzen (Praga) il nuovo metodo di calcolo, corredandolo di esempi applicativi al metodo degli elementi finiti.

Le ricerche di questi autori dimostravano che la nuova quadratura consentiva di riprodurre la medesima precisione del metodo di Gauss-Legendre, ma con un risparmio di tempo di calcolo del 25% per funzioni di 2 variabili, del 50 % per le funzioni di 3 variabili e del 77 % per polinomi di 4 variabili.

Successivamente, in versione ridotta, questi studi vennero pubblicati nella rivista “Advances in Engineering Software” e, come selected papers, nel volume “Software in Engineering Problems”. La diffusione internazionale incontrata dal nuovo metodo, che attualmente è preferito a quelli tradizionali in diversi Paesi del mondo, convinse lo scrivente a riprendere le ricerche sull’argomento.

Fu così che nell’autunno del 2000 venne affrontata dallo scrivente l’integrazione di un polinomio di 9° grado in due variabili, problema questo già studiato senza successo nel 1976 al tempo della tesi di laurea. La difficoltà nasceva dal fatto che il sistema di equazioni era troppo complesso per la risoluzione algebrica già utilizzata per i casi di grado meno elevato e in un primo tempo fu perciò tentata la soluzione con le classiche procedure dell’analisi numerica.

Il metodo di Newton-Raphson, generalizzato a sistemi di equazioni non lineari, si dimostrò del tutto inadeguato per risolvere la questione, perché la sua efficacia dipende dalla posizione iniziale dei punti in cui si ritiene possa trovarsi la soluzione.

Tanto più alto è il grado del sistema, tanto più prossimi alla soluzione reale devono essere collocati i punti di partenza, perché altrimenti si innesca un movimento a serpentina che facilmente determina punti di stallo in grado di arrestare lo sviluppo dei calcoli.

Era quindi necessario individuare metodi alternativi per aggirare il malcondizionamento delle equazioni costitutive e per questo motivo lo scrivente ricorse a sistemi di equazioni integrali che consentivano di ridurre il numero di incognite del problema.

Dopo alcuni mesi, nel febbraio del 2001, fu individuata la soluzione con 21 punti dell’integrazione di polinomi di grado 9 a due variabili, trovando dei coefficienti di peso tutti positivi. Col perfezionamento delle tecniche risolutive, fu in seguito individuata la soluzione che con 20 punti corrispondeva al numero minimo necessario per risolvere il problema al posto dei 25 richiesti dalla quadratura di Gauss-Legendre. Si rendeva così possibile trattare il caso riguardante la quadratura esatta di un polinomio di 11° grado e l’esperienza già acquisita consentì di risolvere con le stesse procedure di calcolo il nuovo problema con 25 punti al posto dei 36 della quadratura di Gauss-Legendre.

La serie fortunata non consentiva tuttavia di proseguire al caso successivo, perché col grado 13 per la prima volta si incontrarono dei coefficienti di peso negativi. Per questo motivo fu necessario modificare la posizione dei punti allo scopo di ottenere soluzioni veramente soddisfacenti. Fu così che il problema fu affrontato a più riprese, ottenendo finalmente soluzioni che risolsero le equazioni con coefficienti di peso interamente positivi.

L’integrazione diretta di funzioni di due variabili si concluse in seguito con la risoluzione della quadratura esatta di polinomi fino al 15° grado, la cui precisione è tale consentire il trattamento di qualsiasi problema riguardante gli elementi finiti a due variabili.

Rimaneva ancora da studiare lo spazio, che in precedenza era stato risolto fino al 7° grado nella tesi di laurea. Si deve aggiungere che in precedenza il passaggio da funzioni di due variabili ad altre di tre variabili non aveva comportato difficoltà insormontabili.

Purtroppo, aumentando il grado, cresce inevitabilmente anche il numero di punti campione, per cui si rende necessario individuare a priori le disposizioni sempre più complesse e varie dei punti base. Più alto è il numero di punti, maggiore è il numero delle possibili alternative, e si deve indovinare quale fra è queste in grado di generare coefficienti di peso positivi. Naturalmente non è possibile conoscere a propri quale fra queste collocazioni comporterà pesi positivi, per cui la risoluzione inevitabilmente richiederà un numero di tentativi tanto maggiore quanto più elevato è il grado del sistema e il numero di variabili.

La soluzione prevede infatti un certo numero di incognite costituite dalle coordinate spaziali e dai coefficienti di peso. I calcoli non producono istantaneamente tutte le soluzioni, ma una sequenza di risultati fra loro interdipendenti. Le prime incognite che vengono fornite sono proprio le coordinate x, y e z dei punti campione. Successivamente in ordine vengono dati uno di seguito all’altro i coefficienti di peso, che all’inizio sono sempre positivi. La serie dei risultati ottenuti non consente di valutare la riuscita della soluzione fino a quando non viene calcolato sulla scorta degli altri anche l’ultimo coefficiente di peso. Solo se pure questo è positivo, l’intera soluzione è accettabile, perché anche un minimo peso negativo può rendere instabile la soluzione generando delle criticità.

Lo svariato numero di collocazioni possibile per la soluzione dei polinomi di 9° grado a tre incognite, ha reso necessari diversi tentativi infruttuosi, fino a quando lo scrivente ha scoperto l’unica che rendeva tutti i coefficienti di peso positivi. La soluzione è stata ottenuta disponendo 62 punti al posto dei 125 richiesti dall’integrazione di Gauss-Legendre, e in questo modo si è reso possibile integrare col nuovo metodo anche le funzioni più complesse che si incontrano nel metodo degli elementi finiti.

Libro

La ricerca, aggravata dal malcondizionamento dei sistemi e ancora di più dalla tassativa condizione di accettare unicamente le soluzioni con coefficienti di peso interamente positivi, ha richiesto diversi anni di studio, ma alla fine i risultati sono stati pubblicati nel testo (di cui lo scrivente è autore): “Quadratura diretta degli integrali multipli”, edito dalla Libreria Internazionale Cortina di Padova con ISBN 88-7784-263-6.

L’opera presenta in modo dettagliato le soluzioni del problema a partire dai casi più semplici che sono stati oggetto della tesi di laurea, fino ai più sofisticati che sono stati affrontati nelle ricerche degli ultimi anni. Il libro è corredato da esempi di calcolo basati su una funzione test uguale per tutte le soluzioni, allo scopo di escludere che la precisione della quadratura fosse dovuta ad un caso fortuito.

Il testo si conclude con la verifica manuale delle soluzioni delle equazioni, risparmiando così al lettore il tempo necessario per il controllo delle medesime.

Il Prof. Vitaliani ha seguito gli esiti delle ultime ricerche mediante simulatori che hanno confermato sperimentalmente la validità della nuova integrazione che, sfruttando in pieno le potenzialità della quadratura di Gauss, la estende a funzioni di più variabili senza ricorrere a palliativi che ne rendono più laboriosa l’applicazione.