Metodi conservativi per equazioni iperboliche

G. Puppo

Riassunto

Metodo di Lax-FriedrichsMetodo di Godunov non entropicoMetodo di Godunov con entropy fixMetodo di Lax-Wendroff conservativo

Preparazione delle condizioni iniziali

La preparazione delle condizioni iniziali è identica al caso delle equazioni iperboliche lineari. Uso quindi la stessa function init.Questa function richiede:

In input:ux = funzione che definisce l’espressione per u(x,t=0);n = numero di intervalli sui quali si divide l’intervallo [a,b]ab (opzionale) = [a,b] (intervallo di integrazione);

In output:u0 = vettore che contiene i valori iniziali di ux(x) sugli n nodi interni + le condizioni al bordo, date dalla variabile globale bc

Listato della funzione 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/2:h:b+h/2; 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

Metodo di Lax-Friedrichs

Il metodo di Lax-Friedrichs può essere scritto nella forma:

dove è un parametro che dà la diffusione artificiale dello schema: per = 1/, trovo il solito schema di Lax-Friedrichs. Lo schema è stabile per:

Implementazione del metodo di Lax-Friedrichs

In input:u0 = vettore che contiene le medie di cella iniziali;flux = funzione che contiene l’espressione del flusso f(u);t = istante finale;cfl = stima di max|f’(u)| (per la condizione di stabilità);ab (opzionale) = [a,b] intervallo di integrazione;

In output:u = vettore che contiene le medie di cella della soluzione all’istante t x = vettore delle ascisse

Listato del metodo di Lax-FriedrichsBlocco con i commenti:

function [u,x]=lf(u0,flux,t,cfl,ab)% [u,x]=LF(U0,FLUX,T,CFL) Calcola la soluzione del problema% u_t+(flux(u))_x=0 su [-1,1] con il metodo di Lax-Friedrichs.% conservativo% U0 e' il vettore che contiene le condizioni iniziali,T e'% l'istante finale, CFL e' una stima di max|f'(u)|% [u,x]=LF(U0,FLUX,T,CFL,AB) Calcola la soluzione del problema% u_t+(flux(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/CFL

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

Blocco con le istruzioni di inizializzazione

Notare il modo in cui vengono calcolati dt e alpha: e’ importante fornire una buona stima di cfl

for kt = 1:nt % Calcola il flusso: memorizza f(i+1/2) in f(i) fl=feval(flux,u0); for i=1:n+1 f(i)=0.5*(fl(i+1)+fl(i))+alpha/2*(u0(i)-u0(i+1)); end % Calcola la soluzione for i=2:n+1 u(i) = u0(i) - landa*(f(i)-f(i-1)); end ………. u0=u; % aggiorna u0 per il prossimo passoendx=linspace(a-h/2,b+h/2,n+2);

Al posto dei ….. devo inserire il blocco delle condizioni al contorno

Calcolo della soluzione (iterazione sul tempo)

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

Condizioni al contorno

Esempio

Applico il metodo di Lax-Friedrichs al seguente problema:

Attenzione: per l’equazione di Burger, f’(u)=u. Quindi la condizione CFL per questo problema èCFL = max|f’(u)| = max |u(x,t)| = max |u(x,t=0)| = 1.2

global bc landa0bc='p'; landa0=0.9;ux=inline('0.2+sin(pi*x)');flux=inline('0.5*x.^2');u0=init(ux,100,[-1,3]);[u,x]=lf(u0,flux,0.1,1.2);plot(x,u0,'Linewidth',2)hold onplot(x,u,'g','Linewidth',2)axis([-1 1 -1 1.2])

Lo script che mi permette di visualizzare la soluzione all’istante t=0.1 è

Ottengo:

Notare che si sviluppa un’onda di espansione dove u(x,t=0) è crescente, mentre si ha un’onda di compressione dove u(x,t=0) è decrescente

Per visualizzare la soluzione ad istanti successivi, posso modificare lo script precedente in questo modo:

