Instabilità dello schema discreto; condizione iniziale.

Messaggioda astrolabio95 » 05/07/2019, 17:03

Salve a tutti,

Sto risolvendo questa equazione di convezione lineare non stazionaria

$ (d\varphi)/dt = a(d\varphi)/dx $

nell'intervallo $ [0,L] $ con condizioni al contorno periodiche e condizione iniziale

$ \varphi(x,0) = (x-L)^2cos((2\pix)/L) $.

La traccia richiede di discretizzare la derivata spaziale in questo modo

$ 1/6\varphi_(i-1)^' + 2/3\varphi_i^' + 1/6\varphi_(i+1)^' = (\varphi_(i+1)-\varphi_(i-1))/(2h) $

e di risolvere il sistema di ODE

$ (d\varphi(t))/dt = aD\varphi $

con un metodo Runge-Kutta su quattro passi il cui array è

$ c = (0,1/2,1/2,1) $

$ A_(ij) = | ( 0 , 0, 0, 0),( 1/2, 0, 0, 0),( 0, 1/2, 0, 0),( 0, 0, 1, 0) | $

$ b = (1/6,1/3,1/3,1/6) $

quindi si ha che

$ \varphi_(n+1) = \varphi_n + dt(1/6K1 + 1/3K2 + 1/3K3 + 1/6K4) $ .

Non ci sono stati problemi nel implementare il metodo RK4. Il problema è che, con quella condizione iniziale, la discretizzazione è instabile.
Ho provato anche a semidiscretizzare la derivata spaziale con una derivata centrata su tre punti, ma ho lo stesso risultato.
L'unica cosa che ho notato è che, se la I.C. è un sinusoide di argomento $ (2\pix) $ allora il tutto funziona e non vi sono problemi.

E' un mio errore nello scrivere la I.C. oppure è lo schema che ha questa "particolarità"?

Ecco il codice

Testo nascosto, fai click qui per vederlo
Codice:
function K = Kappa1 (phi) %Funzione che traduce il sistema di ODE

global W

phi = phi(:);

K = W*phi;

K = K(:)';
end

function PHI = RK4(t,x,Phi0) %Metodo RK4

dt = t(2) - t(1);
Nt = length(t);
N = length(x);

PHI = nan(Nt,N);

Phi0 = Phi0(:)';

PHI(1,:) = Phi0;

phi = Phi0;

for it=2:Nt
   
    phi1 = phi;
    K1   = Kappa1(phi1);
   
    phi2 = phi + (dt/2)*K1;
    K2   = Kappa1(phi2);
   
    phi3 = phi + (dt/2)*K2;
    K3   = Kappa1(phi3);
   
    phi4 = phi + dt*K3;
    K4   = Kappa1(phi4);
   
    phi  = phi + dt*(1/6*K1 + 1/3*K2 + 1/3*K3 + 1/6*K4);
   
    PHI(it,:) = phi;
end
end

clc; clear all; close all;

global W

L = 1;
N = 50;
x = linspace (0, L, N+1); x = x(1:end-1); x = x(:);
h = x(2) - x(1);

a = 1;
T = L/a;
Nt = 600;
t = linspace (0, T, Nt);

Phi0 = ((x-L).^2).*cos(2*pi*x/L);

discretizzazione = 1;

switch discretizzazione
   
    case 1
       
    v1 = zeros (1,N);
    v2 = zeros (1,N);

    v1(1)   = 2/3;
    v1(2)   = 1/6;
    v1(end) = 1/6;

    A = gallery ('circul', v1);

    v2(2)   =  1/2;
    v2(end) = -1/2;

    B = gallery ('circul', (1/h)*v2);

    D = (A\B);

    W = sparse(N,N);

    for i=1:N
       W(i,:) = a*D(i,:);
    end
   
    case 2
       
        v3 = zeros (1,N);
        v3(2)    =  1;
        v3(end)  = -1;
       
        W = gallery ('circul', (a/h)*v3);
       
end


PHI = RK4(t, x, Phi0);

PHI = [PHI PHI(:,1)];

for i=1:Nt
    plot ([x;x(end)+h], PHI(i,:), 'b'); grid on;
    drawnow;
end


Grazie.
astrolabio95
Junior Member
Junior Member
 
Messaggio: 125 di 130
Iscritto il: 27/11/2015, 20:39

Re: Instabilità dello schema discreto; condizione iniziale.

Messaggioda feddy » 05/07/2019, 19:56

Ho letto in fretta... come hai imposto le condizioni al bordo periodiche?

Ad ogni modo, questo è il metodo delle linee... discretizzi gli operatori spaziali (in questo caso con un'accuratzza al secondo ordine) e poi risolvi un sistema di ODEs. Nota che visto che la PDE di partenza è lineare, anche la ODE sarà lineare. E' un semplice problema del trasporto (scalare). La matrice che discretizza la derivata prima tra l'altro si può costruire in un colpo solo come ho già fatto in altre discussioni, così ad esempio:

Codice:
D1 = toeplitz(sparse(1,2,-1/(2*h),1,m),sparse(1,2,1/(2*h),1,m));


dove $h$ è l'ampiezza dell'intervallo, m il numero di nodi.

Non ho tempo di guardare la tua implementazione di RK4... la suppongo corretta comunque :)

Non è che hai preso un passo troppo largo?
Avatar utente
feddy
Advanced Member
Advanced Member
 
Messaggio: 2536 di 2542
Iscritto il: 26/06/2016, 01:25
Località: Italia

Re: Instabilità dello schema discreto; condizione iniziale.

Messaggioda astrolabio95 » 05/07/2019, 23:13

feddy ha scritto:Ho letto in fretta... come hai imposto le condizioni al bordo periodiche?

Ad ogni modo, questo è il metodo delle linee... discretizzi gli operatori spaziali (in questo caso con un'accuratzza al secondo ordine) e poi risolvi un sistema di ODEs. Nota che visto che la PDE di partenza è lineare, anche la ODE sarà lineare. E' un semplice problema del trasporto (scalare). La matrice che discretizza la derivata prima tra l'altro si può costruire in un colpo solo come ho già fatto in altre discussioni, così ad esempio:

Codice:
D1 = toeplitz(sparse(1,2,-1/(2*h),1,m),sparse(1,2,1/(2*h),1,m));


dove $h$ è l'ampiezza dell'intervallo, m il numero di nodi.

Non ho tempo di guardare la tua implementazione di RK4... la suppongo corretta comunque :)

Non è che hai preso un passo troppo largo?


Ciao, ti ringrazio per la risposta!

Le BC periodiche le ho implementate con una matrice circolante.
In pratica, ho creato due matrici.
Supponendo di avere un mesh spaziale di cinque punti ho allora

$ [ ( 2/3 , 1/6, 0, 0 , 1/6),( 1/6 , 2/3, 1/6, 0 , 0),( 0, 1/6, 2/3, 1/6, 0),( 0, 0, 1/6, 2/3, 1/6),( 1/6, 0, 0, 1/6, 2/3) ] * [ ( \varphi^' ),( \varphi^'),( \varphi^'),( \varphi^'),( \varphi^') ] =[ ( 0 , 1/(2h), 0, 0 , -1/(2h)),(- 1/(2h) , 0, 1/(2h), 0 , 0),( 0, -1/(2h), 0, 1/(2h), 0),( 0, 0, -1/(2h), 0, 1/(2h)),( 1/(2h), 0, 0, -1/(2h), 0) ] * [ ( \varphi ),( \varphi),( \varphi),( \varphi),( \varphi) ] $

Quindi ho allocato due vettori

Codice:
v1 = zeros (1,N);
v2 = zeros (1,N);


li ho "riempiti"

Codice:
v1(1)    =   2/3;
v1(2)     =  1/6;
v1(end) =  1/6;

v2(2)    =   1/2;
v2(end) = -1/2;


e poi ho costruito le due matrici circolanti

Codice:
A = gallery ('circul', v1);
B = gallery ('circul', (1/h)*v2);


Quindi, tramite il comando backslash ho creato la matrice che discretizza l'operatore di derivata seconda, moltiplicandola poi per la velocità a

Codice:
D = (A\B);

W = sparse(N,N);

for i=1:N
W(i,:) = a*D(i,:);
end


Per quanto riguarda il passo spaziale, ho calcolato il Courant ed è minore dell'unità... Quindi credo sia corretto. Forse sbaglio le BC?
astrolabio95
Junior Member
Junior Member
 
Messaggio: 126 di 130
Iscritto il: 27/11/2015, 20:39

Re: Instabilità dello schema discreto; condizione iniziale.

Messaggioda feddy » 06/07/2019, 10:56

Dove vedi l'operatore di derivata seconda?
Avatar utente
feddy
Advanced Member
Advanced Member
 
Messaggio: 2537 di 2542
Iscritto il: 26/06/2016, 01:25
Località: Italia

Re: Instabilità dello schema discreto; condizione iniziale.

Messaggioda astrolabio95 » 06/07/2019, 11:15

feddy ha scritto:Dove vedi l'operatore di derivata seconda?

Ho sbagliato a scrivere!
Intendevo di derivata prima!
astrolabio95
Junior Member
Junior Member
 
Messaggio: 127 di 130
Iscritto il: 27/11/2015, 20:39

Re: Instabilità dello schema discreto; condizione iniziale.

Messaggioda astrolabio95 » 11/07/2019, 21:10

Volevo ringraziare qui tutti per il prezioso aiuto.
Oggi ho sostenuto l'esame col massimo dei voti.
Grazie ancora!
astrolabio95
Junior Member
Junior Member
 
Messaggio: 128 di 130
Iscritto il: 27/11/2015, 20:39

Re: Instabilità dello schema discreto; condizione iniziale.

Messaggioda feddy » 11/07/2019, 21:28

Complimenti a te! Contento di esserti stato d'aiuto anche se in piccolissima parte!
Avatar utente
feddy
Advanced Member
Advanced Member
 
Messaggio: 2538 di 2542
Iscritto il: 26/06/2016, 01:25
Località: Italia


Torna a Analisi Numerica e Ricerca Operativa

Chi c’è in linea

Visitano il forum: Nessuno e 4 ospiti