1 / 41

Metodi conservativi per equazioni iperboliche

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.

chalsie
Download Presentation

Metodi conservativi per equazioni iperboliche

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Metodi conservativi per equazioni iperboliche G. Puppo

  2. Riassunto • Metodo di Lax-Friedrichs • Metodo di Godunov non entropico • Metodo di Godunov con entropy fix • Metodo di Lax-Wendroff conservativo

  3. 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

  4. 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].

  5. 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

  6. 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:

  7. 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

  8. 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

  9. 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

  10. 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

  11. 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

  12. 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

  13. 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])

  14. 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

  15. 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

  16. Ottengo: Si stanno sviluppando due onde d’urto

  17. 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?

  18. 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

  19. 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

  20. La soluzione ottenuta a diversi istanti è:

  21. 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

  22. 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

  23. 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

  24. 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

  25. 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.

  26. 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 ...

  27. … 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')

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

  29. 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

  30. 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

  31. 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

  32. 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

  33. Aggiungo la correzione entropica % 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 flux1=f’(u) Uso fzero per calcolare la soluzione di f’(u)=0

  34. 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

  35. Applico questo programma al problema dell’onda quadra: Ora l’onda di rarefazione è corretta.

  36. 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:

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

  38. 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

  39. 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

  40. 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.

  41. 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

More Related