Dubbi su condizione alla Neumann

[john_titor20]

Salve a tutti, sono di nuovo qui a chiedere il vostro immenso aiuto.
Devo risolvere l'equazione di convezione-diffusione \(\displaystyle \frac{\partial{\varphi}}{\partial{t}}+a\frac{\partial{\varphi}}{\partial{x}} -k\frac{\partial^2{\varphi}}{\partial{x^2}}=0\) con:

dominio \(\displaystyle [0, L] \)
condizione iniziale \(\displaystyle \varphi(x, 0)=0 \)
e condizioni al contorno \(\displaystyle \frac{\partial{\varphi}}{\partial{x}}=sin(\frac{2\pi at}{L}) \) a \(\displaystyle x=0 \) e \(\displaystyle \varphi(L, t)=0 \)

su mesh uniforme di ampiezza \(\displaystyle h \) e secondo lo schema \(\displaystyle \varphi^{n+1}_i=\varphi^{n}_i-a\Delta t \frac{\delta_0\varphi^n_i}{2h}+\frac{k\Delta t}{2}( \frac{\delta^2\varphi^n_i+\delta^2\varphi^{n+1}_i}{h^2} )\) considerando che \(\displaystyle \delta_0\varphi_i=\varphi_{i+1}-\varphi_{i-1}\) e \(\displaystyle \delta^2\varphi_i=\varphi_{i+1}-2\varphi_i+\varphi_{i-1}\)

Il mio dubbio ora è applicando il metodo del nodo fantasma dovrei ottenere che \(\displaystyle \frac {\varphi_2-\varphi_0}{2h}=sin(\frac{2\pi at}{L}) \) da cui ricavo che \(\displaystyle \varphi_0=\varphi_2-2h sin(\frac{2\pi at}{L}) \).
Ora la mia domanda è:
ma andando a sostituire nello schema assegnatomi come mi dovrei comportare con la \(\displaystyle \varphi_0^{n+1} \) e la \(\displaystyle \varphi_0^{n} \)? C'è qualche differenza da mettere in evidenza?

[feddy]

Credo ci sia un refuso nel tuo testo: le condizioni al contorno dovrebbero essere $0$ in $x=0$ e $\frac{\partial\varphi}{\partial x}=sin(\frac{2\pi at}{L})$ in $x=L$.

Cosa succede se sostituisci nello schema? Scriviti esplicitamente la prima riga (per un tempo fissato)

[john_titor20]

Innanzitutto grazie per l'aiuto.

Credo ci sia un refuso nel tuo testo

Sì, ho controllato e le condizioni al contorno sono \(\displaystyle \frac{\partial\varphi}{\partial x}=sin(\frac{2\pi at}{L}) \) a \(\displaystyle x=0 \) e \(\displaystyle \varphi(L, t)=0 \) (ho aggiustato per correttezza anche nel primo messaggio)

Cosa succede se sostituisci nello schema? Scriviti esplicitamente la prima riga (per un tempo fissato)

Posto, ora \(\displaystyle i=1 \), ho:
\(\displaystyle \varphi^{n+1}_1=\varphi^{n}_1-a\Delta t (\frac{\varphi^n_2-\varphi^n_0}{2h})+\frac{k\Delta t}{2}( \frac{\varphi^n_2-2\varphi^n_1+\varphi^n_0+\varphi^{n+1}_2-2\varphi^{n+1}_1+\varphi^{n+1}_0}{h^2}) \)

quello che non capisco è come dovrei comportarmi ora, perchè anche fissando \(\displaystyle n=1 \) avrei giustamente delle \(\displaystyle \varphi \) al tempo \(\displaystyle 1 \) e altre al tempo \(\displaystyle 2 \) e perciò come dovrei comportarmi andando a sostituire la \(\displaystyle \varphi_0 \) trovata prima?

[feddy]

Per questa riga che hai scritto, chi è l'incognita? Poni $n=1$, che corrisponde al primo istante temporale, e nota che questo vuol dire $\varphi_i^1=0$ per ogni $i$ (per la condizione iniziale a $t=0$).

[feddy]

Morale della favola: è una equazione in $\varphi_1^{n+1}$

[john_titor20]

Quindi se ho ben compreso, quando ho una situazione di questo tipo devo porre oltre a \(\displaystyle i=1 \) anche \(\displaystyle n=1 \) e di conseguenza considerando lo schema assegnato ho:

\(\displaystyle \varphi^{n+1}_1=\varphi^{n}_1-a\Delta t (\frac{\varphi^n_2-\varphi^n_0}{2h})+\frac{k\Delta t}{2}( \frac{\varphi^n_2-2\varphi^n_1+\varphi^n_0+\varphi^{n+1}_2-2\varphi^{n+1}_1+\varphi^{n+1}_0}{h^2})\)

che diventa

\(\displaystyle \varphi^{2}_1=\varphi^{1}_1-a\Delta t (\frac{\varphi^1_2-\varphi^1_0}{2h})+\frac{k\Delta t}{2}( \frac{\varphi^1_2-2\varphi^1_1+\varphi^1_0+\varphi^{2}_2-2\varphi^{2}_1+\varphi^{2}_0}{h^2})\)

