Soluzione FEM di problemi
parabolici
G. Puppo
Riassunto
•
•
•
•
•
Condizionamento della matrice di massa
Metodi FEM lineari a tratti
Metodo FEM con Eulero esplicito
Metodo FEM con Eulero implicito
Esempi
Condizionamento della matrice
di massa
La matrice di massa è definita dalla seguente equazione
Scrivo una function che calcoli la matrice di massa su un
intervallo [A,B], con N nodi interni.
function massa.m
function b=massa(n,a,b)
% B=MASSA(N,A,B)
% Calcola la matrice di massa sull'intervallo
% [a,b] con n nodi interni
% B=MASSA(N) assume [A,B]=[0,1]
if nargin==1
a=0; b=1;
end
h=(b-a)/(n+1);
d1=1/6*ones(n-1,1);
b=h*(2/3*eye(n)+diag(d1,1)+diag(d1,-1));
Condizionamento della matrice
di massa
Stimo il numero di condizionamento della matrice di massa
sull’intervallo (0,1) in funzione di N
nmax=200;
for n=1:nmax
a=massa(n);
c(n)=cond(a);
end
plot(c)
title('Numero di condizionamento')
Condizionamento della matrice di massa
Soluzione FEM di problemi
parabolici
Consideriamo il seguente problema modello:
FEM lineari a tratti
Usando funzioni di test e di prova lineari a tratti, il problema
precedente si riduce alla risoluzione del sistema di equazioni
differenziali:
Dove c è il vettore delle funzioni incognite, B è la matrice di
massa, A è la matrice di rigidità e F è il vettore dei termini noti
Metodo di Eulero esplicito
Se integro il sistema di equazioni differenziali con il metodo di
Eulero esplicito, ottengo il seguente sistema lineare di equazioni
algebriche:
Qui t è il passo di integrazione, che deve soddisfare la
condizione di stabilità.
Per avere un metodo efficiente, devo memorizzare A e B come
matrici sparse.
FEM lineari a tratti, con Eulero Esplicito
• Input: f, vettore u0 delle condizioni iniziali,
tempo finale T, diffusione NU, numero
passi NT.
• Inizializzazione: calcolare t, A, B, F c0
• Fattorizzare B una volta per tutte
• Risolvere il sistema algebrico, n=0,1,…
• Stampare la soluzione finale
Condizione iniziale
Costruisco una function che calcoli il vettore delle condizioni iniziali
function c0=init(u,m)
% C0=INIT(U,M) Calcola il vettore C0 delle
% condizioni iniziali,C0(I)=U(X(I)),
% sull'intervallo (0,1) con M nodi interni
h=1/(m+1);
x=linspace(h,1-h,m);
c0=feval(u,x)';
N.B. Calcola U solo sui nodi interni
function spdiags
Abbiamo visto che per ottenere un programma efficiente è
essenziale sfruttare la sparsità delle matrici che compaiono nello
schema
Per costruire matrici a struttura diagonale in forma sparsa, posso
usare la function spdiags
A = SPDIAGS(B,d,m,n) creates an m-by-n sparse matrix from the
columns of B and places them along the diagonals specified by d.
Costruzione matrice di massa
Poiché è necessario calcolare diversi prodotti della forma B*x,
conviene memorizzare la matrice B come matrice sparsa
function b=massa_sp(n)
% B=MASSA_SP(N) Calcola la matrice di massa
% in forma sparsa con N nodi interni
h=1/(n+1);
bvect=[1/6*ones(n,1), 2/3*ones(n,1), 1/6*ones(n,1)];
b=spdiags(bvect,[-1,0,1],n,n);
b=h*b;
Costruzione matrice di rigidità
Anche in questo caso, conviene utilizzare matrici memorizzate
in forma sparsa
function a=stiff_sp(n)
% A=STIFF_SP(N) Calcola la matrice di stiffness
% in forma sparsa con N nodi interni
h=1/(n+1);
bvect=[-ones(n,1), 2*ones(n,1), -ones(n,1)];
a=spdiags(bvect,[-1,0,1],n,n);
a=a/h;
Metodo di Eulero Esplicito
Intestazione e commenti:
function [u,x]=fem1_ee(f,u0,t,nu,nt)
% [U,X]=FEM1_EE(F,U0,T,NT,NU) calcola la soluzione con
% FEM lineari a tratti del problema parabolico
% u_t -nu u_xx = f su (0,1), con condizione
% iniziale u0 (vettore colonna).
% integra nel tempo con Eulero esplicito
Questa function lavora sull’intervallo [0,1]
Blocco di inizializzazione
m=length(u0); % calcola il numero di nodi interni
h=1/(m+1);
xx=linspace(h,1-h,m);
% inizializzazione
fx=feval(f,xx); fx=h*fx'; % termine forzante
b=massa_sp(m); % matrice di massa
a=stiff_sp(m); % matrice di rigidita'
dt=t/nt;
% passo di integrazione
ba= b-dt*nu*a; % matrice che moltiplica c^n ad ogni passo
[ll,uu]=lu(b); % fattorizzazione matrice di massa
c=u0;
% condizione iniziale
Iterazione nel tempo e output
%
% iterazione
for n=0:nt-1
noto = ba*c+dt*fx;
c=uu\(ll\noto);
end
% aggiunge le condizioni al contorno
u(1)=0; u(2:m+1)=c; u(m+2)=0;
x(1)=0; x(2:m+1)=xx; x(m+2)=1;
Esempio 1
Supponiamo di voler risolvere il problema:
In questo caso la soluzione esatta è u(x,t) = 0
Supponiamo di voler risolvere questo problema fino a T=1.
Proviamo ad usare M=9 nodi interni, NT=20 passi di integrazione
e NU=2:
Calcolo la condizione iniziale:
>> ux=inline('0*x');
>> u0=init(ux,9)
Preparo la forzante:
>> f=inline('0*x');
Chiamo la function per Eulero esplicito
[u,x]=fem1_ee(f,u0,1,2,20);
Verifico che ho effettivamente ottenuto u=0.
Per esempio, calcolo la norma di u:
>> norm(u)
ans =
0
Da cui si vede che u è identicamente uguale a zero
Esempio 2
Consideriamo un esempio meno banale
Mi aspetto che u(x,t) tenda a zero, per t che tende all’infinito.
Considero NU=1, T=0.01 (la soluzione sarà ancora diversa da
zero) e 99 nodi interni.
Provo con 20 passi di integrazione.
Devo dare questi comandi:
f=inline('0*x');
ux=inline('exp(x).*sin(pi*x)');
u0=init(ux,99);
[u,x]=fem1_ee(f,u0,0.01,1,20);
Per visualizzare il grafico:
plot(x,u)
Non è esattamente quello che volevo...
Diminuendo drasticamente il passo di integrazione, si trova:
Metodo di Eulero implicito
Se integro il sistema di equazioni differenziali con il metodo di
Eulero implicito, ottengo il seguente sistema lineare di equazioni
algebriche:
Qui t è il passo di integrazione, che non deve più soddisfare la
condizione di stabilità.
Per avere un metodo efficiente, devo memorizzare A e B come
matrici sparse.
Listato del metodo implicito
function [u,x]=fem1_ei(f,u0,t,nu,nt)
% [U,X]=FEM1_EI(F,U0,T,NT,NU) calcola la soluzione con
% FEM lineari a tratti del problema parabolico
% u_t -nu u_xx = f su (0,1), con condizione
% iniziale u0 (vettore colonna).
% integra nel tempo con Eulero implicito
Blocco di inizializzazione
m=length(u0); % calcola il numero di nodi interni
h=1/(m+1);
xx=linspace(h,1-h,m);
% inizializzazione
fx=feval(f,xx); fx=h*fx'; % termine forzante
b=massa_sp(m); % matrice di massa
a=stiff_sp(m); % matrice di rigidita'
dt=t/nt;
% passo di integrazione
ba= b+dt*nu*a; % matrice che moltiplica c^(n+1) ad ogni passo
[ll,uu]=lu(ba); % fattorizzazione matrice BA
c=u0;
% condizione iniziale
Iterazione principale
%
% iterazione
for n=0:nt-1
noto = b*c+dt*fx;
c=uu\(ll\noto);
end
% aggiunge le condizioni al contorno
u(1)=0; u(2:m+1)=c; u(m+2)=0;
x(1)=0; x(2:m+1)=xx; x(m+2)=1;
Confronto Eulero esplicito e
implicito
Verifichiamo che il metodo implicito permette di ottenere la
stessa precisione con un numero di passi minore rispetto al
metodo esplicito.
Uso lo stesso problema di prima, quindi:
>> u0=inline('sin(pi*x).*exp(x)');
>> c0=init(u0,19);
>> f=inline('0*x');
>> [ue,x]=fem1_ee(f,c0,0.1,1,240);
>> [ui,x]=fem1_ei(f,c0,0.1,1,240);
Con questi dati, anche il metodo di Eulero esplicito è stabile
e le due soluzioni sono comparabili:
Con il metodo implicito posso diminuire drasticamente il numero
dei passi di integrazione, per esempio:
>>[ui,x]=fem1_ei(f,c0,0.1,1,40);
>> plot(x,ui,'r+','Linewidth',2)
La soluzione cambia in modo
appena percettibile
Viceversa, se diminuisco anche di poco il numero dei passi di
integrazione con il metodo esplicito, per esempio:
>>[ue,x]=fem1_ee(f,c0,0.1,1,200);
>> plot(x,ue,'r','Linewidth',2)
La soluzione cambia drasticamente
Esercizio 1
Considerare la matrice AB=B\A che compare nel sistema
semidiscreto agli elementi finiti. Il metodo di Eulero esplicito e’
stabile per | 1 + Δt λ| < 1, dove λ indica un autovalore di AB.
1) Calcolare gli autovalori di AB per ν fissato e disegnarli . 2)
Calcolare graficamente per quali Δt tutti gli autovalori finiscono
nella regione di stabilita’ e verificare che il Metodo di Eulero
esplicito e’ stabile solo per i valori di Δt che soddisfano la
condizione trovata.
Esercizio 2
Considerare un dato iniziale molto oscillante, tipo:
U0(x) = somme(su k) sin( πx k),
con f(x)=0.
Stampare la soluzione a diversi istanti e osservare quali
componenti si smorzano prima. Confrontare con la stima teorica,
ottenuta con la serie di Fourier
Esercizio 3
Modificare il programma, in modo da poter inserire
una fonte di calore dipendente dal tempo.
In questo caso, f deve essere funzione di due variabili, x e t, e il
vettore del carico deve essere calcolato all’ interno del ciclo sul
tempo.
Esercizio 4
Inserire un termine di convezione, cioè risolvere il
problema:
In questo caso è necessario modificare la matrice di rigidità