global bc landa0bc='p'; landa0=0.9; nx=100;ux=inline('0.2+sin(pi*x)');flux=inline('0.5*x.^2');u0=init(ux,nx,[-1,3]);[u,x]=lf(u0,flux,0.1,1.2);plot(x,u0,'Linewidth',2)hold onplot(x,u,'g','Linewidth',2)axis([-1 1 -1 1.2])[u,x]=lf(u,flux,0.1,1.2);plot(x,u,'r','Linewidth',2)

La seconda volta, ho chiamato la function lf, usando la soluzione calcolata precedentemente come dato iniziale

Ottengo:

Si stanno sviluppando due onde d’urto

Esercizio

1) Continuare l’integrazione, per tempi più lunghi.Come si evolve la soluzione?Che cosa succede se uso una griglia più fitta?

2) A partire dalla soluzione trovata in t=0.2, integrare all’indietro nel tempo, cioè usare come flusso f(u) = -1/2 u 2 .Il problema è reversibile nel tempo?Che cosa succede se uso una griglia più fitta?Cosa cambia se parto dalla soluzione in t=0.4 o t=0.8?

Esempio

Considero lo stesso problema di prima, ma con un dato iniziale ad onda quadra. Per preparare la condizione iniziale, mi serve ora una function definita in un file .M

function u=quadra(x)% U=QUADRA(X) Crea un dato ad onda quadra su Xn=length(x); a=x(1); b=x(n);x2=(a+b)/2; x1=(a+x2)/2;umin=-1; umax=2;for i=1:n if x(i) < x1 u(i)=umin; elseif x(i) < x2 u(i)=umax; else u(i)=umin; endend

La soluzione ottenuta a diversi istanti è:

Metodo di Godunov non entropico

In assenza di onde di rarefazione transoniche, il metodo di Godunov può essere scritto così:

Questo metodo non è entropico, ma è comunque conservativo

Nel metodo appena scritto, non serve il valore numerico di s’, ma solo il suo segno. Posso quindi riformulare lo schema in questo modo:

In questo modo, lo schema funziona anche se U è costante

function [u,x]=god_ne(u0,flux,t,cfl,ab)% ATTENZIONE: questa e' una versione non entropica del metodo % di Godunov% [u,x]=GOD_NE(U0,FLUX,T,CFL) Calcola la soluzione del % problema% u_t+(flux(u))_x=0 su [-1,1] con il metodo di Godunov.% FLUX e' la funzione di flusso.% U0 e' il vettore che contiene le condizioni iniziali,T e'% l'istante finale, CFL e' una stima di max|f'(u)|% [u,x]=GOD_NED(U0,FLUX,T,CFL,AB) Calcola la soluzione% del problema% u_t+(flux(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/CFL

Listato della funzione che implementa il metodo di Godunov

for kt = 1:nt % Calcola il flusso: memorizza f(i+1/2) in f(i) fl=feval(flux,u0); for i=1:n+1 s = (fl(i+1) - fl(i))*(u0(i+1)-u0(i)); if s >= 0 f(i) = fl(i); else f(i) = fl(i+1); end end % Calcola la soluzione for i=2:n+1 u(i) = u0(i) - landa*(f(i)-f(i-1)); end…………….. u0=u; % aggiorna u0 per il prossimo passoend

x=linspace(a-h/2,b+h/2,n+2);

Al posto dei … devo mettere le condizioni al contorno

Esempio

Applico il metodo di Godunov non entropico e il metodo di Lax-Friedrichs, all’equazione di Burger, con condizione iniziale:

u(x,t=0) = 1 + 1/2 sin ( x) su [-1,3]

con condizioni al contorno periodiche.N.B.: in questo problema non ci sono onde transoniche, perchéf’(u) è sempre diverso da zero.

global bc landa0bc='p'; landa0=0.9; nx=100;flux=inline('0.5*x.^2');ux=inline('1+1/2*sin(pi*x)');u0=init(ux,nx,[-1,3]);T1=0.5; T2=1; cfl=3/2subplot(1,2,1)[u,x]=god_ne(u0,flux,T1,cfl,[-1,3]);plot(x,u,'g','Linewidth',2)axis([-1 3 0.5 1.5])hold on[u,x]=lf(u0,flux,T1,cfl,[-1,3]);plot(x,u,'r','Linewidth',2)title('T=0.5')legend('Godunov','Lax-Friedrichs')

