elementi finiti

7/28/2019 elementi finiti

1/58

Metodo degli Elementi Finiti

Matteo Boccafoli

March 27, 2007

1

3

2

3

1

3

2

1

1

2 2

3

1

2

3

u(e)1

u(e)3

u(e)2

u

(e)

e e e

e


2/58

ii

Copyright (c) 2005-2007 Matteo BoccafoliPermission is hereby granted, free of charge, to any person obtaining a copy of this software and

associated documentation files (the Software), to deal in the Software without restriction, includingwithout limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell

copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to thefollowing conditions:

The above copyright notice and this permission notice shall be included in all copies or substantialportions of the Software.

THE SOFTWARE IS PROVIDED AS IS, WITHOUT WARRANTY OF ANY KIND, EXPRESSOR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABIL-ITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALLTHE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OROTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARIS-ING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHERDEALINGS IN THE SOFTWARE.


3/58

Contents

Indice generale iv

Indice figure v

1 Introduzione 1

2 Fondamenti teorici di FEM 1D 32.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2 Richiami e complementi di analisi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.2.1 Integrale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2.2 Integrazione per parti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2.3 Differenziale totale e parziale di una funzione a due variabili . . . . . . . . . . . . . 42.2.4 Cm() . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2.5 Spazio vettoriale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2.6 Metrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2.7 Norma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2.8 Integrale di Lebesgue e spazi Lp . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.2.9 Prodotto interno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2.10 Spazi di Hilbert . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.2.11 Derivata debole . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.2.12 Spazi di Sobolev . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.3 Formulazione variazionale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3.1 Funzionale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3.2 Il simbolo variazionale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.3.3 Formulazione variazionale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.3.4 Minimizzazione di un funzionale . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.3.5 Relazione tra le formulazioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.4 Condizioni al Bordo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.4.1 Alcune formulazioni variazionali . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.4.2 Riassumendo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.5 Metodo variazionale di approssimazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.5.1 Funzioni lineari a tratti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.5.2 Metodo di Ritz-Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.6 Stima dellerrore . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.6.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.6.2 Analisi dellerrore . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.7 Ricapitolando . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3 FEM 1D in pratica 253.1 Formulazione variazionale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.2 Funzione di approssimazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.3 Generiamo la matrice di stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.4 Generiamo il vettore di carico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.5 Risolvere il sistema lineare . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.6 Un esecuzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

iii


4/58

iv CONTENTS

3.7 Stima dellerrore . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.8 Consigli per lo sviluppo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

4 Fondamenti teorici di FEM 2D 31

4.1 Richiami e complementi di analisi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 314.1.1 Insieme aperto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 314.1.2 Derivazione e continuita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 324.1.3 Curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.1.4 Integrazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.1.5 Formula di Green . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.1.6 Teorema della divergenza e del gradiente . . . . . . . . . . . . . . . . . . . . . . . . 40

4.2 Formulazione Variazionale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.3 Alcune formulazioni variazionali . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

4.3.1 Equazione di Poisson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 414.3.2 Variante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 414.3.3 Problema di Neumann . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

4.3.4 Tabella riassuntiva . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 424.4 Metodo variazionale di approssimazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 424.4.1 Triangolazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 434.4.2 Metodo di Ritz-Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

A Riconoscimenti 47


5/58

List of Figures

2.1 Larea rosa rappresenta la distanza d1(f, g) in [a, b] . . . . . . . . . . . . . . . . . . . . . . 62.2 funzione u(x) = |x| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.3 La funzione rosa vh e lineare a tratti. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.4 (A) funzione di base 2(x), (B) una lineare a tratti con relative funzioni di base . . . . . . 20

2.5 Cerchio approssimanto da 6 triangoli (A) e da 9 (B) . . . . . . . . . . . . . . . . . . . . . 222.6 Errore in 1D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.1 Limmagine (A) mostra la decomposizione di in in (B) vengono poste in evidenza lefunzioni che generano un segmento ui . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.2 Funzione un lineare a tratti rossa e u soluzione reale in verde, con 11 nodi . . . . . . . . . 293.3 Funzione un lineare a tratti rossa e u soluzione reale in verde, con 21 nodi . . . . . . . . . 30

4.1 (A),(B) e (C) rappresentano nelle varie metriche B1(0, r), B2(0, r) e B(0, r) . . . . . . . 314.2 Linsieme B e un insieme connesso, mentre A non lo e. . . . . . . . . . . . . . . . . . . . . 324.3 Integrazione di f(x) secondo Riemann, in azzurro i supxIj in rosa gli infxIj . . . . . . . 374.4 Integrazione di f(x, y) secondo Riemann, in azzurro i sup(x,y)IjJi in rosa gli inf(x,y)IjJi 374.5 In rosa abbiamo caratterizzato la funzione f in [a, b]

[c, d], mentre la superficie A e

ristretta alla solo proiezione sul piano di f . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.6 Rappresentazione del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394.7 Alcune triangolazioni (tratta da [BBM05]) . . . . . . . . . . . . . . . . . . . . . . . . . . . 434.8 In rosa le singole (e) relative ad un elemento e, tutte le (e) generano (con opportune

costanti che troveremo) la futura funzione u relativa ad e (u(e)) denotata dal triangolo blu;in azzurro vi sono le (e) che generano u(e) . . . . . . . . . . . . . . . . . . . . . . . . . 43

v


6/58

vi LIST OF FIGURES


7/58

Chapter 1

Introduzione

La neccessita, o leconomicita, di poter simulare degli avvenimenti fisici, o valutare modelli differentimediante un elaboratore elettronico, solleva un interesse nella richiesta di software in settori come aereo-nautica, economia, edilizia, fisica, metereologia, sociologia, etc..

Osservare lelasticita di un corpo, analizzare il comportamento di una struttura edilizia durante unavvenimento sismico, gestire il flusso in una rete, sono tutti problemi problemi fisici che si possono risolvereattraverso una particolare famiglia di equazioni, le equazioni differenziali.

F(t,u,u, u) = 0 u(t) = et+sin(t/ log(t))/2

Focalizziamo la nostra attenzione nel trovare la funzione u(t) la cui derivata i-esima deve essere ugualead una funzione nota.

Dato che queste equazioni quasi sempre - o almeno, nei casi interessanti - non sono risolubili inmaniera algebrica1, ci concentreremo prevalentemente nella loro soluzione numerica, attraverso uno deivari metodi a nostra disposizione, chiamato metodo degli elementi finiti (FEM).

In sintesi il metodo e composto dei seguenti passi:

Riformuliamo il problema in forma variazionale, in modo da focalizzare la nostra attenzione sullavariazione dellincognita del problema.

Fatto cio, andiamo a ricavare delle soluzioni dellequazione differenziale solo in tutti i punti ap-partenenti a quella che chiamiamo discretizzazione del dominio, ovvero un campione del dominio.

Infine, ricostruiamo una funzione approssimata dalle informazioni elaborate, la quale sara piu vicinaalla soluzione reale allaumentare dei punti scelti per rappresentare il dominio.

Osservazione 1.1. E possibile trovare degli esempi di software e applicazioni pratiche che usano FEMin un motore di ricerca di immagini le seguenti stringhe: finite element method, fem simulation, femmagnetic, fem naval, fem medical e fem fluid.

1O hanno forme indeterminate come u = e2t2

oppure e troppo costoso cimentarsi in una formula lunga due pagine

1
http://images.google.com/images?q=finite+element+methodhttp://images.google.com/images?q=fem+simulationhttp://images.google.com/images?q=fem+simulationhttp://images.google.com/images?q=fem+magnetichttp://images.google.com/images?q=fem+magnetichttp://images.google.com/images?q=fem+magnetichttp://images.google.com/images?q=fem+magnetichttp://images.google.com/images?q=fem+navalhttp://images.google.com/images?q=fem+navalhttp://images.google.com/images?q=fem+medicalhttp://images.google.com/images?q=fem+fluidhttp://images.google.com/images?q=fem+fluidhttp://images.google.com/images?q=fem+medicalhttp://images.google.com/images?q=fem+navalhttp://images.google.com/images?q=fem+magnetichttp://images.google.com/images?q=fem+magnetichttp://images.google.com/images?q=fem+simulationhttp://images.google.com/images?q=finite+element+method


8/58

2 CHAPTER 1. INTRODUZIONE


9/58

Chapter 2

Fondamenti teorici di FEM 1D

Nel seguente capitolo tratteremo prevalentemente dei concetti in modo il piu possibile astratto da tutti i

possibili casi di equazioni differenziali. Per dedicare il capitolo successivo ad un esempio pratico.

2.1 Introduzione

Alcune equazioni differenziali con condizioni al bordo possiamo risolverle analiticamente trovando unasoluzione forte, o detta anche soluzione esatta. Quando percorriamo questa strada parliamo di formu-lazione forte del problema.

Definizione 2.1 (formulazione forte). Dato un intervallo I, trovare u C(2)(I) tale che

u = fu(0) = u(1) = 0

e detta formulazione forte.

Osservazione 2.1. La precedente definizione e legata ad un tipo di condizioni al bordo, detta di Dirichlet,ma possiamo anche avere altri tipi come verra mostrato in seguito (eq. 2.24).

Se il procedimento matematico trova ostacoli in funzioni complesse come f(x) = ex2

, il cui integralee indeterminato, possiamo avvicinarci ad una possibile soluzione, detta debole, approssimandola; perprocedere ci sono due punti fondamentali:

Un passaggio alla forma variazionale della equazione differenziale data. Determinare un soluzione approssimata usando un metodo variazionale come il metodo di Ritz, il

metodo di Galerkin, o altri metodi.

Prima di procedere richiameremo delle nozioni presentate in altri corsi, in seguito parleremo di formu-lazione variazionale al fine di trattare il metodo variazionale dapprossimazione di Galerkin. Destiniamo

alle ultime sezioni di questo capitolo temi come la stima dellerrore e gestione delle condizioni al bordo.

2.2 Richiami e complementi di analisi

2.2.1 Integrale f(x) + g(x) dx =

f(x) dx +

g(x) dx (2.1)

ca

f(x) g(x) dx =ba

f(x) g(x) dx +cb

f(x) g(x) dx a b c (2.2)

a f(x) dx = a

f(x) dx (2.3)

3


10/58

4 CHAPTER 2. FONDAMENTI TEORICI DI FEM 1D

2.2.2 Integrazione per parti

Nella formulazione variazionale di equazioni differenziali viene fatto grande uso dellintegrazione per parti.

f(x) g(x) dx = f(x) g(x) dx f(x) g(x) dx (2.4)ba

f(x) g(x) dx = f(b)g(b) f(a)g(a) ba

f(x) g(x) dx (2.5)

2.2.3 Differenziale totale e parziale di una funzione a due variabili

Definizione 2.2 (derivata parziale). Si chiama derivata parziale rispetto ad x della funzione f(x, y)fissato y = y0

limx0

f(x + x,y0) f(x, y0)x

e la indichiamo conf

xo fx.

In modo analogo, la derivata parziale rispetto ad y della funzione f(x, y) fissato x = x0 e

