Metodi numerici per equazioni lineari iperboliche Gabriella Puppo.

Post on 01-May-2015

223 views 2 download

Transcript of Metodi numerici per equazioni lineari iperboliche Gabriella Puppo.

Metodi numerici per equazioni lineari iperboliche

Gabriella Puppo

Riassunto

• Implementazione di uno schema per equazioni lineari iperboliche

• Metodo di Lax-Friedrichs

• Metodo Upwind

• Metodo di Lax-Wendroff

• Esempi con condizioni iniziali regolari ed a gradino

Variabili globali

Ci sono 3 modi per modificare dall’esterno il valore di una variabile all’interno di una function in Matlab:1) impostare la variabile in input;2) modificare la variabile sul listato della function e salvare la modifica;3) definire la variabile come variabile GLOBAL, e modificare la variabile nel workspaceUsero’ questa possibilita’ per impostare dall’esterno le condizioni al contorno nei programmi per le leggi di conservazione.

Intervallo di integrazione

Calcolo la soluzione su un intervallo [A,B], con spaziatura h:Inserisco le condizioni al contorno in x1=A-h/2 e in xN+2=B+h/2.Quindi uso il metodo numerico sui punti di grigliax = A+h, A+2h, … , Bcioe’ per J = 2,…,N+1Nel programma quindi calcolo U(J), per J=2:N+1, mentreh = (B-A) / NAggiungo le condizioni al bordo in J=1 e in J = N+2.In questo modo il vettore soluzione U contiene anche le condizioni al contorno

Condizioni al contornoConsidero due tipi di condizioni al contorno (per ora):

1) Condizione al contorno “free flow” (bc=‘f’): pongo u_x=0 ai bordi del condotto: u(x1) = u(x2) e u(xN+2) = u(x N+1), cioe’ sul programma: U(1) = U(2) e U(N+2) = U(N+1)

2) Condizione al contorno periodica (bc=‘p’)u(A) = u(B) e u(B+h) = u(A+h), cioe’ sul programma U(1) = U(N+1) e U(N+2) = U(2)

Quindi nelle functions che implementano i vari metodi ci sara’ questo blocco di istruzioni

if bc=='f' u(1) = u(2); u(n+2) = u(n+1); elseif bc=='p' u(1) = u(n+1); u(n+2) = u(2); else display('bc non e'' stato impostato') return end

Inoltre all’inizio della function, dichiaro bc come variabile global

global bc

Fissare il passo di integrazione

Il passo di integrazione deve essere fissato in modo da soddisfare la condizione di stabilita’ CFL, che per gli schemi di questo capitolo e’:dt h / |a|, cioè λ 1/ |a|.Per aumentare la flessibilita’ dei nostri programmi, impostero’:dt landa0 h / |a|,Quindi, se landa0 1, soddisfo la condizione CFL. Per studiare il comportamento dei metodi rispetto a λ, landa0 può essere impostata dall’esterno. Ho scelto di dichiarare landa0 come variabile globale, e dunque landa0 deve essere assegnata nel workspace.Inoltre, poiche’ il numero di passi nt deve essere un numero intero, diminuiro’ dt, se occorre, in modo da avere nt intero.

dt=landa0*h/abs(speed);% arrotonda (T/DT) all'intero immediatamente superiorent=fix(t/dt)+1; dt = t/nt; landa=dt/h;

Quindi nelle functions che implementano i vari metodi ci sara’ questo blocco di istruzioni:

Inoltre all’inizio della function, dichiaro landa0 come variabile global

global bc landa0

Condizioni iniziali

Per assegnare le condizioni iniziali, costruisco una function che

1) divide un intervallo [A,B] in N intervalli uguali;

2) calcola l’ampiezza h di questi intervalli, il vettore x dei punti di griglia , e calcola il vettore u0 che contiene il valore di una funzione ux su ognuno dei punti di griglia;

3) aggiunge le condizioni al bordo al vettore u0.

function init

function u0=init(ux,n,ab)% U0=INIT(UX,N) calcola la condizione iniziale per% i metodi per le equazioni iperboliche, applicando% le condizioni al contorno della variabile globale BC% Usa N intervalli su [-1,1], con U0(I)=UX(X(I)).% U0=INIT(UX,N,AB) calcola la condizione iniziale % sull'intervallo AB = [A,B].

