Metodi FEM per problemi ellittici lineari a tratti

Gabriella Puppo

Principali argomenti

Costruzione della matrice di rigiditàCondizionamentoProblemi misti e di DirichletAndamento dell’ errore

Matrice di rigidità

La matrice di rigidità per il filo elastico, con condizioni al bordo di Dirichlet, è data da:

Moltiplicando per h, ottengo una matrice che dipende solo dal numero degli intervalli

function tridiag

Questa function calcola la matrice tridiagonale del filo elastico:

function a=tridiag(n)% A=TRIDIAG(N) Crea la matrice N per N tridiagonale del% filo elasticoa=2*eye(n);for i=1:n-1 a(i,i+1)=-1; a(i+1,i)=-1;end

La stessa matrice può essere costruita in diversi modi. Il listato a sinistra è meno efficiente, perché non usa le istruzioni vettoriali:

function tridiag2

function a=tridiag2(n)for i=1:n for j=1:n if i==j a(i,j)=2; elseif i==j-1 | i==j+1 a(i,j)=-1; else a(i,j)=0; end endend

function tridiag3

Se invece utilizzo le potenzialità vettoriali di Matlab, posso scrivere la function in questo modo:

function a=tridiag3(n)b=ones(n,1);b1=ones(n-1,1);a=2*diag(b,0)-diag(b1,-1)-diag(b1,1);

Condizionamento della matrice di rigidità

Costruire uno script che calcoli il numero di condizionamento della matrice di rigidità, come funzione del numero N di intervalli, e ne disegni il grafico.

nmax=200;for n=1:nmax a=tridiag(n); c(n)=cond(a);endplot(c,'Linewidth',2)title('Numero di condizionamento')

Lo script precedente fornisce il grafico seguente:

Si dimostra che il numero di condizionamento cresce come N 2

Per verificarlo, possiamo interpolare nel senso dei minimi quadrati con un polinomio di grado 2:

p=polyfit(1:nmax,c,2)p = 4.0529e-01 8.1048e-01 -2.5681e-01

Per verificare che questo polinomio rappresenta accuratamente i nostri dati, sovrapponiamo i grafici del polinomio e del numero di condizionamento

Per verificare che questo polinomio rappresenta accuratamente i nostri dati, sovrapponiamo i grafici del polinomio (in verde) e del numero di condizionamento(pallini blu)

… ma se non sapessi m a priori, come potrei determinarlo?

Esercizio.

Calcolare i polinomi di interpolazione dei dati sul condizionamento per diversi valori del grado m e osservare i coefficienti dei polinomi ottenuti.

Problema del filo elastico, con carico uniforme

Per scrivere una function che risolva il problema del filo elastico con condizioni omogenee di Dirichlet al bordo, e con carico uniforme, dobbiamo:

- Impostare in input il numero di intervalli;

-Costruire la matrice di rigidità

- Preparare il vettore di carico

- Risolvere il sistema lineare

- Aggiungere le condizioni al bordo

function filo_unif

function [u,x]=filo_unif(n)% Risolve il problema del filo elastico % per condizioni al bordo omogenee di% Dirichlet, con carico uniforme% e N punti internia=tridiag(n);h=1/(n+1); % spaziatura di grigliab=-h^2*ones(n,1); % ho portato h a sinistrau=a\b;% Ora aggiungo le condizioni al contornou(2:n+1)=u; u(1)=0; u(n+2)=0;% Calcola il vettore delle ascissefor i=0:n+1; x(i+1) = i*h; end

La soluzione del problema del filo elastico ha questo andamento:

Scriviamo ora una function che ci permetta di calcolare la posizione di equilibrio del filo, per carichi f da definire dall’esterno.

Per risolvere un problema specifico, con una funzione carico assegnata, dovremo creare una function come file.m o come oggetto inline che contenga l’espressione del carico

Problema del filo elastico, con carico assegnato

Per scrivere una function che risolva il problema del filo elastico con condizioni omogenee di Dirichlet al bordo, e con carico assegnato dall’esterno, dobbiamo:

- Impostare in input la funzione carico e il numero di intervalli;

-Costruire la matrice di rigidità

- Preparare il vettore di carico

- Risolvere il sistema lineare

- Aggiungere le condizioni al bordo

function [x,u]=filo_dir(f,n)% Risolve il problema del filo elastico su (0,1)% per condizioni al bordo omogenee di% Dirichlet, con carico definito dalla funzione% F e N punti internia=tridiag(n);h=1/(n+1); % spaziatura di griglia% calcola la posizione dei nodi interni:for i=1:n; x(i) = i*h; endb=h^2*feval(f,x); b=b'; % ho portato h a sinistrau=a\b;% Ora aggiungo le condizioni al contornou(2:n+1)=u; u(1)=0; u(n+2)=0;x(2:n+1)=x; x(1)=0; x(n+2)=1;

Ho ottenuto il seguente listato:

Esempio

Risolviamo il problema del filo elastico soggetto ad un caricodefinito dalla funzione:

f(x) = - sin( 2 x)

Devo dare i seguenti comandi:

>> f=inline('sin(2*pi*x)');>> [x,u]=filo_dir(f,50);>> plot(x,u,'Linewidth',2)

Ottengo:

Vedi lo script script_filo