Ho usato il seguente script:

continua ...

subplot(1,2,2)[u,x]=god_ne(u0,flux,T2,cfl,[-1,3]);plot(x,u,'g','Linewidth',2)axis([-1 3 0.5 1.5])hold on[u,x]=lf(u0,flux,T2,cfl,[-1,3]);plot(x,u,'r','Linewidth',2)title('T=1.0')legend('Godunov','Lax-Friedrichs')

… continua

Notare che il metodo di Godunov è meno diffusivo del metodo di Lax-Friedrichs

Ottengo risultati analoghi, se considero un’onda quadra, con i dati umin = 0.2 e umax = 2, in modo da non avere un’onda di rarefazione transonica

Attenzione: qui CFL=2

Se invece scelgo un’onda quadra con umin=-1 e umax=2:

Per questo problema, il metodo di Godunov crea un’onda d’urto non entropica

Metodo di Godunov entropico

Per evitare che il metodo di Godunov dia una soluzione non entropica è necessario aggiungere al flusso numerico di Godunov la correzione entropica

for kt = 1:nt % Calcola il flusso: memorizza f(i+1/2) in f(i) fl=feval(flux,u0); for i=1:n+1 s = (fl(i+1) - fl(i))*(u0(i+1)-u0(i)); if s >= 0 f(i) = fl(i); else f(i) = fl(i+1); end

L’iterazione sul tempo viene modificata in questo modo.Inizialmente, calcolo il flusso di Godunov come prima:

% entropy fix: corregge il flusso se c'e' una rarefazione % transonica fl1=feval(flux1,u0); if fl1(i) <0 & fl1(i+1)> 0 % trova il valore di u per il quale f'(u)=0 if u0(i) >= u0(i+1) us=fzero(flux1,[u0(i+1),u0(i)]); else us=fzero(flux1,[u0(i),u0(i+1)]); end f(i) = feval(flux,us); end end

Aggiungo la correzione entropica

Uso fzero per calcolare la soluzione di f’(u)=0

flux1=f’(u)

% Calcola la soluzione for i=2:n+1 u(i) = u0(i) - landa*(f(i)-f(i-1)); end% Calcola le condizioni al contorno 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 passoend

Infine calcolo la soluzione e aggiungo le condizioni al contorno,come negli altri casi.

Applico questo programma al problema dell’onda quadra:

Ora l’onda di rarefazione è corretta.

Metodo di Lax-Wendroff conservativo

Abbiamo visto che un modo per trasformare il metodo di Lax-Wendroff in uno schema conservativo è il seguente:

Con flusso numerico:

EsempioConsidero l’equazione di Burger, con u(x,t=0)=0.2 + sin(x),su [-1,1]. A t=0.2 ottengo:

Esercizio

Verificare, applicando l’equazione all’indietro nel tempo fino a riottenere la condizione iniziale, che il metodo di Lax-Wendroff è più preciso sia del metodo di Godunov che del metodo di Lax-Friedrichs

Se però continuo l’integrazione fino allo sviluppo dell’onda d’urto,il metodo di Lax-Wendroff comincia ad oscillare, ed il metodo diLax-Friedrichs dà una soluzione migliore

Conclusioni

Abbiamo costruito dei metodi del primo ordine conservativi ed entropici, che però inseriscono una dose elevata di diffusione numerica.Abbiamo costruito uno schema conservativo del secondo ordine, che funziona bene quando la soluzione è regolare, ma che oscilla in prossimità di onde d’urto.

Il prossimo passo è la costruzione di uno schema ibrido, non lineare, che coincida con uno schema del secondo ordine dove la soluzione è regolare, e si riduca ad uno schema del primo ordine vicino alle onde d’urto.

Esercizi

1) Considerare il flusso f(u) = u e verificare che si ottengono gli stessi risultati usando un metodo sviluppato per equazioni lineari e lo stesso metodo in versione conservativa

2) Modificare la funzione del flusso numerico, introducendo il modello del traffico e risolvere il problema del semaforo verde