Metodi numerici per l'ingegneria civile M

2025 Q

UADERNO APPUNTI

S V

OFIA IGNALI Metodi numerici per l'ingegneria civile

1

2

C : M N 'I C -

ORSO ETODI UMERICI PER L NGEGNERIA IVILE

M 2

ODULO

By Michele Ruggeri

Introduzione

Si consideri un problema lineare: () =

Un'approssimazione della soluzione può essere:

() ()

() ≈ = +

Si ottiene un residuo: (; ) )

ℝ ; ; … ; = ( −

In questo modo il problema si sposta da trovare u a trovare i coefficienti a

moltiplicati per funzioni note e questo ci permette di semplificare le

derivate, siccome è un operatore di derivate queste derivate si scaricheranno

sulle funzioni note in modo da ottenere semplici equazioni matematiche,

rappresenta la funzione (nota) che rappresenta le condizioni iniziali del

problema.

Queste equazioni solitamente sono divise in una parte di soluzione omogenea

(approssimazione) e una parte che rappresenta la soluzione particolare del

problema con tutte le condizioni al contorno che possono essere omogenee o

meno. una volta inserite tutte queste componenti nel problema otteniamo un

residuo e attraverso diverse tecniche possiamo minimizzare questo residuo in

modo da ottenere delle soluzioni approssimate che si avvicinino il più

possibile alla soluzione del problema.

E .

SEMPIO

( ) ( )

− + 3 = 0 = − + 3

→

=0

(−1) = (1) = 1

La soluzione esatta () = exp (−1) = (+1) = 0

3

Con i monomi si avrebbe 1, , ,

(1 ) ( ) (1 ) ) (1 )

= ⏟