limy0

f(x0, y + y) f(x0, y)y

e la indichiamo conf

yo fy.

Definizione 2.3 (differenziale parziale). Si chiama differenziale parziale di una funzione a due variabilif(x, y) rispetto alla variabile x relativamente al punto (x, y) ed allincremento x

f

xx

Analogamente rispetto a y

f

yy

Definizione 2.4 (differenziale totale). Si chiama differenziale totale di una funzione a due variabilif(x, y) relativamente al punto (x, y) ed agli incrementi x e y

df =f

x

x +f

y

y

Osservazione 2.2. Per incrementi x e y abbastanza piccoli f df.

2.2.4 Cm()

Definizione 2.5 (funzione continua). Una funzione f e detta continua in un punto c [a, b] e limxc f(x) f(c). Diremo che la funzione e continua in un intervallo I, se e continua per ogni punto di I.

Definizione 2.6 (spazio funzionale). Un insieme di funzioni con un operazione e detto spazio funzionaleo spazio di funzioni.

Esempio 2.1. Linsieme delle funzioni continue in un dominio con operazioni classiche tra funzioni,

denotato daC0() = {u|u continua in }

e uno spazio funzionale.


11/58

2.2. RICHIAMI E COMPLEMENTI DI ANALISI 5

Esempio 2.2. Anche lo spazio con anche la derivata prima continua in

C1() = {u|u, u continue in }

e anche lo spazio con tutte le derivate fino allordine m continue in

Cm() = {u|u, u, . . . , u(m) continue in }sono spazi funzionali.

Osservazione 2.3. Le funzioni elementari ex e sin(x) sono in C.

2.2.5 Spazio vettoriale

Definizione 2.7 (spazio vettoriale). Dato un insieme V e unoperazione + : V V V, (V, +) e ungruppo commutativo, i.e. u, v, w V vale:

1) (u + v) + w = u + (v + w)

2) e + v = v + e = v

3) v + v = 0V = 0

4) v + w = w + v

e data un operazione con un campo K, definita come : KV V, e siano h, k K, valgono le seguenteproprieta:

1) h (v + w) = h v + h w2) (h + k) v = h v + k v3) (h

k)

v = h

(k

v)

4) 1 v = vSe per una struttura (V, +, ) si verifica queste proprieta allora lo chiameremo spazio vettoriale o

lineare su K o K-spazio vettoriale, scrivendolo K s.v..Esempio 2.3. C1 e uno spazio vettoriale

Definizione 2.8. La dimensione di uno spazio vettoriale e il numero di vettori necessari e sufficientiper poterlo rappresentare. Denoteremo la dimensione di V K s.v. con dim V = n. Se n e finito alloraparliamo di spazio finito dimensionale, altrimenti diciamo che ha infinite dimensioni.

2.2.6 Metrica

Definizione 2.9 (distanza). Dato un insieme X chiamiamo distanzao metrica la funzione d : XX Rtale che

1. x, y X d(x, y) 02. d(x, y) = 0 x = y3. x, y X d(x, y) = d(y, x)4. x, y, z X d(x, y) d(x, z) + d(z, y)

Definizione 2.10 (spazio metrico). La struttura algebrica (X, d) dati dalla definizione precedente echiamata spazio metrico.

Esempio 2.4. In un intervallo I = [a, b] definiamo un distanza tra due funzioni come d1(f, g) = ba |f(x)g(x)| dx.

1. f(x), g(x) (I,R) d1(f, g) 0


12/58


f(x)

g(x)

a b

Figure 2.1: Larea rosa rappresenta la distanza d1(f, g) in [a, b]

2. d1(f, g) = 0 f = g

3. f, g (I,R) d1(f, g) = d1(g, f)

4. f, g, h (I,R) d1(f, g) d1(f, h) + d1(h, g)

Intuitivamente sembra tutto ovvio, e bene comunque svolgere almeno una dimostrazione. Prendiamo ilpunto 2 d1(f, g) = 0 f = g.

Dimostriamo un verso. Se d1(f, g) = 0 allora per definizione di d1 abbiamoba

|f(x) g(x)| dx = 0;osservando che

I(x) dx = 0 solo se (x) = 0, quindi |f(x) g(x)| = 0 che e vero se f = g. Laltro

verso e esattamente il percorso inverso, se f = g o anche f g = 0 a maggior ragione |f(x) g(x)| = 0,integrando entrambi i lati dellequazione otteniamo ba |f(x) g(x)| dx = 0 ovvero d1(f, g) = 0.Esempio 2.5. Altri esempi di distanze:

d2(f, g) =

ba

|f(x) g(x)|2 dx1/2

dp(f, g) =

ba

|f(x) g(x)|p dx1/p

d(f, g) = maxxI

|f(x) g(x)|

Osservazione 2.4. La distanza d2(f, g) e la piu usata, perche ha il prodotto scalare tra funzioni.

2.2.7 Norma

Definizione 2.11 (norma). Dato V un K s.v., chiamiamo norma su V un applicazione : V Rtale che:

1. x V x 0

2. x = 0 x = 0

3. x V, a K a x = |a| x

4. x, y V x + y x + y detta disuguaglianza triangolare

Osservazione 2.5. x y definisce una distanza tra due punti x, y


13/58


Esempio 2.6. Norma in C0():

u(x) = d(u, ) = maxx

|u(x)|

dove e usato come elemento nullo dello spazio funzionale C0

.

u(x) u(x) 0u(x) = 0 u(x) = 0

a u(x) = |a| u(x) a Ru(x) + v(x) u(x) + v(x)

Osservazione 2.6. La definizione di norma e strettamente associata ad uno spazio vettoriale a differenzadella distanza.

Definizione 2.12 (spazio normato). Si chiama spazio normato la struttura algebrica data dalla quaterna(V, +, , ), con (V, +, ) uno K spazio vettoriale.Esempio 2.7. Oltre al precedente esempio, anche Cm e uno spazio normato con rispettiva distanza

d(u, v) = max0|a|m

maxx

D|a|u| D|a|v|e norma

u(x) = max0|a|m

maxx

D|a|u|La norma e un concetto astratto che ci consente di misurare o porre a confronto elementi di uno stesso

spazio.

Esempio 2.8. Se abbiamo Rn possiamo definire la norma

Lp = xp = n

i=1|xi|p

1/p

Con p = 1 e detta norma di Manhattan, se p = 2 e detta norma euclidea ed infine abbiamo la norma delmassimo o di Chebichev con p = .Teorema 2.1 (Weierstrass). Uno spazio normato X ha dimensione vettoriale finita se e solo se ad ognisuccessione {xn} X e possibile ricavare una sottosuccessione convergente in X.Definizione 2.13 (successione di Cauchy). Dato uno spazio normato X, una successione {xn} X sichiama di Cauchy se xn xm 0 con n, m Definizione 2.14 (spazio di Banach). Uno spazio normato X e detto spazio di Banach, o completo, setutte le sue successioni di Cauchy a valori in X convergono in X.

2.2.8 Integrale di Lebesgue e spazi Lp

Lintegrazione classica di Riemann-Cauchy, viene presentato nel corsi di base di analisi matematica. Lecarenze di alcuni punti utili per un analisi funzionale ha creato la necessita di introdurre uno strumentopiu potente e flessibile.

Osservazione 2.7. Abbiamo che gli integrali definiti secondo Riemann sono compresi nella definizione diLebesgue ma non e vero il viceversa. O in altre parole una funzione integrabile secondo Riemann lo eanche secondo Lebesgue, il viceversa non vale.

Definizione 2.15 (Spazi Lp). Denotiamo con Lp linsieme delle funzioni I R C il cui modulo hapotenza p-esima integrabile secondo Lebesgue.

Lp(I) = {f : I R|I

|f(x)|p dx < +} (2.6)

Osservazione 2.8. Per p 1 gli spazi sono dotati di norma, quando p = 2 anche di un prodotto scalare.Si presti attenzione allanalogia con gli spazi euclidei [?] la norma generalizza il concetto di modulo, edil prodotto scalare tra funzioni generalizza il concetto di prodotto scalare tra vettori.


14/58


2.2.9 Prodotto interno

Definizione 2.16. Data una qualsiasi coppia di elementi x, y V Ks.v. il prodotto interno, o prodottoscalare, e un applicazione x, y : V V K tale che:

1. x, x 0 x V2. x, x = 0 x = 0

3. x + z, y = x, y + z, y x, y, z V4. x, y = x, y x, y V, K5. x, y = y, x x, y V

Esempio 2.9. Norma indotta dal prodotto scalare in Rn

x = x1, x2, . . . , xnt Rn y = y1, y2, . . . , ynt R

n

x, y = xt y =n

i=1

xi yi (2.7)

Esempio 2.10. Prodotto interno in L2

f, g =I

f(x) g(x) dx

Sono rispettate :

f, f

0

f, f

= 0

f = 0

f, g = g, f f, g = f, g = f, g

f, g + w = f, g + f, w

Osservazione 2.9. I

u(x)v(x) dx =I

f(x) v(x) dxu, v = f, v

Osservazione 2.10. denotiamo u =

u, u = u, u 12 e u2 = u, u.Teorema 2.2 (disugualianza di Cauchy-Schwarz).

|u, v| u, u 12 v, v 12 = u v (2.8)

Dimostrazione 2.2.1. Ci costruiamo la funzione quadratica:

f() = u + v, u + v per def. 2.16 3.= u, u + v + v, u + v per def. 2.16 3.= u, u + u, v + v, u + v per def. 2.16 3.= u, u + u, v + v, u + v, v per def. 2.16 4.= u, u + u, v + v, u + 2v, v per def. 2.16 5.= u, u + u, v + u, v + 2v, v per oss. 2.10= u2 + 2u, v + 2v2 a 2 + 2 b + c


15/58


Ora calcoliamo il discriminante

0b2

4ac

0

(2 u, v)2 4 v, v u, u 04 u, v2 4 v, v u, u 0

4 u, v2 4 v, v u, uu, v2 v, v u, uu, v2 v2 u2u, v2

v2 u2

|u, v | v u

Teorema 2.3. In ogni spazio vettoriale con prodotto scalare ha norma indotta dal prodotto scalare

u

=u, u

12

Dimostrazione 2.3.1. Per definizione di norma e di prodotto interno abbiamo che le prime tre proprietasono identiche, mentre per provare che con il prodotto scalare possiamo indurre una norma dobbiamoprovare

u + v = u + va maggior ragione se

u + v u + v = u + v2= (u + v, u + v 12 )2= u + v, u + v=

u + v, u

+

u + v, v

= u, u + v, u + u, v + v, v= u2 + 2u, v + v2 per Couchy-Swartz u2 + 2uv + v2= (u + v)2= (u + v) (u + v)

2.2.10 Spazi di Hilbert

Definizione 2.17. Uno spazio di Banach con prodotto interno e detto spazio di Hilbert

Osservazione 2.11. Rivedendo lesempio 2.10, possiamo affermare che L2 e uno spazio di Hilbert.

2.2.11 Derivata debole

