Equazione del calore

[sophii1]

Si calcoli la soluzione analitica dell'equazione del calore con sorgente

$\{ ((du)/dt (t,x)=(d^2u)/dx^2 (t,x)+2*e ^t*sin(x)) , (u(t,0)=0), ((du)/dx(t,pi/2)=0) , (u(0,x)=sin(x)) :}$ con $t>0$, $0
usando differenze finite del secondo ordine nello spazio e il metodo dei trapezi nel tempo. Si mostrino gli ordini spaziali e temporali della convergenza alla soluzione analitica al tempo tstar=1.

Qui non so proprio da dove iniziare.

a=0;
b=pi/2;
tstar=1;
m=100;
h=(b-a)/(m-1)
s=@(t,x) 2*exp(t).*sin(x);
A = toeplitz(sparse([1,2],[1,1],[-2,1]/h^2,m,1));

[feddy]

Beh, la soluzione analitica si vede che è $u(t,x)=e^{t} sin(x)$, basta guardare le condizioni iniziali...

Usare differenze finite nello spazio significa sostanzialmente discretizzare gli operatori differenziali spaziali, stando attenti a mettere le corrette condizioni al bordo nello spazio.

Ti verrà una cosa del tipo, stando alla bozza del tuo codice:

$\vec{y'} = A \vec{y(t)} + 2 e^{t} \sin(x)$

, dove $A$ è la matrice che discretizza la derivata seconda. Ora è solamente un sistema differenziale da integrare nel tempo usando trapezi

Ora devi mettere le due condizioni ai bordi e integrare nel tempo. Ti faccio notare che il problema è lineare, perciò il tutto si semplifica:

$(I_{N} - \frac{k}{2}A)y_{n+1} = y_n + \frac{k}{2} A y_n + \frac{k}{2} (b(t_{n}) + b(t_{n+1}))$

, dove $b(t)$ è il termine sorgente

Ora devi risolvere n volte questo sistema lineare.

[feddy]

Una possibile implementazione potrebbe essere questa, con l'aggiunta delle condizioni al bordo, che al tempo $\bar{t} =1$ produce la seguente soluzione.
L'impostazione delle condizione ai bordi è facile in questo caso, perché si procede come si fa sempre per schemi alle differenze finite. Cioè, nel caso della condizione di Dirichlet ti è sufficiente azzerare la prima riga. In questo modo si legge $y_{1}^{'} (t) =0$, da cui $y_{1}(t)=\text{cost}$.

Quella di Dirichlet al bordo $x=\frac{\pi}{2}$ si fa ancora aggiungendo un nodo fantasma, e poiché la condizione di Neumann è omogenea, allora $u(t,x_{m+1}) = u(t,x_{m-1})$. Perciò per la derivata seconda si ha che questa è uguale a $\frac{u(t, x_{m+1}) - 2 u(t,x_m) + u(t,x_{m-1})}{h^2}$, da cui segue $\frac{2}{h^2} (u(t,x_{m-1}) - u(t,x_{m)))$

clear all; close all
m=151;
x=linspace(0,pi/2,m)';
h=(pi/2)/(m-1);
A= toeplitz(sparse([1,1],[1,2],[-2,1]/(h^2),1,m));
A(1,1:2)=[0,0];
A(m,m-1:m)=[2/h^2,-2/h^2];
b=@(t) 2*exp(t)*sin(x);

%% ora ho un sistema differenziale da risolvere ! y'=Ay + b;
%f=@(t,y) A*y + b; funzione da azzerare che è lineare
y0 = sin(x);
ts = 500;
k = 1/ts;
y = y0;
I = speye(m);
t = 0;     
for n = 1:ts 
    y(:,n+1) =(I - k*0.5*A)\(y(:,n)+k*0.5*A*y(:,n)+k*0.5*(b(t)+b(t+k))); %trapezi
    t=t+k;
    plot(x,y(:,n),'o')
    title(sprintf('Soluzione a t = %0.2f',t))
    xlabel('x')
    ylabel('u(t,x)')
    axis([0,pi/2,0,3])
    pause(0.001)
end

A te la verifica degli ordini

[sophii1]

Ma il termine sorgente varia da metodo a metodo??

[feddy]

Dipende dall'equazione. Sicuramente non dipende dal metodo numerico che usi.

[sophii1]

Cioè? Riesci a farmi un esempio con un'altra equazione?

[feddy]

Cioè in questo caso il termine sorgente era $s(t,x)=2 e^(t) \sin(x)$. Come vedi non dipende da $u(t,x)$

[feddy]

Un termine sorgente ha quella forma, può essere fatto come ti pare e piace

[sophii1]

nel caso non ci fosse il termine di sorgente,
$y = ((I+k*(1-\theta)*A)*s)/(I-k*\theta*A)$ con s=sorgente
diventa
$y = ((I+k*(1-\theta)*A))/(I-k*\theta*A)$ ??

[sophii1]

perchè se non c'è s significa che s=0, però in questo caso il numeratore sarebbe nullo,e non avrebbe più senso, y sarebbe 0.

[feddy]

A parte che nel messaggio che hai appena scritto stai usando un theta-metodo generico (ma in trapezi si ha $\theta=\frac{1}{2}$), comunque non specifichi $y_n$ e $y_{n+1}$, quindi io leggendo questo

"sophii":$ y = ((I+k*(1-\theta)*A))/(I-k*\theta*A) $

vedo un "rapporto" tra matrici, cosa che non ha senso. Al limite si moltiplica per l'inversa, ma il rapporto non è definito.

Il ragionamento è questo: se togli il termine sorgente, e guardi l'equazione di partenza, allora vedi che il problema discretizzato risulta subito $y'=Ay(t)$. Basta applicare lo stesso ragionamento fatto nella mia prima risposta, usando trapezi: $(I - \frac{k}{2}A) y_{n+1}=(I+ \frac{k}{2} A ) y_{n}$. Come vedi al primo membro si ha un vettore, e al secondo anche.

[feddy]

Per risolvere questo NON SI MOLTIPLICA per

inv(speye(m) - k*0.5*A)

Ma si risolve il sistema lineare corrispondente, dal momento che l'inversione di matrice è mal condizionata in MatLab

[sophii1]

okok capito..grazie!!

[feddy]

prego !