1 + − + + = − + (1 − + − + ⏟

1

( ) ( )

ℝ = − − 3 = 0

( )

ℝ , , , ≠ 0 = 1,2,3

1

=− 2

=0

→ 1

= 2

+ 3

=

= → = →

+ 3 =

+ 3 4

A questo punto una volta ottenuta l'approssimazione si può scrivere il residuo.

Il residuo deve essere posto pari a 0 per trovare la soluzione. A questo punto

devo utilizzare il metodo di collocazione, cioè, ottenere 3 equazioni per

risolvere il problema se al posto della inserisco cioè, colloco la

funzione invece che in x generico in 3 punti che non possono essere quelli al

bordo perché sennò sarebbe mal condizionata.

A questo punto devo estrarre dall'equazione del residuo i coefficienti che

moltiplicano

E così divido le equazioni in questi tre termini per collocarli nei nodi scelti

devo utilizzare una matrice. Per risolvere l'operazione su MATLAB si può usare

l'operatore backslash

E . Bisogna creare un vettore colonna in perchè sennò non funziona

RRORE

La soluzione vera però è :

5

Per paragonarla con la soluzione esatta:

Le soluzioni sono diverse però si avvicinano molto.

Condizioni al contorno

Se il problema richiede delle condizioni al contorno periodiche nello spazio

seni e coseni della serie di Fourier le soddisfano automaticamente. Quindi le

condizioni al contorno in questo caso sono sempre soddisfatte e per questo non

vanno indicate.

Per problemi non periodici una scelta naturale ricade sui polinomi di Chebichev

che però non soddisfano le condizioni al contorno. Si utilizzano questi, invece

che i monomi, perché hanno diverse proprietà molto interessanti e hanno il

grosso vantaggio di portare a soluzioni molto accurate e per questo spesso

vengono indicati come gli unici da usare.

I primi 5 polinomi di Chebichev sono riportati sotto, definiti in [−1; 1]:

 ()

=

 ()

= 2 − 1

 ()

= 4 − 3

 ()

= 8 − 8 + 1

Queste sono le soluzioni calcolate perché in realtà i polinomi sono definiti da

una formula ricorsiva. È complicato utilizzare questi polinomi nei nostri

problemi, In primis perché i problemi tipici non sono inscritti nell'intervallo

[−1,1]. 6

Un trucco per soddisfare le condizioni al contorno è il seguente: se si

considera con si può imporre

= 0 (1) = , (1)

=

Se si considera la stessa diventa Dunque, la soluzione

∑ (1)

≠ 0 = 0.

deve essere vista come: () = () + ()

 ()

() =

 () =

Le condizioni al contorno sono dunque soddisfatte con:

(1) = (1) = 0

Cioè, la soluzione u è data da una parte di soluzione omogenea e una parte di

soluzione particolare, in questo caso si utilizza la notazione di Einstein la

stessa cosa significa: () = =

Quindi il trucco è scegliere un set di polinomi che soddisfi le soluzioni

omogenee quindi pari a 0 e delle soluzioni che invece soddisfano le condizioni

al contorno. (1 ) (2

= 1 + − ⋅ ⏟

1 + ⋅ ⏟

+ − 1)

È possibile dunque combinando i polinomi di Chebichev ottenere delle funzioni

di base che siano in grado di soddisfare le condizioni omogenee al contorno

(−1) = (1) = 0. () ()

≡ − 1 = 1,2, …

() ()

= −

Questo approccio detto anche ricombinazione della base porta ad una matrice nei

coefficienti con un miglior condizionamento rispetto agli stessi polinomi di

Chebichev poiché gli stessi hanno un gradiente elevato vicino ai bordi e

tendono ad oscillare per derivate di ordine alto. I polinomi di Chebichev sono

molto stabili. 7

N . Il problema del marchio condizionamento inizia a notarsi per oppure

OTA ~100

per derivate di terzo quarto ordine:

= −1 = −1=0

= − = − =0

Bisogna stare attenti a considerare i primi polinomi se risultano poi essere

nulli.

Classi di equazioni differenziali

E (H , P , L ).

QUAZIONI ELLITTICHE ELMHOLTZ OISSON APLACE

+ + (, ) = (, )

∂ u

+ = (, )

∂x

+ =0

E ( ).

QUAZIONI PARABOLICHE CONDUZIONE DEL CALORE

= (, ) +

E ( , ).

QUAZIONI IPERBOLICHE ONDE VIBRAZIONI

= (, ) +

Queste 3 classi di equazioni, scritte in forma generica, rappresentano

tantissime equazioni differenziali presenti nel mondo dell'ingegneria.

Helmholtz ha termine derivato e sorgente, Poisson ha termine sorgente, Laplace

non ha nessuno dei due, annullando la derivato seconda in tutto si ottiene il

termine monodimensionale, utilizzando invece anche la derivata terza di ottiene

il caso tridimensionale.

Le derivate paraboliche sono date dalla derivata prima del tempo mentre quelle

iperboliche sono caratterizzate dalla derivata seconda e sono molto presenti

nel caso ad esempio della sismica delle costruzioni.

Questi sono problemi ben posti ma mancano le condizioni al contorno e queste

equazioni con le condizioni al contorno (ellittiche solo condizioni al contorno

spaziali, mentre le altre sono anche temporali quindi hanno due set di

condizioni al contorno). Queste equazioni rappresentano le forme forti perché

hanno appunto le condizioni al contorno e l'equazione differenziale ed esistono

anche le forme deboli, che viene utilizzata ad esempio negli elementi finiti

(perde in accuratezza ma per praticità per problemi complessi è molto adatta).

8

Metodo dei residui pesati

Nel primo esempio attraverso il metodo di collocazione abbiamo trovato un

risultato a cui appartiene un residuo e da qui partono i metodi dei residui

pesati consistono in una serie di metodi classici atti a minimizzare il residuo

di un problema differenziale o integrale.

Se i polinomi di base A priori le condizioni del contorno questi metodi

()

calcolano i coefficienti incogniti imponendo condizioni nella forma

integrale. ℝ , ,…, = 0

Considerando nel dominio delle funzioni test/peso

[−1,1] ().

∈

M . Il metodo di collocazione in questo caso risulta una

ETODO DI COLLOCAZIONE

tipologia di metodo dei residui pesati attraverso il di Diack (?).

Se la funzione peso è il residuo diventa:

),

= ( −

( )

ℝ , , , … , = 0 = 1, … ,

Il numero di parametri deve essere uguale al numero di punti di collocazione

. Questo approccio è stato utilizzato nel primo esempio.

M G . Se la funzione peso è approssimazione

()con ()

ETODO DI ALERKIN =

il problema diventa =

= =

N . È più difficile risolvere un problema alla Galerkin rispetto al metodo di

OTA

collocazione/pseudo spettrale, dato che ho una moltiplicazione della funzione

poiché è una funzione è più semplice valutarla che integrarla. Tuttavia, per un

algoritmo nel tempo (risoluzione interattiva) il metodo Galerkin risulta più

veloce poiché le matrici sono a banda.

E . equazione differenziale del secondo ordine con

SEMPIO − = − cos − cos 2

condizioni al contorno periodiche.

1 3 9

() = = − → = − cos − cos 2

2 2 2

Se si ha così si sceglie una base trigonometrica con

= 2 → = cos + cos 2

le giuste condizioni al contorno. Quello che si deve fare ora per applicare

9

Galerkin è valutare innanzitutto per un generico polinomio e poi valutarlo

per la nostra funzione test e integrarlo. Le derivate sono:

= cos = cos 2

= = = 1,2

Che diventa: 1 1

) )

( = − ( = −

2 2

) ) ) )

