testo - Dipartimento di Matematica

PROBLEMA 1
LABORATORIO di CALCOLO
Laurea Magistrale in Matematica Pura ed Applicata
Facolt`
a di SMFN
Considerare il problema:
(
−αu00 + βu0 + γu = f, a < x < b
u(a) = 0,
u(b) = 0
con f ∈ L2 (a, b), α > 0, γ ≥ 0, β ∈ IR . La forma debole del problema pu`o essere scritta come
trovare u ∈ V : a(u, v) = F (v), ∀v ∈ V
dove F : V → IR funzionale lineare continuo e a(u, v) : V × V → IR `e la forma bilineare
Z b
a(u, v) =
αu0 v 0 dx +
Z b
a
βu0 vdx +
a
Z b
γuvdx
a
e V = H01 (a, b) `e un sottospazio dello spazio di Sobolev H 1 (a, b).
Consideriamo una partizione su (a, b) :
a = x0 < x1 < · · · < xN +1 = b
con
xi = a + ih,
i = 0, ..., N + 1, h = (b − a)/(N + 1)
e consideriamo l’approssimazione ad elementi finiti del problema sopra:
trovare uh ∈ Vh : a(uh , vh ) = F (vh ), ∀vh ∈ Vh
(1)
dove Vh `e un sottospazio di V a dimensione finita. In particolare consideriamo
Vh0 = {vh ∈ Xh1 : vh (a) = vh (b) = 0}.
con
Xh1 = {vh ∈ C 0 (Ω) : vh |[xi−1 ,xi ] ∈ IP1 , ∀i = 1, 2, . . . , N + 1}.
data una base per Xh1 : {ϕi (x); i = 0, 1, . . . , N + 1} il problema 1, si pu`o scrivere in forma di
sitema lineare:
Au = f
dove
A = αK + βH + γM
con
Z b
(K)ij :=
a
ϕ0j (x)ϕ0i (x)
Z b
dx,
(H)ij :=
a
ϕ0j (x)ϕi (x)
1
Z b
dx,
(M )ij :=
ϕj (x)ϕi (x) dx.
a
IMPLEMENTAZIONE MATLAB
Calcolo di A
1.a) Scrivere una function per creare una mesh uniforme su (a, b) :
a = x0 < x1 < · · · < xN +1 = b
con h = xi − xi−1 per i = 1, 2, . . . , N + 1. La function deve prendere in entrata i valori a, b
e N e restituire un vettore mesh che contiene i valori dei nodi {xi , i = 0, 1, . . . , N + 1} e il
valore di h.
1.b) Scrivere una function che dato h calcoli le matrici locali Ki , Hi , Mi , i = 1, . . . , N + 1
definite da
 R xi

R xi
0
2
0
0
xi−1 (ϕi−1 (x)) dx
xi−1 ϕi (x)ϕi−1 (x) dx


Ki =  R

R
xi
xi
0 (x)ϕ0 (x) dx
0 (x))2 dx
ϕ
(ϕ
i
i
xi−1 i−1
xi−1
 R xi
xi−1
R xi
ϕ0i−1 (x)ϕi−1 (x) dx
Hi =  R
xi
xi−1
ϕ0i (x)ϕi−1 (x) dx

xi−1

ϕ0i−1 (x)ϕi (x) dx
R xi
2
xi−1 (ϕi−1 (x))
Mi =  R
xi
R xi
dx
xi−1
R xi
xi−1
ϕ0i (x)ϕi (x) dx
ϕi (x)ϕi−1 (x) dx

xi−1
ϕi−1 (x)ϕi (x) dx
R xi
2
xi−1 (ϕi (x))
dx






considerando la definzione esplicita della base lagrangiana per Xh1 ,
ϕj (xi ) = δij
.
1.c) Utilizzare le matrici locali definite al punto precedente e scrivere una funzione MATLAB
per costruire le matrici K, H, M corrispondenti, rispettivamente, alla approssimazione di
Galerkin dei termini:
Z b
0 0
Z b
u v dx,
a
0
Z b
uv dx,
u v dx,
a
a
definire quindi
Ae = αK + βH + γM
ed estrarre da Ae la sottomatrice A corrispondente ai gradi di libert`a del problema.
1.d) Analizzare per ogni N + 1 = 10, 20, 40, 80 la struttura della matrice A (con il comando
spy) e determinare i valori h e κ2 (A).
Quale relazione intercorre fra h e κ2 (A)? Riportarli in un grafico cartesiano.
(Si ricorda che κ2 (A) indica il numero di condizionamento in norma 2 e per matrici
simmetriche definite positive risulta
κ2 (A) =
.
2
λmax (A)
λmin (A)
Calcolo di f e uh (x)
1.e) Scrivere una function per calcolare i valori locali fi , i = 1, . . . , N + 1 definiti da:
 R xi
xi−1
f (x)ϕi−1 (x) dx
fi =  R
xi

xi−1
f (x)ϕi (x) dx



dove f (x) `e una funzione integrabile su (a, b).
Per calcolare gli integrali sopra, considerare il metodo dei trapezi che permette di approssimare l’integrale di un funzione data g su un intervallo (xa , xb ) con
Z xb
xa
g(x) dx '
xb − xa
(g(xa ) + g(xb )) .
2
1.f ) Utilizzando i vettori fi definiti al punto precedente, scrivere una funzione MATLAB per
costruire il vettore f corrispondente all’approssimazione di Galerkin agli elementi finiti del
termine
Z
b
f v dx.
a
1.g) Scrivere una funzione MATLAB che data una funzione v ∈ L2 (a, b) approssimi
kvk2L2 (a,b) =
Z b
(v)2 dx
a
utilizzando la formula dei trapezi composita su una mesh pi`
u fitta di quella usata precedentemente.
1.h) Considerare adesso f (x) in modo che la soluzione del problema suddetto sia u(x) =
xsin(2πx) per (a, b) = (0, 1) e risolvere il sistema lineare Au = f associato all’approssimazione
di Galerkin agli elementi finiti utilizzando lo spazio Vh con N + 1 = 10 e tracciare il grafico
della soluzione approssimata uh .
1.i) Calcolare la norma L2 dell’errore u − uh in corrispondenza di N + 1 = 10, 20, 40, 80, per
ciascuno dei valori di N suddetti riportare in uno stesso grafico la soluzione esatta u e
quella approssimata uh .
Stabilire quale relazione intercorre fra ku − uh kL2 (a,b) e h.
Tracciare il grafico della norma dell’errore rispetto ad h utilizzando una scala logaritmica
su entrambi gli assi (loglog).
3

Download Report