Definizione 2.18 (derivata debole). Sia u : I = [a, b] R e C1(I) dico che u e derivabile in sensodebole su I se

g L2(I) |I

u = I

g

C1(I) (a) = (b) = 0 (2.9)in questo caso g si chiama derivata in senso debole, o derivata debole, di u.

Osservazione 2.12. Nella precedente definizione abbiamo che g = u

I

u dx = u(b)(b) u(a)(a) I

u dx = I

u dx

Osservazione 2.13. Sebbene non e stato specificato u C1([a, b]).


16/58


Proposizione 2.4. Se u e derivabile in senso classico allora e derivabile in senso debole.

Dimostrazione 2.4.1. Supponiamo u C1(I) u in senso classico. Moltiplichiamo per una funzione(x) C1(I) e con (0) = (1) = 0, e integrando:b

au dx = u(b)(b) u(a)(a) b

au dx per parti eq. 2.5 e eq. 2.9

ba

u dx = ba

u dx per definizione di derivata debole

g = u quindi u e derivabile anche in senso deboleOsservazione 2.14. Quindi abbiamo dedotto che se u e derivabile in senso classico, allora lo e anche insenso debole, il viceversa e falso.

Esempio 2.11. u(x) = |x| nellintervallo I =] 1, 1[, non ha derivata ovunque dato che u(0), vedi fig.2.2

f(x) = |x|

Figure 2.2: funzione u(x) = |x|

Seguendo la definizione 2.18, prendiamo C1([1, 1]]) con le condizioni al bordo (1) = (1) = 0.11

u(x)(x) dx =11

|x|(x) dx

=

01

(x)(x) dx +10

(x)(x) dx per eq. 2.2 (2.10)

= [x(x)]01 01

(x)(x) dx + [x(x)]10 10

(x)(x) dx per parti eq. 2.5

= 0

1(x)(x) dx 1

0 (x)(x) dx (1) = (1) = 0=

01

1 (x) dx 10

1 (x) dx svolto la derivata

= 01

(x) dx +

10

(x) dx per eq. 2.3

= 0

1(x) dx

10

(x) dx

raccolto (2.11)

se prendiamo

g(x) =

1 x [0, 1]

1 x

[

1, 0[

oppure

g(x) =

1 x ]0, 1]1 x [1, 0]


17/58

2.3. FORMULAZIONE VARIAZIONALE 11

allora lequazione 2.11 e equivalente a 11

g(x)(x) dx

Osservazione 2.15. Il passaggio dellequazione 2.10 e ottenuto per definizione di valore assoluto e pre-stando attenzione al dominio dove integriamo. Se integriamo in un intervallo positivo, x [0, 1], sappiamoche la x e positiva, viceversa in [1, 0] la x e negativa.

2.2.12 Spazi di Sobolev

Introduciamo i seguenti spaziH0(I) = {u L2(I)} = L2(I)

Con la derivata prima debole

H1(I) = {u L2(I)|u debole}=

{u

L2(I)

|g

L2(I) per cui Iu = Ig C1(I) (a) = (b) = 0}

Con valore nullo al bordoH10 (I) = {u H1(I)|u(a) = u(b) = 0}

Osservazione 2.16. se prendiamo un insieme di funzioni la cui derivata prima e continua a tratti, questoe un sottoinsieme di

V = {u C(I)|u continua a tratti} H1(I)

2.3 Formulazione variazionale

2.3.1 Funzionale

Useremo il termine funzionale per descrivere funzioni i cui argomenti sono funzioni, ma non e da con-fondere con la funzione composta dove loutput di una funzione diventa linput di una altra. Nel funzionalela funzione stessa e concepita come una variabile. Nel nostro caso useremo quasi sempre funzionali definitida integrali.

I(u) =

ba

F(x,u,u) dx

I(u) restituira un valore al variare della funzione u.

Esempio 2.12. Data la seguente equazione differenziale:

ddx

a(x) dudx

f(x) = 0 (2.12)

con condizioni al bordo di Dirichletu(0) = 0 u(1) = 0 (2.13)

possiamo vederla come

F : R3 R F(x, u(x), u(x)) = ddx

a(x) dudx

f(x) u : [0.0, 1.0] R

da cui il funzionale

I(u) =

10

F(x,u,u) dx

Definizione 2.19 (funzionale lineare). Un funzionale l(u) e detto essere lineare se e solo se soddisfa la

relazione:l( u + v) = l(u) + l(v)

per ogni scalare e in R e per ogni u e v in di uno spazio vettoriale.


18/58


Esempio 2.13. Per verificare che il seguente funzionale

l() =

1

0

f(x) (x) dx

e lineare, procediamo come segue:

l( + ) =

10

f(x) ((x) + (x)) dx

=

10

f(x)(x) + f(x)(x)) dx

=

10

f(x)(x) +

10

f(x)(x)) dx

=

10

f(x)(x) +

10

f(x)(x)) dx

= l() + l()

Definizione 2.20 (funzionale bilineare). Un funzionale B(u, w) e detto essere bilineare se e lineare inognuno degli argomenti u e w:

B( u + v, w) = B(u, w) + B(v, w)B(u, w + z) = B(u, w) + B(u, z)

Esempio 2.14. Un esempio di funzionale bilineare e dato da

B(, u) =

10

(x) u(x) dx

Verifichiamo la linearita per largomento

B( + , u) =10

((x) + (x)) u(x) dx

=

10

(x) u(x) + (x) dx

=

10

(x) u(x) + 10

(x) dx

= B(, u) + B(, u)

Verifichiamo la linearita per largomento u(x)

B(, u + w) =10

(x) (u(x) + w(x)) dx

=

10

(x)u(x) + (x)w (x)) dx

=

10

(x)u(x) + 10

(x)w(x)) dx

= B(, u) + B(, w)

Definizione 2.21 (funzionale simmetrico). Un funzionale B(u, v) e detto essere simmetrico se per ogniargomenti u e v:

B(u, v) = B(v, u)

Queste definizioni sono utili per poter riformulare il problema in altre forme, come puo mostrarelesempio seguente.


19/58


Esempio 2.15. Se moltiplichiamo lequazione 2.12 per una funzione di test v(x), la quale e due voltedifferenziabile e rispetta le condizione (al bordo) v(0) = 0 e v(1) = 0, otteniamo

0 = v(x) d

dx a(x) du

dx f(x) dx0 = v d

dx

a dudx

vf dx

integrando lequazione ottenuta e usando le condizioni al bordo, otteniamo la forma variazionale (eq.2.17) 1

0

v d

dx

a dudx

vf

dx = 0 per eq. 2.1 (2.14)1

0

v ddx

a dudx

dx

10

vf dx = 0 per parti eq. 2.5 (2.15)

10

adv

dx

du

dx dx + v(1)a(1)u(1) v(0)a(0)u(0) 1

0vf dx = 0 v(0) = v(1) = 0 (2.16)1

0

adv

dx

du

dxdx

10

vf dx = 0 (2.17)

Possiamo ora ricavare

B(v, u) =

10

adv

dx

du

dxdx

l(v) = 10

vf dx

Osservazione 2.17. La funzione u e un incognita, ovvero non e una funzione data a priori e dobbiamo

ricavarla, mentre conosciamo f(x) e a(x) e v la scegliamo opportunamente.Osservazione 2.18. Abbiamo riformulato la forma variazionale 2.17 in termini di un funzionale bilineare edi uno lineare; quanto svolto ci permette di usare delle proprieta relative alla forma bilineare simmetricaal fine di risolvere il problema.

2.3.2 Il simbolo variazionale

Mostriamo ora un punto di vista diverso, tratto da [Red84], in cui studiamo il comportamente del fun-zionale F che denotiamo con F = F(x,u,u). Dato un x arbitrariamente fissato, F dipendera da u e u.Denotiamo con u la variazione di u, il cambiamento ottenuto sara uguale a un v con costante

u = v

e chiamato simbolo variazionale. In altre parole, con u identifichiamo una variazione ammissibile nellafunzione u(x) fissato un valore della variabile indipendente x.

Quando variera u di un u, abbiamo anche una variazione di F (dipendendo da u). Possiamo alloradire che al variare di u e u avremo una variazione di F, in analogia al differenziale totale per una funzionea due variabili (def. 2.4),

F =F

uu +

F

uu

possiamo notare una somiglianza tra la precedente variazione e il differenziale totale dF

dF =F

xdx +

F

udu +

F

udu

Essendo la variabile indipendete x fissata durante la variazione da u a u + u, possiamo porre dx = 0rendendo evidente il legame tra F e dF. agisce come un operatore differenziale; le leggi di somma, prodotto, rapporto, potenza e altro sono

completamente simili alle corrispondenti leggi di differenziazione.


20/58


Esempio 2.16. Se prendiamo F1 = F1(u) e F2 = F2(u), allora

1. (F1 F2) = F1 F2

2. (F1 F2) = F2F1 + F1F23.

F1

F2

=

F2F1 F1F2F22

4. (F1)n = nFn11 F1

Inoltre, loperatore variazionale puo commutarsi con gli operatori di integrazione e differenziale.

d

dx(u) =

d

dx(v) =

dv

dx= v = u =

du

dx(2.18)

b

a

u(x) dx = b

a

u(x) dx (2.19)

2.3.3 Formulazione variazionale

Definizione 2.22 (formulazione variazionale). Trovare us tale che C(1)([0, 1]) (0) = (1) = 010

u(x) (x) dx =10

f(x) (x) dx

e detta formulazione o problema variazionale, o debole.O denotata in forma di funzionale al variare di u come

F(u) = 10

u(x) (x) dx 10

f(x) (x) dx = B(u, v) l(u) (2.20)

Ottenuta con lo stesso procedimento dellequazione 2.17.

Osservazione 2.19. Il nome formulazione debole e calzante1 se si nota che nellequazione precedente lacontinuita richiesta per la u e di un ordine inferiore rispetto def. 2.1. Inoltre si osservi lanalogia con laderivata debole 2.18.

2.3.4 Minimizzazione di un funzionale

Definizione 2.23 (problema di minimizzazione). Lo spazio lineare

V :

{v

|v

C0([0, 1])

v

H10

}e la funzione lineare F : V R definita come

F(v) =1

2

10

v vdx 10

vf =1

2B(v, v) l(v) (2.21)

Trovare u V tale che F(u) F(v) per ogni v V

2.3.5 Relazione tra le formulazioni

Ora osserviamo le analogie tra le formulazioni del problema di partenza 2.1 nello stile seguito in [Joh90]p. 16-18.

Soluzione forte e variazionale Nellesempio 2.15 abbiamo mostrato che la soluzione u del problemain forma forte 2.1 e anche soluzione del problema in forma variazionale 2.22.

1come fa notare [Red84] a pag. 26


21/58


Soluzione variazionale e del problema di minimizzaione Ora mostriamo che i problemi in formavariazionale 2.22 e problemi di minimizzazione 2.23 hanno le stesse soluzioni.

Assumiamo che us sia una soluzione di 2.22, consegue che abbiamo H10 e possiamo porre w = ustale che = us + w e w

H10 .

F() = F(us + w)

=1

2

10

(us(x) + w(x)) (us(x) + w(x)) dx

10

(us(x) + w(x))f(x) dx

=1

2

10

us(x)us(x) + 2u

s(x)w

(x) + w(x)w(x) dx 10

us(x)f(x) + w(x)f(x) dx

=1

2

10

us(x)us(x) + u

s(x)w

(x) + w(x)us(x) + w(x)w(x) dx

10

us(x)f(x) dx 10

w(x)f(x) dx

=1

2

10

us(x)us(x) dx +

1

2

10

2us(x)w(x) dx +

1

2

10

w(x)w(x) dx 10

us(x)f(x) + w(x)f(x) dx

=1

210

us(x)us(x) dx 10

us(x)f(x) F(us)

+1

210

2us(x)w(x) dx +1

210

w(x)w(x) dx + w(x)f(x) dx

F(us)

Soluzione problema di minimizzaione e variazionale Ora vediamo che se um e soluzione di unproblema di minimizzazione allora se prendiamo una funzione H10 e R sicuramente

F(um) F(um + ) =1

2

10

(um + ) (um + ) dx

10

(um + )f dx

= 1210

(umum + 2um + 2) dx 10

umf + f dx

=1

2

10

umum dx +

10

um dx +1

2210

dx 10

umf dx + f dx

=1

2

10

umum dx +

10

um dx +1

2210

dx 10

umf dx 10

f dx

g()Studiamo il comportamento di g al variare di . Possiamo vedere che ce minimo con = 0 per cui

g(0) = 0 (2.22)

Osservando anche che

g(0) =10

um dx + 10

dx 10

f dx =

10

um dx

10

f dx (2.23)

valendo 2.22 e 2.23 allora abbiamo10

um

10

f dx = 010

um dx =

10

f dx

che corrisponde al nostro problema in forma variazionale 2.22.

Soluzione variazionale e unica Supponiamo che esistano due soluzioni variazionale e vediamo sesono la stessa.

Siano us1 , us2 H10 tali che10

us11 dx =

10

f 1 dx 1 H10 e10

us22 =

10

f 2 2 H10 .


22/58


Osservazione 2.20. Essendo per ogni 1, 2, us1 , us2 H10 puo anche essere che 1 = 2 = us1 us2 .Osservazione 2.21. Se a = b e c = d allora a c = b d

Nel primo passaggio uso losservazione 2.21 e nel secondo 2.20:

10

us11 dx

10

us22 dx =

10

f 1 dx 10

f 2 dx10

us1(us1 us2) dx

10

us2(us1 us2) dx =

10

f(us1 us2) dx 10

f(us1 us2) dx10

us1us1 us1us2 us2us1 + us2us2 dx =

10

f us1 f us2 f us1 + f us2 dx10

us1us1 2us2us1 + us2us2 dx = 0

1

0

(us1

us2)2 dx = 0

Sappiamo che us1(x) us2(x) = (us1 us2)(x) e che10

(us1 us2)2 dx = 0 se e solo se us1 us2 = 0.Ne segue che us1 us2 = 0.

Ma cosa significa che la derivata prima e zero? ne deduciamo che la differenza tra le due funzionius1 us2 e costante nellintervallo [0, 1], diciamo us1 us2 = k. Sapendo che appartengono entrambe aH10 allora us1(0) = us2(0) = 0 abbiamo che k = 0, ovvero us1 = us2 .

Soluzione variazionale e soluzione forte Supponiamo che us sia abbastanza regolare2 per soddisfare

1

0

us dx =

1

0

f dx

10

us dx =10

f dx per parti 2.5 e (0) = (1) = 0

10

us dx 10

f dx = 0

10

us + f dx = 0

10

(us + f) dx = 0

Abbiamo come conseguenza che 1

0(us + f) dx = 0 solo se u

s + f = 0 in [0, 1].

Dimostrazione 2.4.2. Assumiamo che non sia vero, quindi w = (us + f) continua in [0, 1] e w > 0 in unsotto intervallo [x1, x2] e

10

(us + f) dx = 0.Ora scegliamo tra le possibili funzioni (x) = (xx1)2 (xx2)2, continua nellintervallo e omogenea

rispetto alle condizioni al bordo.

Per lipotesi abbiamo:

0 =

10

(u + f) dx avendo (x) = 0 fuori da [x1, x2]

=

x2x1

(u + f) dx. avendo lipotesi (us + f) > 0

> x2

x1

dx essendo > 0 e > 0

> 0 contraddizione

2 quindi che esiste us e sia continuo


23/58

2.4. CONDIZIONI AL BORDO 17

Ricapitolando Abbiamo mostrato come la soluzione della formulazione forte 2.1 e anche soluzione dellaformulazione debole 2.22 e del problema di minimizzazione 2.23. Infine abbiamo visto che la soluzionedel problema variazionale 2.22 e, se sufficientemente continua, anche una soluzione per la formulazioneforte 2.1.

Osservazione 2.22. I passaggi sono strettamente relativi al problema posto allinizio 2.1 con condizioni diDirichlet. Analogamente per altri problemi possiamo intraprendere lo stesso tipo di sentiero per mostrarelequivalenza tra le tre formulazioni, relativamente allequazione differenziale e alle condizioni al bordo.

2.4 Condizioni al Bordo

In precedenza ci siamo riferiti quasi sempre a condizioni al bordo di Dirichlet; possiamo anche riferirciad altri tipi come

u(0) = 0 u(1) = 0 (di Dirichlet) (2.24)

u(0) = 0 u(1) = k2 (di miste) (2.25)

u(0) = k1 u(1) = k2 (di Neumann) (2.26)Il tipo di condizioni al bordo e fondamentale per poter creare la formulazione variazionale ed in seguito

per giustapporle in quella che verra chiamata matrice di stiffness.Mostriamo in seguito alcune varianti.

2.4.1 Alcune formulazioni variazionali

1) u = fu(0) = u(1) = 0

x I =]0, 1[ f C(I) f L2(I) u C2(I) (2.27)

Abbiamo visto (esempio 2.15) che H10 (I)10

us dx =

10

f dx (2.28)

se e vero allora e lo stesso che dire che f e la derivata debole di us. Essendo f continua ha derivata,ovvero f = (us) quindi us C2(I) ossia us e la soluzione forte.

2) u + u = fu(0) = u(1) = 0

x I =]0, 1[ f C(I) f L2(I) u C2(I) (2.29)

Per ogni H10 (I)u + u = f1

0

u + u dx = 10

f10

(u + u) dx =10

f

10

u dx +10

u dx =

10

f

u(1)(1) + u(0)(0) +10

u dx +10

u dx =

10

f 10

u dx +10

u dx =

10

f

se e vero allora e lo stesso che dire che f e la derivata debole di us. Essendo f continua ha derivata,ovvero f = (us) quindi us C2(I) ossia us e la soluzione forte. abbiamo che e equivalente a trovareu H10 (I) per


24/58


3) u + u = fu(0) = u(1) = 0

x I =]0, 1[ f C(I) f L2(I) u C2(I) (2.30)

Osservazione 2.23. Il problema e in forma simile a quello di 2), ma a differenza del precedente ci sonodifferenti condizioni al bordo.

Supponiamo che u e una soluzione forte

u + u = f10

u + u dx =10

f10

(u + u) dx =10

f

1

0

u dx + 1

0

u dx = 1

0

f

u(1)(1) + u(0)(0) +10

u dx +10

u dx =

10

f 10

u dx +10

u dx =

10

f

Abbiamo usato le condizioni al bordo u(0) = u(1) = 0,Se prendiamo il problema di trovare la soluzione forte allora equivale a trovare u H1(I) : 1

0u +

u dx =10

f .Viceversa, se u H1 soddisfa la formulazione variazionale allora u e un solo elemento della formu-

lazione forte.

10

u dx =10

(f u) dx C(I)

per definizione di derivata debole (2.18)(u) = f + u

quindi u C(I) e u C2(I).

4) u = fu(0) = u(1) = 0

x I =]0, 1[ f C(I) f L2(I) u C2(I) (2.31)

soluzione non unicau = fu(0) = u(1) = 0

x I =]0, 1[ f C(I) f L2(I) u C2(I) (2.32)

Se f = 0 allora u = 0 con u(0) = u(1) = 0 sono tutte le infinite soluzioni u(x) = k.

2.4.2 Riassumendo

Osservazione 2.24. Nelle precedente formulazioni variazionali ci sono molti punti di contatti e analogiacon cui vengono trattati i problemi.

Osservazione 2.25. Tutte le formulazioni presentato sono di equazioni differenziali di secondo ordine, oanche di ordine pari, con il funzionale B bilineare e l lineare.

Osservazione2.26. Uniformita con gli appunti e altri testi denoteremo la forma bilineare B(u, ) = a(u, )


25/58

2.5. METODO VARIAZIONALE DI APPROSSIMAZIONE 19

Problema u B(u, ) =1) H10 (I)

10

u dx2) H10 (I)

1

0u + u dx

3) H

1

(I) 10 u + u dx4) {u H1(I)|u(0) = 0} 10

u dx

Table 2.1: Spazi funzionali ove giace la soluzione al problema in forma variazionale e funzionale bilineare

2.5 Metodo variazionale di approssimazione

Lobiettivo di un metodo variazionale di approssimazione e trovare, da uno spazio infinito dimensionale,una soluzione in uno spazio finito dimensionale.

Prendiamo in generale un qualsiasi spazio dimensionale finito Vh = span(1, . . . , m) con le funzionidi base 1im(x) V.Osservazione 2.27. Lo spazio Vh e un sottospazio finito dimensionale di V, questultimo uno spazio

infinito dimensionale. Vh = span(1, . . . , m) = {vh| vh = mi=0 ii} V.Osservazione 2.28. In base alla tabella 2.1, possiamo prendere V = H10 ([0, 1]) o V = H

1([0, 1]), inrelazione al problema trattato

Al fine di aver riferimenti pratici, o meglio visivi, con cio che stiamo trattando, introduciamo ilconcetto di funzione lineare a tratti, per dedicarci in seguito al metodo di Ritz e quello di Galerkin.

2.5.1 Funzioni lineari a tratti

Prendiamo un sottospazio finito dimensionale Vh di uno spazio infinito dimensionale V, i.e. Vh V.Possiamo prendere come Vh un insieme di funzioni lineari a tratti partizionando lintervallo [0, 1] in

intervallini Ij di lunghezza hj come mostrato in figura 2.3. Piu precisamente essendo xi i punti della

x6x4x3 x5x2x1 x7 = 10 = x0

vh

Figure 2.3: La funzione rosa vh e lineare a tratti.

partizione, con i = 0, . . . , m + 1, abbiamo Ij = [xj1, xj ] e hj = xj xj1 con j = 1, . . . , m + 1.Prendiamo Vh come linsieme delle funzioni vh lineari in ogni intervallo Ij, continue in [0, 1] e con

vh(0) = vh(1) = 0.

