Metodi conservativi per
equazioni iperboliche
G. Puppo
Riassunto




Metodo di Lax-Friedrichs
Metodo di Godunov non entropico
Metodo di Godunov con entropy fix
Metodo 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 bc
if nargin < 3
ab=[-1,1];
end
a=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 contorno
if 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')
return
end
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-Friedrichs
Blocco 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
Blocco con le istruzioni di inizializzazione
global bc landa0
if nargin < 5
ab=[-1,1];
end
n=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 superiore
dt = t/nt; landa=dt/h; alpha=abs(cfl); %parametro di diffusione
Notare il modo in cui vengono calcolati dt e alpha: e’ importante
fornire una buona stima di cfl
Calcolo della soluzione (iterazione sul tempo)
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 passo
end
x=linspace(a-h/2,b+h/2,n+2);
Al posto dei ….. devo inserire il blocco delle condizioni al contorno
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
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
Lo script che mi permette di visualizzare la soluzione
all’istante t=0.1 è
global bc landa0
bc='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 on
plot(x,u,'g','Linewidth',2)
axis([-1 1 -1 1.2])
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 landa0
bc='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 on
plot(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 X
n=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;
end
end
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
Listato della funzione che implementa il metodo di Godunov
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
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 passo
end
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 LaxFriedrichs, 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.
Ho usato il seguente script:
global bc landa0
bc='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/2
subplot(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')
continua ...
… 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')
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
L’iterazione sul tempo viene modificata in questo modo.
Inizialmente, calcolo il flusso di Godunov come prima:
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
Aggiungo la correzione entropica
% entropy fix: corregge il flusso se c'e' una rarefazione
% transonica
flux1=f’(u)
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
Uso fzero per calcolare la soluzione di f’(u)=0
Infine calcolo la soluzione e aggiungo le condizioni al contorno,
come negli altri casi.
% 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 passo
end
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 LaxWendroff in uno schema conservativo è il seguente:
Con flusso numerico:
Esempio
Considero 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 LaxFriedrichs
Se però continuo l’integrazione fino allo sviluppo dell’onda d’urto,
il metodo di Lax-Wendroff comincia ad oscillare, ed il metodo di
Lax-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