dove considerando che dalla condizione iniziale ricavo che \(\displaystyle \varphi_i^1=0 \) ottengo che:

\(\displaystyle \varphi^{2}_1=\frac{k\Delta t}{2}( \frac{\varphi^{2}_2-2\varphi^{2}_1+\varphi^{2}_0}{h^2})\)

e sostituendo, quindi ora, che \(\displaystyle \varphi_0=\varphi_2-2h sin(\frac{2\pi at}{L}) \), ottengo:

\(\displaystyle \varphi^{2}_1 (1+\beta)-\beta\varphi_2^2=\beta hsin(\frac{2\pi at}{L})\) dove \(\displaystyle \beta=\frac{k\Delta t}{h^2} \)

e così ho i coefficienti da andare a sostituire nella matrice su MatLab

[feddy]

Non ho controllato i conti, ma l'idea è quella: per ogni istante temporale, hai un sistema lineare con una matrice alla quale dovrai cambiare la prima e l'ultima riga per tenere conto delle condizioni ai bordi e un termine noto al quale in prima e ultima entrata dovrai aggiungere eventualmente dei termini

[feddy]

Di fatto, il nodo fantasma lo usi quando ti è richiesto dallo schema il valore puntuale della soluzione in un punto che non appartiene al dominio. In questo tipo di schemi alle differenze finite, capita al bordo ovviamente.

[john_titor20]

Scusa se continuo a tediarti, però credevo di aver capito ma a quanto pare no.
Ho modificato il codice in base a quanto detto qui, però ottengo che le matrici \(\displaystyle A \) e \(\displaystyle B \) sono singolari e di conseguenza la matrice \(\displaystyle T \) è una matrice di \(\displaystyle NaN \)
Metto sotto il codice:

Testo nascosto, fai click qui per vederlo

Codice:

%% Pulizia
clc;
clear;
close all;

%% Variabili
L=2*pi;
k=1;
a=10;
Tf=20;

beta=0.48;
N=20;

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

dt=beta*dx^2/k;
nt=round(Tf/dt);
dt=Tf/nt;

beta=k*dt/(dx^2);
disp(['beta modificato ', num2str(beta)]);
c=a*dt/dx;

t=linspace (0, Tf, nt);
x=x(:);
p=zeros(N, 1);
phi0=0;

%% Plot condizione iniziale
phi=phi0;
figure(1);
plot (x, phi0, '-rx');
title (['it =0 (out of ', num2str(nt), ')']);

%% Operatori differenziali
I=eye(N);
D1=gallery('tridiag', N, 1, 0, 1);
D2=gallery('tridiag', N, 1, -2, 1);
A=I-(beta/2)*D2;
B=I-(c/2)*D1+(beta/2)*D2;

%% Condizioni al contorno
A(N, 1:N)=0; %Dirichlet
B(N, 1:N)=0; %Dirichlet

T=A\B;

T(1, 1)=1+beta; %Neumann
T(1, 2)=-beta; %Neumann

for it=1:nt
    figure(2);
    p(1)=((k*dt)/dx)*sin((2*pi*a*t(it))/L);
    p(N)=0;
    P=A\p;
    phi=T*phi+P;
    plot (x, phi0, '-xr', x, phi, '-ob');
    grid on;
    title(['it= ', num2str(it), ' of ', num2str(nt)]);
    drawnow;
end

[john_titor20]

Ho modificato il codice, correggendo le condizioni al contorno, in particolare le Dirichlet. Ora sembra funzionare, ma vorrei una conferma da voi, in quanto non sono sicuro sulle condizioni alla Neumann. Ve lo chiedo per favore

Testo nascosto, fai click qui per vederlo

Codice:

%% Pulizia
clc;
clear;
close all;

%% Variabili
L=2*pi;
k=1;
a=8;
Tf=20;

beta=0.48;
N=20;

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

dt=beta*dx^2/k;
nt=round(Tf/dt);
dt=Tf/nt;

beta=k*dt/(dx^2);
disp(['beta modificato ', num2str(beta)]);
c=a*dt/dx;
disp(['il numero di Courant è ', num2str(c)]);

t=linspace (0, Tf, nt);
x=x(:);
p=zeros(N-1, 1);
phi0=0*x;
phi=phi0(1:end-1);

%% Plot condizione iniziale
figure(1);
plot (x, phi0, '-rx');
title (['it =0 (out of ', num2str(nt), ')']);

%% Operatori differenziali
I=eye(N-1);
D1=gallery('tridiag', N-1, 1, 0, 1);
D2=gallery('tridiag', N-1, 1, -2, 1);
A=I-(beta/2)*D2;
B=I-(c/2)*D1+(beta/2)*D2;

%% Condizioni al contorno

A(1, 1)=1+beta; %Neumann
A(1, 2)=-beta; %Neumann

T=A\B;

for it=1:nt
    figure(2);
    p(1)=-((k*dt)/dx)*sin((2*pi*a*t(it))/L);
    P=A\p;
    phisol=T*phi+P;
    phi=phisol;
   
    plot (x, phi0, '-xr', x, [phi; 0], '-ob');
    grid on;
    axis([ 0 L -0.2 0.2]);
    title(['it= ', num2str(it), ' of ', num2str(nt)]);
    drawnow;
end


Posto il plot: