Instabilità dello schema discreto; condizione iniziale.

[astrolabio95]

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)^{\prime} + 2/3\varphi_i^{\prime} + 1/6\varphi_(i+1)^{\prime} = (\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

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.

[feddy]

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:

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?

[astrolabio95]

"feddy":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:

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^{\prime} ),( \varphi^{\prime}),( \varphi^{\prime}),( \varphi^{\prime}),( \varphi^{\prime}) ] =[ ( 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

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

li ho "riempiti"

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

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

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?

[feddy]

Dove vedi l'operatore di derivata seconda?

[astrolabio95]

"feddy":Dove vedi l'operatore di derivata seconda?

Ho sbagliato a scrivere!
Intendevo di derivata prima!

[astrolabio95]

Volevo ringraziare qui tutti per il prezioso aiuto.
Oggi ho sostenuto l'esame col massimo dei voti.
Grazie ancora!

[feddy]

Complimenti a te! Contento di esserti stato d'aiuto anche se in piccolissima parte!