Risolviamo il problema del filo elastico soggetto ad un caricopuntiforme, cioe’:

f(x) = - delta( x-0.5 )

Esempio 2

function che descrive un carico puntiforme:

function f=delta(x)% F=DELTA(X): Costruisce un vettore riga, che contiene% un carico a forma di delta, quindi F(N/2)=-N, dove% N e' la lunghezza di X, cioe' F(N/2) = -N su un intervallo% di ampiezza 1/Nn=length(x);f=zeros(1,n);if mod(n,2)==0 % controlla se n e' pari f(n/2)= -n/2; f(n/2 +1)= -n/2;else f((n+1)/2)= -n;end

Risolviamo il problema per diversi valori di N:

[x,u]=filo_dir(@delta,9);plot(x,u,'r','Linewidth',2)hold on[x,u]=filo_dir(@delta,99);plot(x,u,'g','Linewidth',2)[x,u]=filo_dir(@delta,999); plot(x,u,'c','Linewidth',2)

Si ottiene il seguente grafico:

Problema misto

Il problema di Poisson, con condizioni al bordo miste, dà matrici che hanno un numero di condizionamento più elevato rispetto al caso con condizioni di Dirichlet.Costruire un programma che costruisca un grafico con:

1) L’andamento del numero di condizionamento per il problema di Dirichlet, in funzione del numero di incognite, N.

2) L’andamento del numero di condizionamento per il problema misto, in funzione del numero di incognite, N.

Esempio (qui la condizione di Neumann è posta in x=0)

% Paragona il numero di condizionamento della matrice del filo% elastico, con condizioni di Dirichlet, e con condizioni miste, a% parita’ di numero di nodi.nmax=100;for n=1:nmax a=tridiag(n); cd(n)=cond(a); a(1,1)=1; cn(n)=cond(a);endni=1:nmax;plot(ni,cd,'go'); hold on; plot(ni,cn,'b+')legend('Dirichlet','Problema misto')title('Numero di condizionamento')

Ottengo questo risultato

Filo vincolato e filo libero ad un estremo

Abbiamo visto come le condizioni al contorno influenzano il numero di condizionamento della matrice di rigidità.

Vorrei ora vedere come le condizioni al contorno influenzano la soluzione, cioè vorrei vedere come la forma del filo dipende dalle condizioni al contorno, a parità di h.

Considero due carichi diversi:1) Carico uniforme: f(x) = - 1;2) Carico puntiforme.

Matrice di rigidità

Per il problema misto, la matrice di rigidità ha una riga in più, einoltre a(1,1) = 1.

function a=tridiag_mista(n)% A=TRIDIAG(N) Crea la matrice (N+1) per (N+1) tridiagonale % del filo elastico, con condizione di Neumann in x=0% e condizione di Dirichlet in x=1a=2*eye(n+1);for i=1:n a(i,i+1)= -1; a(i+1,i)= -1;enda(1,1)=1;

Il programma che impone condizioni miste al bordo è:

function [x,u]=filo_misto(f,n)% Risolve il problema del filo elastico su (0,1)% per condizioni al bordo omogenee di% Neumann in x=0, e di Dirichlet in x=1, % con carico definito dalla funzione% F e N punti interni

continua ...

a=tridiag_mista(n);h=1/(n+1); % spaziatura di griglia% calcola la posizione dei nodi interni:for i=1:n+1; x(i) = (i-1)*h; end %considera anche x=0b=h^2*feval(f,x); b=b'; % ho portato h a sinistrau=a\b;% Ora aggiungo la condizione al contorno% di Dirichlet a destrau(n+2)=0;x(n+2)=1;

…continua

Posso ora paragonare le due soluzioni, con f(x) = 1:

>> f=inline('-1+0*x');>> [x,u]=filo_misto(f,100);>> [xd,ud]=filo_dir(f,100);>> plot(x,u,'Linewidth',2)>> hold on>> plot(xd,ud,'g','Linewidth',2)>> legend('Condiz. miste','Condiz. di Dirichlet')>> title('Filo elastico, f(x) = -1')

Ottengo questa figura:

Per confrontare le due soluzioni nel caso del carico puntiforme,devo dare questi comandi:

[x,u]=filo_misto('delta',100);[xd,ud]=filo_dir('delta',100);plot(x,u,'Linewidth',2)hold onplot(xd,ud,'g','Linewidth',2)legend('Condiz. miste','Condiz. di Dirichlet')title('Filo elastico, f(x) = -delta(0.5)')

Ottengo questa figura

Esercizio 1

Costruire una function che accetti anche condizioni al bordo di Dirichlet non omogenee.Suggerimento: è sufficiente modificare la prima e l’ ultima componente del vettore dei termini noti.

Esercizio 2

Costruire una function che accetti anche una condizione al bordo di Neumann non omogenea.Suggerimento: è sufficiente modificare la prima e l’ ultima componente del vettore dei termini noti, per il problema di Neumann.

Esercizio 3

Calcolare un problema in modo tale da conoscere la soluzione esatta: per esempio, assegno:

Calcolo ora il vettore di carico che darebbe quella soluzione, cioè:

Infine chiamo la function filo_dir, con il carico f(x).Calcolare la velocità di convergenza, al variare di N.