Osservazione 2.29. Si tenga presente che la proprieta vh(0) = vh(1) = 0 e scelta appositamente da noiper trattare un problema con condizioni al bordo di Dirichlet. Si puo variare questa condizione in mododa avere un funzione lineare a tratti omogenea con il problema dato, ovvero con stesso tipo di gestionedelle condizioni al bordo (vedi figura 2.4 (B)).

Se prendiamo una funzione di base j

V con j = 1, . . . , m definita come

j(xi) :

1 i = j

0 j = i i, j = 1, . . . , m


26/58


possiamo rappresentare vh Vh come combinazione lineare delle funzioni di base i

vh

=

m

i=1

i(x) x

[0, 1]

dove e il valore vh(xi).

x3x2x10 = x0 x4 = 1

vh

x3x2x10 = x0 x4 = 1

1 12(x)

Figure 2.4: (A) funzione di base 2(x), (B) una lineare a tratti con relative funzioni di base

Abbiamo introdotto il concetto di funzione lineare a tratti. Luso di questo strumento nel metododegli elementi finiti in dimensione 1 sara concentrato nel trovare i valori di = u(xi) ricostruendo unapprossimazione uh della soluzione forte u.

2.5.2 Metodo di Ritz-Galerkin

Consideriamo il problema variazionale di trovare la soluzione u tale che

a(u, v) = l(v) (2.33)

per le funzioni v differenziabili due volte e omogenee rispetto alle condizioni al bordo poste nel problema.Quando il a e bilineare e simmetrica e l e lineare allora il problema e equivalente alla minimizzazione delfunzionale quadratico

I(u) =1

2a(u, u) l(u) (2.34)

Il metodo di Ritz e unapprosimazione dellequazione 4.17 nella forma della serie finita

un = 0 +n

j=1

cjj (2.35)


27/58


dove le costanti cj sono dette coefficienti di Ritz, scelti in modo che per v = i=1...n

a(un, v) = l(v) v = i

a(un, i) = l(i) eq. 4.18

a(0 +n

j=1

cjj , i) = l(i) esplicito

a(0 + c11 + + cnn, i) = l(i) per linearita 2.19a(0, i) + c1a(1, i) + + cna(n, i) = l(i) ovvero

a(0, i) +

nj=1

cja(j , i) = l(i) cio che e noto a destra

nj=1

cja(j , i) = l(i) a(0, i) i {1, . . . , n}

Quindi abbiamo n equazioni, una per ogni in

j=1 cja(j , 1) = l(1) a(0, 1)...n

j=1 cja(j , n) = l(i) a(0, n)esplicitando

c1a(1, 1) + + cna(n, 1) = l(1) a(0, 1)...

c1a(1, n) + + cna(n, n) = l(n) a(0, n)ovvero in forma matriciale

a(1, 1) . . . a(n, 1)... . . . ...a(1, n) . . . a(n, n)

c1...cn

= l(1) a(0, 1)...l(n) a(0, n)

A C = F (2.36)

Chiamiamo A matrice di stiffness ed F il vettore di carico3.

Osservazione 2.30. Trovare i valori di ci che soddisfano il sistema equivale a trovare i valori di u(xi) erelativa funzione lineare a tratti un 4.18. Per questo motivo talvolta scriveremo A U = F.

Il metodo di Ritz puo essere visto come trovare i ci che sono soluzione dellequazione 4.18 determi-nandoli come minimizzazione del funzionale I(u) dellequazione 2.34. Ci sono una sequenza di passaggi(vedi [Red84] pag. 37-38) i quali ci consentono mi mettere in evidenza la condizione necessaria perche ilmetodo di Ritz funzioni correttamente, ovvero che il funzionale a sia bilineare e simmetrico.

Inoltre la funzione di approssimazione i deve: essere omogenea rispetto alle condizioni al bordo essere sufficientemente differenziabile in modo da soddisfare la forma bilineare linsieme {i} sia completo (vedi 2.14).

Queste condizioni ci garantiscono la convergenza della soluzione un, ricavata applicando il metodo, rela-tivamente alla soluzione forte u, con laumentare di n.

La sezione e intitolata Metodi di Ritz-Galerkin per il seguente motivo. Il metodo di Galerkin e uncaso particolare, o meglio una delle tante varianti, del method of wieghted residuals4 o conosciuto anchecome metodo di Petrov-Galerkin. Quando lequazione differenziale si presenta di grado pari allora ci siriduce al metodo di Ritz ottenendo una matrice simmetrica.

3la lettera F e usata per la presenza della funzione nota f in l(i)4 traducibile nel metodo dei residui pesati. Si tratta nel prendere il problema posto nella forma forte prendendo un

come u; essendoci una differenza, o meglio errore, tra le due funzioni allora avremmo che E un f = 0. A questo punto

moltiplichiamo E per una funzione di residuo i la quale consente che la loro integrazione deve essereRIi(x)E(x, ci) = 0.

Il metodo di Galerkin e il caso particolare in cui i = i.


28/58


Osservazione 2.31. Rivedendo losservazione 2.25, possiamo dire che i problemi presenti nella tabella 2.1sono tutti risolubili con il metodo che abbiamo presentato.

Osservazione 2.32. Libri diversi hanno lievi differenze nel presentare il metodo, convergono tutti nel

considerare la ricerca del problema di minimizzazione di in uno spazio finito dimensionale come metododi Ritz

trovare uh Vh : F(uh) F(v) v Vhe la ricerca di una soluzione al problema posto in forma variazionale in uno spazio finito dimensionalecome metodo di Galerkin

trovare uh Vh : a(uh, v) = l(u) v VhInoltre la varieta dei problemi e delle approssimazioni consente un panorama ampio di scelta e com-

binazioni.

2.6 Stima dellerrore

2.6.1 Introduzione

Possiamo commettere errori di approssimazione del dominio (bordo), errori di integrazione numerica (for-mulazione variazionale), aritmetici (soluzione del sistema di equazioni) ed infine errori nellapprossimarela soluzione.

Ma quando ci riferiamo alla analisi o stima dellerrore e esclusivamente tra la soluzione approssimatae quella forte.

Prima di proseguire, e interessante riportare un esempio tratto da Reddy pag. 3 [Red84]

Figure 2.5: Cerchio approssimanto da 6 triangoli (A) e da 9 (B)

Esempio 2.17. Determinare larea di un cerchio di raggio r rappresentandolo come una collezione ditriangoli.

Possiamo osservare nella figura 2.5, senza entrare nel dettaglio, che allaumentare dei triangoli miglioralapprossimazione del poligono e diminuisce lerrore denotato dallarea rosa.

Con lo stesso concetto di base possiamo analizzare la bonta di un approssimazione.

2.6.2 Analisi dellerrore

Se prendiamo un equazione differenziale di cui possiamo trovare una soluzione esatta u e vi applichiamoil metodo otterremo una soluzione approssimata u

Come visto nellesempio precedente, lerrore e dato dalla soluzione esatta, o reale, e lapprossimata;

per stimarla ci sono varie vie, tra queste le piu usate in generale sono la norma, o distanza, in L2:

Definizione 2.24 (norma L2). norma L2 : u u0 =b

a|u u|2dx

12


29/58

2.7. RICAPITOLANDO 23

xn = ba = x0

Figure 2.6: Errore in 1D

Abbiamo due tipi di impiego della stima dellerrore.

Una relativa alla bonta del metodo conoscendo la soluzione esatta, ci serve per analizzare la correttezzadel software realizzato o per attuare stime sul numero di intervallini che posso dare buoni risultati.

Laltra e nella realizzazione di un metodo adattivo che ad ogni passo controlla se, rispetto alle soluzioniprecedenti, ce una riduzione sensibile dellerrore; ovvero confronta le due soluzioni approssimate a stepi differenti dove al crescere di i aumenta il numero degli intervalli, i.e. |ui+1 ui|.

Per esempio nella figura 2.6 potrebbe essere ripartizionato solo il primo intervallino per poter diminuirelerrore rappresentato dallarea rosa.

2.7 Ricapitolando

Abbiamo un equazione differenziale con condizioni al bordo da risolvere, riformuliamo il problema informa variazionale; suddividiamo lintervallo in elementi, anche in modo non uniforme, al fine di generarelequazione in forma matriciale relativa ai punti di discretizzazione. Il sistema finale risolto ci restituiscele soluzioni esatte nei punti della partizione e ci consente di delineare una approssimazione che sar a vicinaalla soluzione reale al crescere della raffinazione della partizione.


30/58



31/58

Chapter 3

FEM 1D in pratica

Includiamo in questa sezione un esempio di implementazione in maxima di un risolutore di problema.

Una sua generalizzazione per vari problemi e lasciata al lettore.

3.1 Formulazione variazionale

Scegliamo un problema, per esempio prendiamo il problema visto in 2.4.1,u = fu(0) = u(1) = 0

Prendiamo f(x) = x + 3 di cui possiamo ricavare la soluzione forte

u(x) = x3 + 9x2

6 +5x

3 (3.1)

Osservazione 3.1. Il problema sopra esposto equivale a trovare il comportamento di una corda elasticasottoposta ad una forza f, generalmente consideriamo u e non u per una questione di convenzionerelativa al fatto che la forza (tipo di gravita) si usa indicarla dallalto verso il basso.

Nel caso di una svista di non considerare il segno, si avr a un equazione simmetrica rispetto allasse x:

u(x) = +x3 + 9x2

6 5x

3

Ne ricaviamo la forma variazionale come visto in precedenza

10

us dx = 10

f dx

Osservazione 3.2. Ai fini dellimplementazione e piu comodo porla in termini di un funzionale per potergenerare la matrice di stiffness e il vettore forza da funzioni che possiamo sostituire in base al problematrattato.

a(us, ) = l()

con

a(us, ) =

10

us dx

e

l() =10

f dx

25


32/58

26 CHAPTER 3. FEM 1D IN PRATICA

3.2 Funzione di approssimazione

Possiamo partizionare lintervallo [0, 1] in vari modi come

x : SORT(APPEND(MAKELIST(RANDOM(1.0),n,1,16),[0.0,1.0]));

x : SORT(APPEND(MAKELIST(n/10,n,1,9),[0.0,1.0]));

Le parti precedenti di codice sono dei parametri che daremo in pasto alla funzione che ora andiamo acreare.

/* Il fantastico mondo dei commenti annidati di maxima

KILL(all); /* ripuliamo tutte le precedenti associazioni fatte */

LOAD("fem1d_dirichlet.mc"); /* file dove salviamo quanto segue */

f(x) :=1; xx: SORT([ 0.0,1.0,0.5,0.25,0.75]);

u: fem1d_dirichlet(f,xx); */

fem1d_dirichlet(f,lx) := BLOCK(

/* dentro a BLOCK il separatore e la virgola , */

[dist,genK,genF,x,h], /* tutte le variabili locali */RATPRINT:FALSE, /* evitiamo un output verboso */

Inoltre per controllare eventuali errori creiamo una funzione di proiezione x(i) e una funzione h(i) perottenere li-esima distanza:

N : LENGTH(x)-1 ,

x(i):=IF iN+1 THEN PRINT("x(",i,") range error") ELSE RATEXPAND(lx[i]),

h(i):=IF iN THEN PRINT("h(",i,") range error") ELSE RATEXPAND(lx(i+1)-lx(i)),

Tra le possibili funzioni di approssimazione che possiamo usare, decidiamo di usare quella trattata nellasezione 2.5.1.

x3x2x10 = x0 x4 = 1

1 1

x3x2x10 = x0 x4 = 1

u = un

2(x)

(0) (1)

(1)1

(1)2

(1)

u1(1)2

(2)1

Figure 3.1: Limmagine (A) mostra la decomposizione di in in (B) vengono poste in evidenza lefunzioni che generano un segmento ui

Se osserviamo meglio linfluenza di nella creazione della soluzione approssimata, vedi figura 4.2;

possi notare che per un singolo tratto ue relativo ad un elemento (e) vi sono due funzioni che lo generano,che chiamiamo

(e)1 e

(e)2 .

ue(x) = ce1(e1)2 + ce

(e)1 (3.2)


33/58

3.3. GENERIAMO LA MATRICE DI STIFFNESS 27

infine possiamo esplicare (e)1,2

(e)1 =

xe+1 xxe+1

xe

(e)2 =

x xexe+1

xe

(3.3)

quindi

e(x) :

(e1)2 xe1 < x < xe

1 x = xe

(e)1 xe < x < xe+1

(3.4)

Possiamo o fare delle osservazioni sulla funzione di base, affermando che i = 1/hi oppure lasciarlocalcolare allelaboratore.

/* implementate per sport, visto che usiamo direttamente 1/h(i) */

phi(x):= MAKELIST(

IF x(e)


34/58


Osservazione 3.4. Nel codice compaiono p[1] e p[2] che rappresentano uno strumento per gestire ilcomportamento delle derivate negli estremi.

Ai fini della soluzione di trovare un di questo problema non sono importanti e vengono in molti libriomesse, focalizzando il risultato solo nelle equazioni dove compaiono le incognite che ci interessano.

Al fine di creare una applicazione flessibile che tiene conto di tutte le condizioni puo dare soddisfazioneprenderle in considerazione sin da subito.

3.5 Risolvere il sistema lineare

Ora abbiamo tutto per poter risolvere il sistema lineare, potremmo risolverlo nello stile di Matlab, oOctave, usando U: K-1 . F, oppure U: INVERSE(K) . F,.

Preferiamo limitare lerrore numerico e risparmiare tempo di calcolo usando la funzione SOLVE:

/* creiamo la lista di incognite $u_i$ o $c_i$ */

list_u : MAKELIST(u[i],i,1,N+1),

/* creiamo la lista di equazioni $K.U=F$ */

list_eq : MAKELIST(K[i].list_u=F[i],i,1,N+1),

/* applichiamo le condizioni al bordo */

list_u[1] : p[1], /* incognita u(0)=p_1 */

list_u[LENGTH(list_u)] : p[2], /* incognita u(1)=p_2 */

list_eq:EV(list_eq,u[1]=0), /* u(0)=0 */

list_eq:EV(list_eq,u[LENGTH(list_u)]=0), /* u(1)=0 */

/* finalmente solve KU=F */sol:SOLVE(list_eq,list_u),

/* creiamo il vettore delle soluzioni relative alla partizione list_x */

u_tilde : MAKELIST(IF i=1 OR i>N THEN 0 ELSE RHS(sol[1][i]) ,i,1,N+1),

RETURN(FLOAT(u_tilde) /* e lo restituiamo al chiamante * /

) /* chiusura di BLOCK( */

3.6 Un esecuzione

Avviamo maxima caricando la funzione che abbiamo scritto, scegliamo una funzione f(x) = x + 3 epartizioniamo in 11 punti equidistanti.

KILL(all); LOAD("fem1d_dirichlet.mc");

f(x):=x+3; x:SORT(APPEND(MAKELIST(n/10,n,1,9),[0.0,1.0])); c:fem1d_dirichlet(f,x);

TEX(u); FPPREC: 4; TEX(K);

otterremo come risultato

c : [0.0, 0.1515, 0.272, 0.3605, 0.416, 0.4375, 0.424, 0.3745, 0.288, 0.1635, 0.0]

Tratto dalla matrice di stiffness che osserviamo essere trilineare e simmetrica come ci aspettavamo:


35/58

3.7. STIMA DELLERRORE 29

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0 0.2 0.4 0.6 0.8 1

fem

u_nu

Figure 3.2: Funzione un lineare a tratti rossa e u soluzione reale in verde, con 11 nodi

K :

10.0 10.0 0 0 0 0 0 0 0 0 010.0 20.0 10.0 0 0 0 0 0 0 0 0

0 10.0 20.0 10.0 0 0 0 0 0 0 00 0 10.0 20.0 10.0 0 0 0 0 0 00 0 0 10.0 20.0 10.0 0 0 0 0 00 0 0 0 10.0 20.0 10.0 0 0 0 00 0 0 0 0

10.0 20.0

10.0 0 0 0

0 0 0 0 0 0 10.0 20.0 10.0 0 00 0 0 0 0 0 0 10.0 20.0 10.0 00 0 0 0 0 0 0 0 10.0 20.0 10.00 0 0 0 0 0 0 0 0 10.0 10.0

F : [p1 + 0.1517, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, p2 + 0.1983]t

Se raffinando ulteriormente la partizione con 21 elementi

un : [0.0, 0.0795625, 0.1515, 0.2156875, 0.272, 0.3203125, 0.3605, 0.3924375, 0.416, 0.4310625, 0.4375,

0.4351875, 0.424, 0.4038125, 0.3745, 0.3359375, 0.288, 0.2305625, 0.1635, 0.0866875, 0.0]

otteniamo maggiore vicinanza alla soluzione forte (eq. 3.1), come si puo direttamente osservare in figura3.3.

3.7 Stima dellerrore

Volendo si possono automatizzare le stime degli errori, per molte equazioni, creando un modulo basatosulla seguente intuizione.

phi(x,i):= /* gestirla */

c:u,

u_tilde(x) := phi(x,0) + sum(c[i]*phi(x,i),i,2,N),

/* u reale */

de : -diff(y,xi,2)=f(xi),

gsol : ode2(de,y,xi),

bc1 : EV(gsol, [y=0, xi = 0]), /* u(0.0)=0 */


36/58


0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0 0.2 0.4 0.6 0.8 1

fem

u_nu

Figure 3.3: Funzione un lineare a tratti rossa e u soluzione reale in verde, con 21 nodi

bc2 : EV(gsol, [y=0, xi = 1]), /* u(1.0)=0 */

csol : solve([bc1, bc2], [%K1, %K2]),

u_real(xi):= RHS(EV(gsol, csol)),

/* con unopportuna phi converge */

err:SQRT(ROMBERG( (u_real(xi)-u_tilde(xi))2 , xi, 0.0, 1.0)),

PRINT("Errore in L_2 : ", err),

3.8 Consigli per lo sviluppo

Come visto nel corso di Analisi Numerica, un implementazione di un algoritmo puo essere soggetta aderrori di diversa natura. Gli errori numerici piu classici fatti nella progettazione del metodo sono iseguenti:

Conversioni da classi di numeri ad altre, come float e double, con linguaggi come C++. Integrazione numerica. Esistono vari metodi e parametri per rendere minimo lerrore.

Risolvere un sistema di equazione con luso di una matrice inversaA U = F U = A1 F

al posto di usare il metodi di fattorizzazione LUo QR od altre funzioni stabili rese disponibili dallapiattaforma di sviluppo.

Osservazione 3.5. La rappresentazione su un monitor attuale avviene mediante luso dei pixel, i qualipermettono una rappresentazione discreta delle nostre funzioni. Se raffiniamo molto la nostra partizionea video, la funzione reale e la poligonale ottenuta applicando il metodo degli elementi finiti sarannosovrapposte, sebbene non siano identiche.


37/58

Chapter 4

Fondamenti teorici di FEM 2D

Nel seguente capitolo vediamo dei concetti su funzioni a due o piu variabili, per trattare in seguito la

formulazione variazionale di un problema e la sua approssimazione nel metodo visto in una dimensione.

4.1 Richiami e complementi di analisi

Ora richiamiamo e presentiamo concetti relativi a funzioni a due variabili reali o anche definite in sot-toinsiemi di R2

f : A R2 Re funzioni a piu variabili reali.

Come notazione useremo x = (x1, . . . , xn) Rn.

4.1.1 Insieme aperto

Definizione 4.1 (Boccia). Dato uno spazio vettoriale V e una distanza d, quindi (X, d) e uno spaziometrico, x V, r R+. Diciamo boccia aperta, o sfera aperta, linsieme dei punti x1 compresi entro unraggio r da un centro x

B(x, r) = {x1 V : d(x, x1) < r}

Figure 4.1: (A),(B) e (C) rappresentano nelle varie metriche B1(0, r), B2(0, r) e B(0, r)

Definizione 4.2 (Insieme aperto). A R2

si dice aperto se x A B(x, ) A : B(x, ) = {(y1, y2) R2|(y1 x1)2 + (y2 x2)2 }Definizione 4.3 (insieme chiuso). Un insieme A e detto chiuso se e il complementare di uno aperto A.

31


38/58


Definizione 4.4 (punto interno). Un punto x A si dice interno ad un insieme E A se esisteB(x, > 0) contenuta in E. Indichiamo linsieme dei punti interni ad E con Int(E).

Definizione 4.5 (punto esterno). Un punto x A si dice esterno rispetto ad un insieme E A seappartiene a Int(E).Definizione 4.6 (punto di frontiera). Un punto si dice di frontiera se per ogni > 0 se e solo se x A:

B(x, ) B(x, ) E= B(x, ) E= Indichiamo con F r(E), o E, linsieme dei punti di frontiera.

Osservazione 4.1. Un insieme A e chiuso se e solo se F r(A) A.Esempio 4.1. Prendiamo come esempio un ellisse

A = {(x, y)|x2 + 2 y2 3}Avremo

F r(A) = {(x, y)|x2 + 2 y2 = 3}Int(A) = {(x, y)|x2 + 2 y2 < 3}

A = {(x, y)|x2 + 2 y2 > 3}Definizione 4.7 (insieme limitato). Un insieme A R2 e detto limitato se e contenuto da una bocciaB(0, r).

Definizione 4.8 (insieme compatto). Quando un insieme e chiuso e limitato lo chiamiamo insiemecompatto.

Definizione 4.9 (insieme convesso). Quando la congiunzione di due punti qualsiasi di un insieme nonesce dallinsieme stesso, parliamo di un insieme convesso.

Definizione 4.10 (insieme connesso). Se esiste una spezzata a tratti tra due punti qualsiasi di un insieme,tale che ogni segmento non esca dallinsieme stesso, allora parliamo di un insieme connesso.

BA

Figure 4.2: Linsieme B e un insieme connesso, mentre A non lo e.

4.1.2 Derivazione e continuita

Riprendiamo il concetto di continuita su R.

Definizione 4.11. Diciamo che una funzione f : I R e continua su x I se e solo se limxx f(x) =f(x) ovvero > 0 > 0 : x I |x x| = |f(x) f(x)| < Definizione 4.12. In R2 parliamo di continuita se data f : A R2 R, con A aperto, abbiamo chelimx

x f(x) = f(x) ovvero

> 0

> 0 :

x

B(x, ) f(x)

f(x)

|<

Definizione 4.13. Sia A un aperto e f : A R2 R, diciamo che f C1(A) se esistono e sono continuef

x1(x),

f

x2(x) x A.


39/58


Abbiamo che se f C1(A) allora fv

(x) per ogni versore v e per ogni punto x A con

limh0 f(x + hv) f(x)h = fv(x) = f(x), v (4.1)

con

=

x1, . . . ,

xn

tf(x) =

f(x, y)

x,

f(x, y)

y

t(4.2)

v = (1, 0) f(x), v = f(x, y)x

(x, y) v = (0, 1) f(x), v = f(x, y)y

(x, y)

Esempio 4.2. Studiamo il comportamento dif

v(0, 0) e con il versore v = (

2

2,

2

2) di

f(x, y) = sin(x + y) + 3x2 (4.3)

f(x, y) = x cos y + exy (4.4)

f(x, y) = ex2y2 (4.5)

Prendiamo 4.3, abbiamo due modi per muoverci. In uno, ne osserviamo il limite come mostrato nel

primo elemento dellequazione 4.1:

f

v(x) = lim

t0f(x + tv) f(x)

t

= limt0

f((0, 0) + (t22

, t22

)) f(0, 0)t

= limt0

sin((0 + t22

) + (0 + t22

)) + 3(0 + t22

)2

t

= limt0

sin(2 t 22

) + 3(t22

)2

t

= limt0sin(t

2) + 3 t224t

= limt0

sin(t 2) + 32

t2

t

Abbiamo che per Taylor

sin(t 2)t

=t 2

t+ O(t) =

2 + O(1) =

2

quindi

limt0

sin(t 2) + 32

t2

t= lim

t0

2 +32

t2

t

2


40/58


Laltra forma di risoluzione e dato usando lultima parte dellequazione 4.1:

f

v(x) = f(x), v

= f(0, 0), (22

,2

2) per eq. 4.2

=

f(0, 0)

x,

f(0, 0)

y

t, (

2

2,

2

2)

= ((cos(x + y) + 6x)(0,0), (cos(x + y))(0,0))t, (

2

2,

2

2)

= (cos(0 + 0) + 6 0, cos(0 + 0))t, (

2

2,

2

2)

= (1, 1), (

2

2,

2

2)

= 22

+ 22

= 222

= 2

Prendiamo 4.4 usando solo il metodo delle derivate parziali

f

v(x) =

f(0, 0)

x,

f(0, 0)

y

t, (

2

2,

2

2)

= ((1cos y + y exy)(0,0), (x( sin y) + x exy)(0,0))t, (

2

2,

2

2)

= (0, 0)t, (

2

2,

2

2)

= 0

Prendiamo 4.5 usando solo il metodo delle derivate parziali

f

v(x) =

f(0, 0)

x,

f(0, 0)

y

t, (

2

2,

2

2)

= ((2x ex2y2)(0,0), (2y ex2y2)(0,0))t, (

2

2,

2

2)

= (0, 0)t, (

2

2,

2

2)

= 0

Esempio 4.3. f(x, y) = sin y + 3x + exy studiaref

v(x) con x = (1, 0) e v = (

3

2, 12

)

f

v(x) =

f(1, 0)

x,

f(1, 0)

y

t, (

3

2,

1

2)

= ((3 + y exy)(1,0), (cos y + x exy)(1,0))t, (

3

2,

1

2)

= ((3, 1)t, (

3

2,

1

2)

= 3

3

2+

1

2=

3

3 + 2

2

Definizione 4.14. Sia f : A R2 R, diciamo che f C2 se esistono e sono continue2f

x2,2f

xy,

2f

yx,

2f

y2, in ogni punto di A


41/58


Teorema 4.1 (Schwartz). Dato un aperto A Rn e una funzione f : A R abbiamo2f

xixj(x) =

2f

xjxi(x)

Definizione 4.15 (Hessiana). Sia A un aperto di Rn e f C2, chiamiamo Hessiana la seguente matrice:

Hessf(x) =

2fx1x1

. . . 2f

xnx1...

. . ....

2fx1xn

. . . 2f

xnxn

Teorema 4.2 (Formula di Taylor). Per f(x) C2(]a, b[) e x ]a, b[

f(x) = f(x) + f(x) +f(x)

2(x x)2 + o(|x x|2)

Per un aperto A Rn

conx, x A e f C2

allora

f(x) = f(x) + f(x), x x + 12

Hessf(x)(x x), (x x) + o(

(x x)2 + (x x)2)

4.1.3 Curve

Definizione 4.16 (parametrizzazione). Sia I = [a, b], : I R2 con = (1.2), diciamo che e unaparametrizzazione se

1. C1(I) con 1, 2 che sono in C1

2. (t) = 0 t IEsempio 4.4. Possiamo rappresentare una circonferenza con una funzione parametrica (t) = (cos t, sin t)con t [0, 2].

In gnuplot e possibile visualizzarla con set parametric; set trange [-pi:pi]; plot sin(t),cos(t);

Definizione 4.17. Se I = [a, b] e J = [, ] con : I R2 e : J R2, allora se f : I Jcon f C1(I) iniettiva tale che

1. f(I) = J e f(t) > 0 t2. (t) = (f(t)) t

Esempio 4.5. Le seguenti curve sono equivalenti:

(t) = (s, 1 s2) con s [1, 1] (t) = (cos t, sin t) con t [0, ]

In gnuplot e possibile visualizzarla con set parametric; set trange [0:pi] ; plot cos(t),sin(t)e set parametric; set trange [-pi:pi]; plot -t, sqrt(1-t**2)

Definizione 4.18. Chiamiamo F linsieme delle parametrizzazioni C1. Con F/ indichiamo la classedelle relazioni di equivalenza.

Definizione 4.19. Definiamo curva regolare di classe C1 (orientata) un elemento del la classe di equiv-alenza F/.Definizione 4.20. Se = [] F/ , con : I Rn diremo che e una parametrizzazione del lacurva .Definizione 4.21 (curva semplice). Se abbiamo una parametrizzazione e iniettiva di una curva, questae detta semplice


42/58


Definizione 4.22 (curva semplice aperta). Una curva semplice parametrizzata da : [a, b] = I Rne (a) = (b) allora la curva e detta curva semplice aperta.Definizione 4.23 (curva semplice chiusa). Una curva semplice parametrizzata da : [a, b] = I Rne (a) = (b) allora la curva e detta curva semplice chiusa.

Definizione 4.24. Sia g : A R2 R con g C(A) e sia = [] F/ allora

g ds =

ba

g((t))(t) dt

=

ba

g((t))(

1(t)2 + 2(t)

2 dt

Osservazione4.2. In particolare con g(x, y) = 1 e () =

1 ds =ba

1(t)2 +

2(t)

2 dt con

1(t)2 + 2(t)

2

norma.

Esempio 4.6. Con g(x, y) = x + y e (cos t, sin t) allora

g ds =

20

g(cos t, sin t)(

( sin)2(t) + cos2(t) dt

=

20

g(cos t, sin t)(

1 dt

=

20

cos t + sin t dt

=

20

cos t dt +

20

sin t dt

= [ sin t]20 + [cos t]20= sin2 sin 0 + cos 2 cos0= 0 0 + 1 1 = 0

Osservazione 4.3. Si osservi la lunghezza della curva del precedente esempio e20

( sin)2(t) + cos2(t) dt = 2

in altre parole, quando g = 1 abbiamo la lunghezza della curva. inoltre () = ()Osservazione 4.4. Se abbiamo due parametrizzazioni equivalenti (t) = (f(t)) allora si osservi che

g ds = b

a

g((t))(1(t)

2 + 2(t)2 dt (t) = (f(t)) (t) = (f(t)) f(t)

=

ba

g((f(t)))(

1(f(t))2 + 2(f(t))

2 f(t) dt f(t) = s ds = f(t)dt

=

g((s))(

1(s)2 + 2(s)

2 ds

4.1.4 Integrazione

Riprendiamo lintegrazione secondo Riemann per funzioni con una variabile.Definiamo una partizione di un itervallo I = [a, b] come

= {xi R | a = x1 < x2 < < xn = b}

Somme inferiori

s(f, ) =n

j=1

infxIj

f(x)(xj xj1)


43/58


xi1 xi

f(x)

a = x1 xn = b

Figure 4.3: Integrazione di f(x) secondo Riemann, in azzurro i supxIj in rosa gli infxIj

Somme superiori

S(f, ) =n

j=1

supxIj

f(x)(xj xj1)

Definizione 4.25 (Riemann integrabile 1D). Diciamo che f e Riemann integrabile su [a, b], e scriveremo

f R([a, b]), se sup s(f, ) = infS(f, ) =

b

af(x) dx o anche se > 0 : S(f, ) s(f, ) <

Teorema 4.3. In un intervallo I si ha C(I) R(I).

d

c

a

b

f(x, y)

Figure 4.4: Integrazione di f(x, y) secondo Riemann, in azzurro i sup(x,y)IjJi in rosa gli inf(x,y)IjJi

Nel caso di funzioni a due variabili f : [a, b] [c, d] R con [a.b] partizionata come

x = {xi R | a = x1 < x2 < < xn = b}[c, d] partizionata come

y = {yi R | a = y1 < y2 < < yp = b}


44/58


e = x y.Somme inferiori

s(f, ) =n

j=1p

i=1inf

(x,y)

Ij

Ji

f(x, y) (xj xj1) (yi yi1)

Somme superiori

S(f, ) =n

j=1

pi=1

sup(x,y)IjJi

f(x, y) (xj xj1) (yi yi1)

Definizione 4.26 (Riemann integrabile 2D). Diciamo che f e Riemann integrabile su [a, b] [c, d],e scriveremo f R([a, b] [c, d]), se sup s(f, ) = infS(f, ) =

[a,b][c,d] f(x, y) dx dy o anche se > 0 : S(f, ) s(f, ) <

A

Figure 4.5: In rosa abbiamo caratterizzato la funzione f in [a, b] [c, d], mentre la superficie A e ristrettaalla solo proiezione sul piano di f

Definizione 4.27. Se A R2 con A [a, b] [c, d] limitato e una funzione f : A R allora possiamocrearci una funzione

f(x, y) = f(x, y) (x, y) A0 (x, y) [a, b] [c, d] \ Atale che f C([a, b] [c, d]) e

A

f(x, y) dx dy =

[a,b][c,d]

f(x, y) dx dy

vedi figura 4.5. Diremo che f R(A) se f R([a, b] [c, d])Definizione 4.28 (x-semplice e y-semplice). Dato un insieme A R2 si dira x-semplice se e del tipo

A = {(x, y) R2|x [a, b].g1(x) y g2(x)}con g1, g2 : [a, b] R continue.

Dato un insieme A R2

si dira y-semplice se e del tipoA = {(x, y) R2|y [c, d].h1(x) x h2(x)}

con h1, h2 : [c, d] R continue.Definizione 4.29 (insieme semplice). Se A e o x-semplice o y-semplice allora e detto semplice.

Definizione 4.30 (insieme regolare). Un insieme che e unione di insiemi semplici e detto regolare.

Teorema 4.4. Se A R2 e x-semplice e B A conf C(B) allora f R eA

f(x, y) dx dy =

ba

g2(x)g1(x)

f(x, y) dy

dx

Se A

R2 e y-semplice e B

A conf

C(B) allora f

Re

A

f(x, y) dx dy =

dc

h2(x)h1(x)

f(x, y) dx

dy


45/58


Esercizio 4.1. Dato A = {(x, y)|(x 1)2 + y2 1, y2 x2} calcoliamo A 1 dx dy.Osservazione 4.5. dobbiamo calcolare un area di un cerchio di raggio 1 e centro (1, 0) intersecata conla parabola y2 x3 che equivale a y x3/2.Vogliamo gli y

x3/2 e che sia dentro la circonferenza. conosciamo il valore di 3/4 di circonferenza00000000000000000000000000000000000000000000000000000000000000000000001111111111111111111111111111111111111111111111111111111111111111111111B A B

Figure 4.6: Rappresentazione del problema

(A B) essere 3/4 e cercheremo di calcolare B.Ricaviamo il punto di intersezione

(x 1)2 + y2 = 1y2 = x3

(x 1)2 + y2 = 1y = x3/2

(x 1)2 + (x3/2)2 = 1y = x3/2

x2 2x + 1 + x3 = 1y = x3/2

x(x2 + x 2) + 1 = 1y = x3/2

x(x2 + x 2) = 0y = x3/2

x(x + 2)(x 1) = 0 vera per x = 0, x = 2 , x = 1; accettiamo la positiva x = 1 sostituiamo a y2 = x3 eotteniamo il punto di intersezione (1, 1), vedi figura 4.6.

Osservazione 4.6. Abbiamo quindi che larea di A = 3/4 +B con B = {(x, y)|0 x 1, 0 y < x3/2}

B = {(x, y)|0 x 1, 0 y < x3/2}

=

10

(

x3/20

1 dy) dx

=10

x3/2 dx =2

5x2510

= 25

abbiamo quindi che A = (A B) + B = 3/4 + 25Esempio 4.7. Prendiamo come sottinsieme di R2 linsieme A = {(x, y)|x2 + y2 4y + 2x + 4 0} m(A) = A 1 dx dy

Svolgiamo il completamento dei quadrati ottenendo:

x2 + y2 4y + 2x + 4 = (x2 + 2x) + (y2 4y) + 4= (x2 + 2x + 1) 1 + (y2 4y + 4) 4 + 4= (x2 + 2x + 1) + (y2 4y + 4) 1

= (x + 1)

2

+ (y

2

2)2

1Abbiamo (x + 1)2 + (y2 2)2 = 1 che sappiamo essere un cerchio con centro (1, 2) e raggio 1Deduciamo che larea e /2 senza usare lintegrale


46/58


4.1.5 Formula di Green

Lequivalente del teorema per parti in Rn

Teorema 4.5 (Formula di Green). Sia un insieme regolare diR2 e siano v

C1(R2) e u

C2(R2).

Allora

u(x, y)v(x, y) dx dy +

u(x, y), v(x, y) dx dy =

v un

ds (4.6)

conn il versore normale esterno al bordo di A denotato da .

Il simbolo e detto Laplaciano e equivale

u(x, y) =2u

x2+

2u

y2

da riguardare e vedere se porla dopo la divergenza

4.1.6 Teorema della divergenza e del gradiente

Un equivalente del teorema per parti inTeorema 4.6 (della divergenza).

div(F) dx dy =

F dx dy =

nF ds (4.7)

dove n denota la normale unitaria direzionale sul bordo di , denotato da , e loperatore

=

x,

y

t

4.2 Formulazione Variazionale

Definizione 4.31 (formulazione forte). Dato un aperto regolare e una funzione f C(), chiamiamolequazione di Poisson: u = f

u = 0 (4.8)Trovare u C(2)() e detta formulazione forte.Osservazione 4.7. Tranne nel caso in cui sia una circonferenza, non possiamo scrivere una soluzioneesplicita.

Nello stesso stile della formulazione variazionale in 1D data una qualsiasi funzione (x, y) che siacontinua in con x e

y continue a tratti in e (x, y) = 0 sul bordo .

2u(x, y)

x2 2u(x, y)

y2 = f(x, y)u = f moltiplicando per

u = f

u =

f

u =

f per la formula di Green 4.6

u

n

u,

=

f = 0 sul bordo

u, =

f

Abbiamo quindi ricavato

a(u, ) =

u, l() =

f (4.9)


47/58

4.3. ALCUNE FORMULAZIONI VARIAZIONALI 41

Definizione 4.32 (formulazione variazionale). Trovare u H10 () tale che soddisfi a(u, ) = a() H10 () e detta formulazione variazionale.

Osservazione 4.8. Analogamente al caso 1D possiamo dimostrare che una soluzione della formulazione

forte e anche una soluzione debole, ed il viceversa se abbiamo u abbastanza regolare.

4.3 Alcune formulazioni variazionali

Analogamente al caso monodimensionale abbiamo alcuni problemi che possiamo trattare.

4.3.1 Equazione di Poisson

Abbiamo gia visto lequazione di Poisson:

u = fu = 0 (4.10)

la cui formulazione variazionale e data da

2u(x, y)

x2

2u(x, y)

y2= f(x, y)

u = f moltiplicando per u = f

u =


u

n

u,

=

f = 0 sul bordo

u, =

f


a(u, ) =

u, l() =

f (4.11)

4.3.2 Variante

Troviamo la formulazione variazionale di

u + u = f

u = 0 (4.12)abbiamo che

2u(x, y)

x2

2u(x, y)

y2+ u = f(x, y)

u + u = f moltiplicando per u + u =

f

u +

u =


u

n u, + u = f = 0 sul bordo

u, +

u =

f


48/58



a(u, ) =

u, +

u l() =

f (4.13)

4.3.3 Problema di Neumann

Troviamo la formulazione variazionale di u + u = fun = g (4.14)

abbiamo che

2u(x, y)

x2

2u(x, y)

y2

+ u = f(x, y)

u + u = f moltiplicando per u + u =

f

u +

u =


u

n

u,

+

u =

f un

= g sul bordo

u, +

u =

f +

g


a(u, ) =

u, +

u l() =

f +

g (4.15)

4.3.4 Tabella riassuntiva

Nelle precedente formulazioni variazionali ci sono molti punti di contatti e analogia con cui vengonotrattati i problemi.

Problema u a(u, ) =Poisson H10 (I)

u, Variante H1

0(I) u, + u Neumann H1(I)

u, +

u

Table 4.1: Spazi funzionali ove giace la soluzione al problema in forma variazionale e funzionale bilineare

4.4 Metodo variazionale di approssimazione

Prendiamo in generale un qualsiasi spazio dimensionale finito Vh = span(1, . . . , m) con le funzioni dibase 1im(x) V.Osservazione 4.9. Lo spazio Vh e un sottospazio finito dimensionale di V, questultimo uno spazio infinitodimensionale. Vh = span(1, . . . , m) = {vh| vh =

mi=0 ii} V.

Osservazione 4.10. Per il problema trattato V = H10 ().

Seguendo le orme di quanto presentato in 1D vediamo lanalogo della funzione interpolante in 2D, perdedicarci in seguito al metodo di Ritz e quello di Galerkin per il problema posto.


49/58


4.4.1 Triangolazione

Prendiamo un sottospazio finito dimensionale Vh di uno spazio infinito dimensionale V, i.e. Vh V.Possiamo prendere come Vh un insieme di triangoli K, partizionando laperto regolare . Piu pre-

cisamente, chiamiamo triangolazione Tn = {K1, . . . , K n} con Ki Kj = , e denotiamo i nodi internicon N1, . . . , N m. Otteniamo come unione di triangoli Kj Tn

=n

i=1

Ki = K1 Kn

Possiamo vedere alcuni esempi come mostrato in figura 4.7.

0

0.2

0.4

0.6

0.8

1

0 0.2 0.4 0.6 0.8 1

ciclo3

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

-1 -0.5 0 0.5 1

cerchio-10-8

0

0.2

0.4

0.6

0.8

1

0 0.2 0.4 0.6 0.8 1

grid-3

Figure 4.7: Alcune triangolazioni (tratta da [BBM05])

Per descrivere tutte le funzioni che entrano in gioco dobbiamo limitarci a considerare un solo elementoKi e. Se interpoliamo la funzione con lo spazio {1, x , y}

u(x, y) = c1 + c2x + c3y (4.16)

in un singolo triangolo allora abbiamo tre punti (xi, yi) i = 1, 2, 3. Possiamo ricavare u(xi, yi) = ui,vedi figura 4.8. Questo ci mette nelle condizioni di dover determinare le costanti c1, c2, c3 in un sistema

1 1

2

31

1 1

2 2 2 2

33 3 3

u(e)1u(

e

)3

u(e)2

u(e)(x, y)

Figure 4.8: In rosa le singole (e) relative ad un elemento e, tutte le (e) generano (con opportune costantiche troveremo) la futura funzione u relativa ad e (u(e)) denotata dal triangolo blu; in azzurro vi sono le (e) che generano u(e)

di equazioni

u1 = u(x1, y1) = c1 + c2x1 + c3y1

u2 = u(x2, y2) = c1 + c2x2 + c3y2

u3 = u(x3, y3) = c1 + c2x3 + c3y3u1u2

u3

=

1 x1 y11 x2 y2

1 x3 y3

c1c2

c3

Risolvendo lequazione troviamo i ci i quali sostituiti allequazione 4.16 ci permette di trovare (e)i

(e)1 =

1

2Ae((x2y3 x3y2) + (y2 y3)x + (x3 x2)y)

(e)2 = 12Ae((x3y1 x1y3) + (y3 y1)x + (x1 x3)y)

(e)3 =

1

2Ae((x1y2 x2y1) + (y1 y2)x + (x2 x1)y)


50/58


Con 2Ae area del triangolo e

2Ae =

1 x1 y11 x2 y21 x3 y3

Osservazione 4.11. La numerazione e in senso antiorario, altrimenti il valore ottenuto e opposto di segno1.

Osservazione 4.12. Abbiamo che

i(xj, yj) = i(Nj) = ij

1 1 = j

0 i = j

Inoltre3

i=1 i(Nj) = 1.

Possiamo rappresentare vh Vh come combinazione lineare delle funzioni di base i

vh =m

i=1 i(x, y) (x, y)

dove e il valore vh(xi, yi) = vh(Ni).Abbiamo introdotto una delle possibili discretizzazioni di . Ce anche la possibilita di usare altre

forme d

elementi finiti

Documents

Transcript of elementi finiti