global bcif nargin < 3 ab=[-1,1];enda=ab(1); b=ab(2); h=(b-a)/(n);x=a:h:b; u0=feval(ux,x);% Aggiunge le condizioni al contornoif bc=='f' u0(1)=u0(2); u0(n+2)=u0(n+1);elseif bc=='p' u0(1)=u0(n+1); u0(n+2)=u0(2);else display('bc non e'' stato impostato') returnend

Costruzione di una function per integrare un’equazione lineare

• u0 (vettore con le condizioni iniziali),

• la velocita’ di propagazione (speed),

• il tempo finale (T)

• l’intervallo di integrazione ab = [a,b] (se questo parametro manca, il programma assume [a,b] = [-1,1] ).

In input devo assegnare:

In output devo avere:

• Il vettore u, con la soluzione numerica al tempo t

• il vettore x delle ascisse

In questo modo per stampare il grafico bastera’

>> plot(x,u)

All’interno di ogni function:

• Un ciclo sul tempo, nel quale si calcola u all’istante t n;

• all’interno del ciclo sul tempo, ci vuole un ciclo su tutti i punti di griglia, possibilmente costruito vettorialmente

• inoltre all’interno del ciclo sul tempo devo: fissare le condizioni al contorno e copiare la soluzione u appena calcolata sul vettore dei dati all’istante precedente

Metodo di Lax-Friedrichs

Nelle prossime pagine, segue il testo della function lf_linche applica il metodo di Lax-Friedrichs all’equazione lineare

function lf_lin

function [u,x]=lf_lin(u0,speed,t,ab)% [u,x]=LF_LIN(U0,SPEED,T) Calcola la soluzione del problema% u_t+speed*u_x=0 su [-1,1] con il metodo di Lax-Friedrichs.% U0 e' il vettore che contiene le condizioni iniziali,T e'% l'istante finale.% [u,x]=LF_LIN(U0,SPEED,T,AB) Calcola la soluzione del problema% u_t+speed*u_x=0 sull'intervallo AB=[A,B]% Le condizioni al contorno sono contenute nella variabile globale% BC. I valori possibili sono: % BC='p' Condizioni al contorno periodiche% BC='f' Condizioni al contorno free-flow% LANDA0 (global) e' lo scalare t.c. dt=LANDA0*h/SPEED

Continua ...

global bc landa0if nargin < 4 ab=[-1,1];endn=length(u0)-2; a=ab(1); b=ab(2);h=(b-a)/n;dt=landa0*h/abs(speed);nt=fix(t/dt)+1; % arrotonda (T/DT) all'intero immediatamente superioredt = t/nt; landa=dt/h;

Blocco di inizializzazione

continua ...

for kt = 1:nt for i=2:n+1 u(i)=0.5*(u0(i+1)+u0(i-1))-speed*landa/2*(u0(i+1)-u0(i-1)); end if bc=='f' u(1) = u(2); u(n+2) = u(n+1); elseif bc=='p' u(1) = u(n+1); u(n+2) = u(2); else display('bc non e'' stato impostato') return end u0=u; % aggiorna u0 per il prossimo passoendx=linspace(a,b+h,n+2);

Iterazione sul tempo e output

Costruendo il ciclo sui nodi di griglia in modo vettoriale, l'iterazione sul tempo si modifica in questo modo:

fisico=2:n+1;for kt = 1:nt u(fisico)=0.5*(u0(fisico+1)+u0(fisico-1)) - speed*landa/2* ... (u0(fisico+1) – u0(fisico-1));

.... condizioni al contorno ....

u0=u; % aggiorna u0 per il prossimo passoendx=linspace(a,b+h,n+2);

Esempio

Risolviamo il problema seguente:

Con a = 1, e condizioni al contorno di periodicita’

Devo dare questi comandi:

clearglobal bc landa0bc='p'; landa0=0.9;ux=inline('exp(-5*x.^2)');ab=[-2,2];u0=init(ux,100,ab);[u,x]=lf_lin(u0,1,0.5,ab);plot(x,u0)holdplot(x,u,'g')legend('T = 0','T = 0.5')title('Metodo di Lax-Friedrichs')

I comandi sono memorizzati nello script script_lin.m

Ottengo questa figura:

La soluzione iniziale si e’ spostata a destra di 0.5

Se integro per un tempo piu’ lungo, per esempio t=9, ottengo:

