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 condizioni iniziali regolari ed a gradino
Variabili globali
Intervallo di integrazione
Condizioni al contorno Considero 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(x 1) = u(x 2) e u(x. N+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 bc
Fissare il passo di integrazione
Quindi nelle functions che implementano i vari metodi ci sara’ questo blocco di istruzioni: dt=landa 0*h/abs(speed); % arrotonda (T/DT) all'intero immediatamente superiore nt=fix(t/dt)+1; dt = t/nt; landa=dt/h; Inoltre all’inizio della function, dichiaro landa 0 come variabile global bc landa 0
Condizioni iniziali
function init
Costruzione di una function per integrare un’equazione lineare In input devo assegnare: • • u 0 (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 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_lin che applica il metodo di Lax-Friedrichs all’equazione lineare
function lf_lin function [u, x]=lf_lin(u 0, speed, t, ab) % [u, x]=LF_LIN(U 0, SPEED, T) Calcola la soluzione del problema % u_t+speed*u_x=0 su [-1, 1] con il metodo di Lax-Friedrichs. % U 0 e' il vettore che contiene le condizioni iniziali, T e' % l'istante finale. % [u, x]=LF_LIN(U 0, 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 % LANDA 0 (global) e' lo scalare t. c. dt=LANDA 0*h/SPEED Continua. . .
Blocco di inizializzazione global bc landa 0 if nargin < 4 ab=[-1, 1]; end n=length(u 0)-2; a=ab(1); b=ab(2); h=(b-a)/n; dt=landa 0*h/abs(speed); nt=fix(t/dt)+1; % arrotonda (T/DT) all'intero immediatamente superiore dt = t/nt; landa=dt/h; continua. . .
Iterazione sul tempo e output for kt = 1: nt for i=2: n+1 u(i)=0. 5*(u 0(i+1)+u 0(i-1))-speed*landa/2*(u 0(i+1)-u 0(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 u 0=u; % aggiorna u 0 per il prossimo passo end x=linspace(a, b+h, n+2);
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*(u 0(fisico+1)+u 0(fisico-1)) - speed*landa/2*. . . (u 0(fisico+1) – u 0(fisico-1)); . . condizioni al contorno. . u 0=u; % aggiorna u 0 per il prossimo passo end x=linspace(a, b+h, n+2);
Esempio Risolviamo il problema seguente: Con a = 1, e condizioni al contorno di periodicita’
Devo dare questi comandi: clear global bc landa 0 bc='p'; landa 0=0. 9; ux=inline('exp(-5*x. ^2)'); ab=[-2, 2]; u 0=init(ux, 100, ab); [u, x]=lf_lin(u 0, 1, 0. 5, ab); plot(x, u 0) hold plot(x, u, 'g') legend('T = 0', 'T = 0. 5') title('Metodo di Lax-Friedrichs') I comandi sono memorizzati nello 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 u 0(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_lin fisico=2: n+1; for kt = 1: nt if speed >= 0 u(fisico)=u 0(fisico) -landa*speed*(u 0(fisico)-u 0(fisico-1)); else u(fisico)=u 0(fisico) -landa*speed*(u 0(fisico+1)-u 0(fisico)); end ………. u 0=u; % aggiorna u 0 per il prossimo passo end 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-Wendroff Considero 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)=u 0(fisico)-speed*landa/2*(u 0(fisico+1)-u 0(fisico-1)). . . + (speed*landa)^2/2 *(u 0(fisico+1)-2*u 0(fisico) +u 0(fisico-1)); end ……. u 0=u; % aggiorna u 0 per il prossimo passo end 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 0 n=length(x); for i=1: n if x(i) <= 0 y(i) = 2; else y(i) = 0; end
Devo dare questi comandi: clear global bc landa bc='f'; landa='0. 9'; ab=[-1, 4]; u 0=init('step', 100, ab); [u, x]=upwind_lin(u 0, 0. 9, 3, ab); plot(x, u 0) plot(x, u 0, 'Linewidth', 2) hold plot(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
Esercizi Studiare il comportamento dei metodi considerati rispetto a variazioni nei parametri computazionali 1) Con la condizione iniziale u 0(x)=exp(-5 x^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, tipo u 0(x) = sin( x) + sin(10 x) e studiare l'andamento dei metodi al variare di Nx
- Slides: 39