(1,1) = ( (1,2) = ( (2,1) = ( (2,2) = (

3

− 0 3 1 0

2

= =−

9 0 3

2

0 −

2

Si ha così una matrice diagonale poiché le funzioni approssimanti

trigonometriche sono funzioni ortogonali, poiché all’interno dell’integrale si

hanno due funzioni trigonometriche uguali di ordine diversi e quindi

l’integrale del loro prodotto è sempre nullo.

Il vettore risulta:

(1,1) ( )

= ⋅

3 9 1

= − cos − cos 2 → →= = → = cos + cos 2

1

2 2

(2,1) ( )

= ⋅

Principi variazionali

La forma debole rientra nel concetto di principi variazionali ed è il cuore

della formulazione variazionale del problema ed è quella di moltiplicare il

residuo pesato per una funzione test ed integrare per parti.

N . Si parla di formulazione variazionale solo per problemi di meccanica

OTA

delle strutture, dove è possibile definire un funzionale iper-elastico, mentre

per tutti gli altri ambiti della meccanica si parla di forma debole.

− () ⋅ + () = () (−1) = (1) = 0

La forma debole è:

() ⋅ + () − () = 0 (−1) = (1) = 0

10

D .

IMOSTRAZIONE

− () ⋅ + () = () (−1) = (1) = 0

Si riportano tutti i termini a sinistra e si impone la funzione nulla poi si

integra su tutto il dominio, trovando così una funzione dei residui pesati:

− () ⋅ +⋅− =0

Qui l’obiettivo è di indebolire la funzioni e lo si fa integrando per parti:

⋅⋅ + ⋅ ⋅ − ⋅ − ⋅ ⋅ =0

 variabile primaria

→

 variabile secondaria

→

 funzione nota

→

Così si passa dalla forma forte alla forma debole in pochi passaggi.

Considerando il metodo di Galerkin (non sempre coincide con il metodo di Ritz

come in questo caso): () = =

() ⋅ + () − () = 0

In forma compatta: =

Ma può essere risolto per inversione matriciale una volta valutati gli

integrali. 11

Interpolazione

Un polinomio approssimato può essere definito come:

() ) ()

= (

Questa è anche detta funzione interpolante di Lagrange dove prendono il

()

nome di funzioni cardinali che soddisfano la relazione e hanno valore

=

unitario nel singolo punto e sono nulle negli altri punti. Scritte anche come:

−

()

= Π , −

P R . Questa interpolanti funzionano bene con gli zeri dei polinomi

ROBLEMA DI UNGE

di Chebichev, cioè quando il polinomio interseca l’asse e questi zeri sono

ricorsivi, come griglia (definizione più comoda per utilizzare MATLAB):

(2 − 1)

= − cos = 1,2, … , + 1

2( + 1)

V MATLAB: _01_R . ; _02_R .

EDERE I CODICI EX UNGE M EX UNGE M

Polinomi di Chebichev

I polinomi di Chebichev, come altri polinomi ortogonali, sono definiti

nell’intervallo per questo i problemi andranno adimensionalizzati in

[−1; 1],

questo dominio. I punti della griglia sono distribuiti secondo:

[−1;

= − cos = 0, … , ∈ 1]

Questa griglia considera anche i punti di bordo (detta di Gauss-Lobato). Le

equazioni cardinali sono: ()

1− 2 1

() (−1) ()

= ⋅ ⋅ = ⋅ ⋅

−

 o

= 2 = 1

±

 o

= =2 =1

I polinomi di Chebichev hanno una definizione ricorsiva:

() () ()

=1 = = 2 −

12

I coefficienti che permettono il calcolo della derivata prima dei polinomi

risultano: 1 + 2

⎧ − ==0

6

⎪ 1 + 2

⎪ ==

⎪ 6

( )

|

= =

− = ∧ 0<<

⎨ 2 1−

⎪

⎪ (−1)

⎪ ≠

−

⎩

Le derivate successive vengono calcolate con la moltiplicazione matriciale:

( ) ( ) ( )

=

Ed è anche questo uno dei vantaggi nell’utilizzo dei polinomi di Chebichev.

Perché, se il polinomio interpolante di vale:

()

) () ) () ()

= ( = ( → ≈ ()

Allora la derivata della funzione può essere scritta come: ()

= ⋅

( )

è

Ma se si vuole calcolare la derivata della funzione in un determinato punto

e non in maniera analitica si ha: ()

| |

= ⋅ → = ⋅

Ciò significa che, se io ho i valori della funzione in certi punti e la

moltiplico in questa matrice trovo la funzione in tutti gli altri punti.

S MATLAB: _03_ _ .

I VEDA IL CODICE EX BASIS CHEB M

Polinomi di Legendre

I punti della griglia (con punti esterni = Gauss-Lobato) sono, in questo caso

con definizione non ricorsiva: [−1,1]

= .1 = 1 = =0 ∈

Le funzioni cardinali sono: ()

1−

()

= − ⋅

( + 1) −

13

I coefficienti della derivata sono:

( + 1)

⎧− ==0

4

⎪ ( + 1)

⎪

==

( ) | 4

= =

⎨ 0 = ∧0 < <

⎪ ( )

⎪ ≠

−

⎩

S MATLAB: _04_ _ .

I VEDA IL CODICE EX BASIS LEG M

Integrazione di Gauss

I metodi di collocazione sono detti pseudo-spettrali perché, se si utilizzano

griglie di punti interpolanti ottimali allora il metodo di collocazione

coincide con il metodo di Galerkin i cui integrali sono valutati numericamente

alla Gauss: ) ()

() ≈ ( ← =

Quindi il tutto consiste in una sommatoria della funzione nei punti

moltiplicata per i pesi, facili da calcolare. Solitamente gli integrali sono

valutati numericamente con i punti di griglia che corrispondono agli zeri dei

polinomi utilizzati per l'approssimazione funzionale (griglia di Chebichev –

polinomi di Chebichev, griglia di Legendre – polinomi di Legendre).

C'è una relazione tra il metodo di collocazione e il metodo di Galerkin.

Considerando le ascisse della quadratura gaussiana come punti di collocazione

si ottiene il metodo pseudo-spettrale. Considerando il sistema discreto secondo

il metodo di collocazione dato un problema lineare:

) ( )

() = ( =

=

 ( )

=

 )

= (

Se si pre-moltiplica tutto per ( ):

( ) ( )

=

Considerando: ( ) () () ( ) ()()

= ≈ = ≈

14

Si ottiene il sistema risolvente secondo il metodo di Galerkin:

() () ()ℝ

− () =

Si è ottenuto il sistema discreto secondo Galerkin che risulta analogo

=

al metodo pseudo-spettrale . Lo stesso si può ottenere se si introducono

=

in questa dimostrazione anche le condizioni al contorno del problema.

P L . Per i polinomi di Legendre (con punti estremi) se si ha:

OLINOMI DI EGENDRE

= −1 = 1 =0

Con pesi: 2

()

= = ( + 1)

Questi sono lo standard per gli elementi finiti spettrali poiché questi

polinomi sono ortogonali rispetto al prodotto interno con valore peso unitario.

Quindi comode e stabili per il calcolo.

I polinomi di Legendre (con punti interni) se hanno:

= = 0

Allora hanno pesi pari a: 2

= 1−

Si utilizza per gli elementi finiti standard con codici commerciali e no, ed è

chiamata griglia di Gauss-Legendre, non Lobato poiché i punti estremi non ci

sono.

P C . I polinomi di Chebichev (con punti interni) hanno pesi:

OLINOMI DI HEBICHEV 2 1 − cos

= sin ⋅

+1

Con punti interni della griglia:

= = 1,2, … ,

+1 15

Questi pesi sono validi solo se si considera il cambio di coordinate e

= cos()

l’integrale infatti diventa:

() = (cos()) sin() = cos sin

S MATLAB: _05_ _ _ .

I VEDA IL CODICE EX INT CHEB LEG M

Soluzione di un problema

Si consideri il problema al contorno: (−1)

+ () = () → = (1) =

Si consideri la griglia (di Chebichev):

= cos = 1,2, … , − 2

−1

Con approssimazione: ()

≈ = = 1,2, … ,

Le equazioni discrete sono:

( ) ) ( )

= + (

, = 1, … , − 1 = 1,2, … ,

)

= (

Le equazioni sono state discretizzate perché:

() = = + () = ()

) ( ) ) () )

