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;endh=(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);endplot(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 internih=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 internih=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 internih=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 internih=1/(m+1);xx=linspace(h,1-h,m);% inizializzazionefx=feval(f,xx); fx=h*fx'; % termine forzanteb=massa_sp(m); % matrice di massaa=stiff_sp(m); % matrice di rigidita'dt=t/nt; % passo di integrazioneba= b-dt*nu*a; % matrice che moltiplica c^n ad ogni passo[ll,uu]=lu(b); % fattorizzazione matrice di massac=u0; % condizione iniziale

Iterazione nel tempo e output

%% iterazionefor n=0:nt-1 noto = ba*c+dt*fx; c=uu\(ll\noto);end% aggiunge le condizioni al contornou(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 integrazionee 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);

Da cui si vede che u è identicamente uguale a zero

>> norm(u)ans = 0

Verifico che ho effettivamente ottenuto u=0.

Per esempio, calcolo la norma di u:

Esempio 2

Consideriamo un esempio meno banale

Mi aspetto che u(x,t) tenda a zero, per t che tende all’infinito.

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);

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:

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

m=length(u0); % calcola il numero di nodi internih=1/(m+1);xx=linspace(h,1-h,m);% inizializzazionefx=feval(f,xx); fx=h*fx'; % termine forzanteb=massa_sp(m); % matrice di massaa=stiff_sp(m); % matrice di rigidita'dt=t/nt; % passo di integrazioneba= b+dt*nu*a; % matrice che moltiplica c^(n+1) ad ogni passo[ll,uu]=lu(ba); % fattorizzazione matrice BAc=u0; % condizione iniziale

Blocco di inizializzazione

%% iterazionefor n=0:nt-1 noto = b*c+dt*fx; c=uu\(ll\noto);end% aggiunge le condizioni al contornou(1)=0; u(2:m+1)=c; u(m+2)=0;x(1)=0; x(2:m+1)=xx; x(m+2)=1;

Iterazione principale

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 è stabilee 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à