Implementazione delle condizioni al contorno per un'equazione non stazionaria

[astrolabio95]

Salve a tutti,
Sono uno studente di ingegneria aerospaziale alle prese con l'esame di metodi numerici.
In questo esame si tratta la risoluzione di PDE con il metodo delle differenze finite, sia stazionarie che non stazionarie.

Purtroppo la materia, che trovo interessante, è piuttosto complessa per me in quanto manca l'adeguato materiale didattico; esistono solo dei codici compilati da qualche collega, senza un po' di commento.

La teoria l'ho studiata qua e la da qualche libro, e pertanto trovo difficoltà, non tanto nel campo dello stazionario, quanto in quello della non stazionarietà ed in particolare nell'implementare le opportune condizioni al contorno.
Anche perché molti testi sono molto specifici, mentre il mio esame consta di soli 6 CFU.

Fatto questo preambolo, vengo al problema.
Ho, ad esempio, questa equazione di diffusione non stazionaria

$ (d\varphi)/dt=k(d^2\varphi)/dx^2 $

Il testo chiede di risolverla con un metodo FTCS, ovvero Eulero esplicito.
A questo punto so che la derivata temporale la si svolge tra due intervalli di tempo mentre quella spaziale viene effettuata per i soli punti interni del mesh.
In particolare, se suppongo un mesh eqispaziato di passo h avrò la seguente discretizzazione

$ \varphi_i^(n+1)=\beta\varphi_(i-1)^n+(1-2\beta)\varphi_i^n+\beta\varphi_(i+1)^n $

dove $ \beta=(k\Deltat)/(h^2) $.

Cioè in questo caso, essendo il mesh spaziale uniforme, la matrice che discretizza la derivata seconda nei punti interni del mesh è tridiagonale ed ha sulle diagonali 1, -2, 1.

A questo punto, però, se ho una condizione alla Neumann non so come procedere.
Il testo mi consiglia una discretizzazione della derivata prima sui due nodi estremanti asimmetrica e del primo ordine; inoltre mi suggerisce di applicarla solo ai nodi all'intervallo di tempo n+1.

Qualcuno sa darmi qualche suggerimento sul da farsi?

Grazie a tutti.

[feddy]

Non ho molto tempo, comunque si chiama metodo delle linee: prima discretizzi gli operatori differenziali spaziali con differenze finite (in questo caso centrate del secondo ordine), così da ridurti ad una ODE, poi integri nel tempo, stando attento al metodo numerico utilizzato, in quanto il problema è in generale "stiff", a causa della matrice che hai menzionato sopra.


Per semplicità, chiamo $\varphi=u$, e supponi di avere una condizione tale PDE con condizione al bordo $x=b$ data da $\frac{\partial u}{\partial x}(t,b)=0$.

Introduci la variabile $y_{m+1}(t) \approx u(t,x_{m+1})$, dove $x_{m+1}=b+h$, e imponi $y_{m+1}(t)=y_{m-1}(t)$.

Quindi, quando vai ad approssimare la derivata seconda nel tuo problema, avrai al bordo:

$\frac{\partial^2 u}{\partial x^2}(t,b) \approx \frac{u(t,x_{m+1})-2*u(t,x_m)+u(t,x_{m+1} )} {h^2}=\frac{y_{m+1}(t)-2y_m(t)+y_{m-1}(t)}{h^2}=\frac{2*y_{m-1}(t)-2y_m(t)}{h^2}$

Quindi nella matrice tridiagonale dovrai cambiare le entrate $(m,m-1)$ e $(m,m)$ con rispettivamente $2$ e $-2$, ricordandoti di dividere per $h^2$!

[astrolabio95]

Ciao, ti ringrazio per la celere risposta.

Il problema che mi ritrovo ad affrontare è il modo con cui implementare questa condizione in MatLab.
Cioè, per i nodi interni non c'è problema... il problema mi si pone quando vado a discretizzare per i nodi di bordo.

Ad esempio, se ho N nodi che discretizzano il mesh, devo allocare una matrice NXN?

[feddy]

Rileggi l'ultima frase nella mia risposta... è quello che c'è da fare nella pratica per nua condizione di Neumaann omogenea. Insomma, come per quelle di Dirichlet, vanno cambiate la prima e l'ultima riga della matrice che discretizza il problema.

Ovvio che devi allocare una matrice NxN, e tale matrice sarà chiaramente sparsa, con, a priori, $N+2*(N-1)=3N-2$ elementi diversi da 0

[astrolabio95]

Ti ringrazio per la risposta, ho risolto il quesito!

Nel caso avessi una derivata prima (quindi dovessi discretizzare Neumann) da discretizzare non su due punti ma, ad esempio, su quattro?

Ad esempio, se mi metto in x = b, posso effettuare una discretizzazione su quattro punti asimmetrica?
Il procedimento è lo stesso, solo che devo cambiare gli ultimi quattro posti della matrice, giusto?

Grazie!

[feddy]

Felice di esserti stato utile!
Sostanzialmente vuoi usare differenze finite non centrate. Sì, sostanzialmente userai i punti a sinistra. Di solito ci sono già le formule, ad esempio, per discretizzare la derivata prima al bordo di destra usando due punti a sinistra, o cose così.