Il segnale dovrebbe rimanere immutato, invece il picco si e’ abbassato tantissimo

Esercizio

Studiare il comportamento del metodo di Lax-Friedrichs rispetto ai suoi parametri

1) Con le stesse condizioni iniziali ed al contorno, e mantenendo la stessa griglia dell’esempio precedente, studiare il comportamento della soluzione per =0.3, 0.5, 0.7 0.9 e poi = 0.95, 0.99, 1, 1.01

2) Come nell’esempio precedente, ma con = 0.9, provare a cambiare Nx.

3) Inserire un dato iniziale oscillante, tipo u0(x) = sin( x) + sin( 10 x)

Metodo Upwind

Il metodo upwind e’ definito dalla relazione:

Nella pagina seguente, metto solo la parte centrale della function upwind_lin

function upwind_linfisico=2:n+1;for kt = 1:nt if speed >= 0 u(fisico)=u0(fisico) -landa*speed*(u0(fisico)-u0(fisico-1)); else u(fisico)=u0(fisico) -landa*speed*(u0(fisico+1)-u0(fisico)); end ………. u0=u; % aggiorna u0 per il prossimo passoend

dove, al posto dei puntini, devo mettere le condizioni al contorno

Confrontando i risultati del metodo di Lax-Friedrichs con quelli del metodo Upwind sullo stesso problema di prima, ottengo:

Il metodo upwind da’ una soluzione migliore del metodo di Lax-Friedrichs, tuttavia c’e’ ampio spazio per migliorare...

Esercizio

Risolvere il problema precedente con il metodo di Lax-Friedrichs e con il metodo Upwind, per diversi valori di N. Calcolare l’errore sul picco, sapendo che la soluzione esatta, sul picco, vale 1.

1) Disegnare l’errore in funzione di N per tutti e due gli schemi;

2) A parita’ di N, disegnare l’errore in funzione di landa

Metodo di Lax-WendroffConsidero ora un metodo del secondo ordine, il metodo di Lax-Wendroff:

Anche questa volta, aggiungo soltanto la parte centrale della function lw_lin

function lw_lin

for kt = 1:nt for i=2:n+1 u(fisico)=u0(fisico)-speed*landa/2*(u0(fisico+1)-u0(fisico-1)) ... + (speed*landa)^2/2 *(u0(fisico+1)-2*u0(fisico) +u0(fisico-1)); end ……. u0=u; % aggiorna u0 per il prossimo passoend

dove, al posto dei puntini, devo mettere le condizioni al contorno

Confronto ora i risultati del metodo di Lax-Friedrichs, del metodo Upwind e del metodo di Lax-Wendroff, sempre sullo stesso problema di prima. Ottengo:

Se integro per tempi molto lunghi, ottengo

Esempio

Risolviamo ora il problema

Con condizioni al contorno di free flow e a = 0.9

function step

Questa function crea una funzione con un gradino in x=0

function y=step(x)% Y=STEP(X) Crea una funzione a gradino in 0n=length(x);for i=1:n if x(i) <= 0 y(i) = 2; else y(i) = 0; endend

Devo dare questi comandi:

clearglobal bc landabc='f'; landa='0.9';ab=[-1,4];u0=init('step',100,ab);[u,x]=upwind_lin(u0,0.9,3,ab);plot(x,u0)plot(x,u0,'Linewidth',2)holdplot(x,u,'g','Linewidth',2)axis([-1 4 -0.1 2.1])

Metodo Upwind

Metodo di Lax-Wendroff

Commenti

• I metodi del primo ordine hanno una elevata diffusione artificiale, che diminuisce aumentando N

• I metodi del secondo ordine sono molto più precisi, e il termine principale dell’errore e’ di dispersione

• I metodi del secondo ordine oscillano in prossimita’ dei fronti a gradino: in questo caso, aumentare N non da’ nessun miglioramento

EserciziStudiare il comportamento dei metodi considerati rispetto a variazioni nei parametri computazionali

1) Con la condizione iniziale u0(x)=exp(-5x^2) e condizioni al contorno di periodicita', studiare la diminuzione dell'errore al diminuire di h

2) Implementare il metodo di Beam Warming e studiare il suo comportamento sia su un dato iniziale regolare che su un dato a gradino.

3) Inserire un dato iniziale oscillante, tipou0(x) = sin( x) + sin(10 x) e studiare l'andamento dei metodi al variare di Nx