( = + () → = + ( = ( → =

Questi polinomi così non stanno a significare nulla e vano aggiunte le

condizioni al contorno. Si sono saltate le righe iniziali per le condizioni al

contorno, perché sono miste e non omogenee, cioè non sono nulle e una riguarda

la funzione e l’altra la derivata (non è sbagliato è un approccio diverso):

(−1)

= =

, (1)

= =

,

Quindi la soluzione si trova come: =

16

Soluzione analoga al problema può essere ottenuta considerando una

approssimazione utilizzando le funzioni cardinali:

()

= =

S MATLAB: _06_ _ . ( B (?)

I VEDA IL CODICE EX BVP CHEB M FUNZIONE PREPARATA DA TESTO DI OYD PER

C )

IMPLEMENTARE LE FUNZIONI CON I POLINOMI DI HEBICHEV IN AUTOMATICO

Ricombinazione della base

Se il problema differenziale ha condizioni al contorno non omogenee si può

effettuare la ricombinazione della base scegliendo una funzione che soddisfa le

condizioni non omogenee al contorno ed una nuova variabile che soddisfa solo le

condizioni al contorno omogenee.

Dato con condizioni e considerando

(−1) (1)

() = = , = −(1) = () = () +

.

(), ()

() = () −

Il problema da risolvere è e consider