[astrolabio95]

Ti ringrazio ancora per l'aiuto ed il prezioso supporto...

Il fatto è che all'esame il professore di diverte coi più disparati metodi di interpolazione; stencil asimmetrici, con punti iniziali particolari e condizioni al contorno periodiche la cui implementazione faccio fatica a comprendere...

Non riesco a trovare nemmeno un testo che tratti, in maniera chiara, le equazioni del calore (diffusione e convezione).

[astrolabio95]

Ciao, perdona ancora il disturbo..

Ho provato a svolgere questa traccia, ma i risultati sono deludenti...





Questo è il mio codice...

Testo nascosto, fai click qui per vederlo

Codice:

clc; clear all; close all;

N = 30;
L = 1;
x = linspace(0, L, N);
h = x(2) - x(1);

k = 1;
beta = .4;
Dt = (h*h*beta)/k;
Tf = .5;
Nt = round(Tf/Dt);


Dtilde = gallery('tridiag', N, 1, 2, 1);
I = eye(N);

A = I + beta*Dtilde;

full(A)


full(A)

Itilde = eye(N);
Itilde([1,end], :) = 0;

full(Itilde)

ptilde = zeros(N, 1);
qtilde = ptilde + Dt*(1-exp(-1))/(1+exp(-1))
Phi0 = zeros(N,1);

Phinm = Phi0;

Phin = (A\Itilde)*Phinm + A\qtilde;

full(Phin)

plot (x, Phin, 'r--'); grid on;

Astar = I + 2*beta*I;

full(Astar)

Astar([1,end],:) = 0;

Astar(1,1) = 3/2*h;
Astar(1,2) = -2/h;
Astar(1,3) = 1/2*h;
Astar(end,end) = 1;

full(Astar)

B = gallery('tridiag', N, 2*beta, 0, 2*beta);

full(B)

B([1,end],:) = 0;

full(B)

C = I - 2*beta*I;

C([1,end],:) = 0;

full(C)


p = zeros(N,1);

for it = 2:Nt
    p = p + 2*Dt*(1-exp(-it))/(1+exp(-it));
   
    Phinp = (Astar\B)*Phin + (Astar\C)*Phinm + (Astar\p);
   
    Phin = Phinp;
   
    Phnm = Phin;
   
    plot (x, Phin, 'b'); grid on;
    drawnow;
end


Non riesco a capire dove sbaglio, ho seguito le indicazioni di qualche collega ma il mio risultato oscilla...fa schifo insomma

Ti ringrazio

[astrolabio95]

Ho fatto una "versione" senza quella strana forzante e con condizione iniziale di tipo sinusoidale...e pare venga.
Ho, inoltre, rivisto alcune cosette...
Se puoi buttarci un occhio te ne sarei infinitamente grato!

Testo nascosto, fai click qui per vederlo

Codice:

clc; clear all; close all;

N = 23;
L = 1;
x = linspace(0, L, N);
h = x(2) - x(1);

k = 1;
beta = .4;
Dt = (h*h*beta)/k;
Tf = .2;
iTf = ceil(Tf/Dt);
Dt = Tf/iTf;

A = sparse(N,N);
D = gallery('tridiag', N, 1, 2, 1);
I = eye(N);

A = I + beta*D;

full(A)

A([1,end], :) = 0;

A(1,1) = 3/2*h;
A(1,2) = -4/2*h;
A(1,3) = 1/2*h;
A(end,end) = 1;

full(A)

Itilde = eye(N);
Itilde([1,end], :) = 0;

full(Itilde)

p0 = zeros(N, 1);
q0 = Dt*p0;

x = x(:);
Phi0 = sin(2*pi*x);
Phinm = Phi0;

Phin = (A\Itilde)*Phinm + (A\q0);

DD =sparse(N,N);
DD = gallery('tridiag', N, 0,2,0);

B = I + beta*DD;

B([1,end],:) = 0;

B(1,1) = 3/2*h;
B(1,2) = -4/2*h;
B(1,3) = 1/2*h;
B(end,end) = 1;

full(B)

C = sparse(N,N);
C = gallery('tridiag', N, 2*beta,0,2*beta);

C([1,end],:) = 0;

full(C)

D = sparse(N,N);
D = I - beta*DD;

D([1,end],:) = 0;

full(D)

p = zeros(N, 1);

q = 2*Dt*p;

dPhi = 1;

while dPhi > 0.0100
   
    Phinp = (B\C)*Phin + (B\D)*Phinm + (A\q);
   
    dPhi = max(abs(Phinp-Phinm))/Dt;
   
    Phin = Phinp;
   
    Phinm = Phin;
   
    plot (x, Phi0, 'r--', x, Phin, '-b');
    axis([0 L -1.5 1.5]); grid on;
    title(['Evoluzione temporale. dPhi = ', num2str(dPhi)]);
    drawnow;
end


[feddy]

Onestamente non ho guardato bene il codice, anche perché ora non ho molto tempo. L'imposizione della condizione di Neumann con FD decentrate va bene: hai cambiato la prima riga della matrice di discretizzazione.