Equazione del calore.

Alvise Sommariva

Università degli Studi di Padova

Dipartimento di Matematica

6 giugno 2012

Equazione del calore.

Consideriamo l’equazione del calore

\{ \begin{matrix}\frac{∂u}{∂t} = \frac{∂^{2}u}{∂x^{2}} + G, 0 < x < 1, t > 0 \\ u(0, t) = d_{0}(t), u(1, t) = d_{1}(t), t ≥ 0 \\ u(x, 0) = f (x), 0 ≤ x ≤ 1 \end{matrix}  (1)

Sia m > 0 intero e sia h_{x} = 1/m ed x_{j} = jh_{x} con j = 0, 1, …, m. Si può mostrare che per j = 1, 2, …, m − 1 e ξ_{j} ∈ (x_{j−1}, x_{j+1})

\frac{∂^{2}u}{∂x^{2}}(x_{j}, t) = \frac{u(x_{j+1}, t) − 2u(x_{j}, t) + u(x_{j−1}, t)}{h_{x}^{2}} − \frac{h_{x}^{2}}{12} \frac{∂^{4}u}{∂x^{4}}(ξ_{j}, t)

(2)

Alvise Sommariva  Equazione del calore.  2/ 24

Equazione del calore.

Tralasciando il termine finale e posto u_{j}(t) := u(x_{j}, t) otteniamo quindi per j = 1, …, m − 1 il sistema di equazioni differenziali

u′_{j}(t) = \frac{u_{j+1}(t) − 2u_{j}(t) + u_{j−1}(t)}{h_{x}^{2}} + G(x_{j}, t)  (3)

Risolto (3), si avrà una approssimazione della soluzione dell’equazione del calore per x_{j} = jh_{x} e t ≥ 0. Il procedimento appena descritto è noto in letteratura come metodo delle linee. Nel risolvere il sistema dobbiamo far attenzione alle condizioni sul bordo

u_{0}(t) = d_{0}(t), u_{m}(t) = d_{1}(t)

e ricordare che la condizione iniziale del sistema di equazioni differenziali è

u_{j}(0) = f (x_{j}), j = 1, …, m − 1.

Il sistema differenziale (3) può essere riscritto matricialmente. Posto

u(t) := [u_{1}(t), …, u_{m−1}(t)]^{T}
u_{0}(t) := [f (x_{1}), …, f (x_{m−1})]^{T}

g(t) := [\frac{1}{h_{x}^{2}}d_{0}(t), …, \frac{1}{h_{x}^{2}}d_{1}(t)]^{T}+[G(x_{1}, t), 0, …, 0, G(x_{m−1}, t)]^{T}

Λ = \frac{1}{h_{x}^{2}} [ \begin{matrix}−2 & 1 & 0 & 0 & … & 0 \\ 1 & −2 & 1 & 0 & … & 0 \\ 0 & 1 & −2 & 1 & … & 0 \\ … & … & … & … & … & … \\ 0 & 0 & 0 & 0 & 1 & −2 \end{matrix} ]  (4)

otteniamo che (3) è equivalente al sistema di equazioni differenziali (lineari)

u′(t) = Λu(t) + g(t), u(0) = u_{0}(t)  (5)

Equazione del calore.

Tra i metodi più comuni nel risolvere il problema differenziale (di Cauchy)

\{ \begin{matrix}u′(t) = F (t, u(t)) \\ u(0) = u_{0} \end{matrix}  (6)

citiamo il metodo di Eulero esplicito, per u_{n+1} ≈ u(t_{n+1})

\{ \begin{matrix}u_{n+1} = u_{n} + hF (t_{n}, u_{n}) \\ u_{0} assegnato \end{matrix}  (7)

e quello di Eulero implicito

\{ \begin{matrix}u_{n+1} = u_{n} + hF (t_{n+1}, u_{n+1}) \\ u_{0} assegnato \end{matrix}  (8)
Alvise Sommariva  Equazione del calore.  5/ 24

Equazione del calore.

Nel nostro caso

F (t, v(t)) := Λv(t) + g(t)

e quindi il metodo di Eulero esplicito genera la successione

\{ \begin{matrix}v_{n+1} = v_{n} + h_{t}(Λv_{n} + g(t_{n})) \\ v_{0} assegnato \end{matrix}  (9)

mentre Eulero implicito determina

\{ \begin{matrix}v_{n+1} = v_{n} + h_{t}(Λv_{n+1} + g(t_{n+1})) \\ v_{0} assegnato \end{matrix}  (10)

o equivalentemente

\{ \begin{matrix}(I − h_{t}Λ)v_{n+1} = v_{n} + h_{t}g(t_{n+1}) \\ v_{0} assegnato \end{matrix}  (11)
Alvise Sommariva  Equazione del calore.  6/ 24

Osserviamo che a differenza del metodo esplicito, poichè

\{ \begin{matrix}(I − h_{t}Λ)v_{n+1} = v_{n} + h_{t}g(t_{n+1}) \\ v_{0} assegnato \end{matrix}  (12)

ad ogni iterazione si richiede la soluzione di un sistema (che nel nostro caso è lineare). Usando i primi due teoremi di Gerschgorin, si può mostrare che la matrice (I − h_{t}Λ) è definita positiva (e quindi non singolare). A partire da Eulero esplicito ed Eulero implicito si definiscono i cosidetti θ−metodi in cui

v_{n+1} = (1 − θ) (v_{n} + h_{t}(Λv_{n} + g(t_{n})))
+θ (v_{n} + h_{t}(Λv_{n+1} + g(t_{n+1})))  (13)

con v_{0} assegnato. Per θ = 0 si ottiene il metodo di Eulero esplicito, per θ = 1 si ha il metodo di Eulero implicito. Crank-Nicolson è il caso θ = 0.5.

Alvise Sommariva  Equazione del calore.  7/ 24

Esperimento numerico.

Studiamo numericamente l’equazione del calore

\{ \begin{matrix}\frac{∂u}{∂t} = \frac{∂^{2}u}{∂x^{2}} + G, 0 < x < 1, t > 0 \\ u(0, t) = d_{0}(t), u(1, t) = d_{1}(t), t ≥ 0 \\ u(x, 0) = f (x), 0 ≤ x ≤ 1 \end{matrix}  (14)

per

\{ \begin{matrix}G(x, t) = (−0.1 + π^{2}) (exp(−0.1 · t) sin(π x)) \\ d_{0}(t) = 0, d_{1}(t) = 0 \\ f (x) = sin(π x) \end{matrix}  (15)

avente quale soluzione

u(x, t) = exp(−0.1 · t) sin(π x).

Alvise Sommariva  Equazione del calore.  8/ 24

Esperimento numerico.

Figura: Grafico della soluzione dell’equazione del calore (14).

Alvise Sommariva  Equazione del calore.  9/ 24

Esperimento numerico.

f u n c t i o n gt=g ( t , x , delta , d0 , d1 , G )
% CALCOLA
% g ( t ) =(1/hx ˆ 2) ∗(d0 ( t ) , 0 , …, 0 , d1 ( t ) ) +[G( x1 , t ) , …, G( x
(m−1) , t ) ] .
gt=f e v a l ( G , x , t ) ;
gt ( 1 )=gt ( 1 ) +(1/delta ˆ 2) ∗ f e v a l ( d0 , t ) ;
gt ( l e n g t h ( gt ) )=gt ( l e n g t h ( gt ) ) +(1/delta ˆ 2) ∗ f e v a l ( d1 , t ) ;

 f u n c t i o n [ V_hist , x_mid , t_mid , hx , ht]=cnheat ( theta , tfin , m
     , f , d0 , d1 , G , time_step_factor )

% INPUT .
% t h e t a  : PARAMETRO DEL THETA METODO (TRA 0 ED 1 ,
POSSIBILMENTE 0 ,1 INCLUSI ) .
% t f i n  : ISTANTE FINALE .
% m  : DETERMINA IL PASSO SPAZIALE CON PASSO h x
=1/m.
% f  : u ( x , 0 )=f ( x )
% d0 , d1 : u ( 0 , t )=d0 ( t ) , u ( 1 , t )=d1 ( t ) .
% G  : ( D t ) u=(Dˆ 2 x ) u+G .
% t i m e s t e p f a c t o r : IL PASSO TEMPORALE E ’
t i m e s t e p f a c t o r ∗(hx ˆ 2) / 2 .
% OUTPUT.
% V h i s t : SOLUZIONE ISTANTE PER ISTANTE . RIGHE
I N I Z I A L I −> BASSO t .
% x mid  : PUNTI INTERNI NELLA VARIABILE x .
% t mid  : PUNTI INTERNI NELLA VARIABILE t .
% hx  : PASSO SPAZIALE .
% ht  : PASSO TEMPORALE.
Alvise Sommariva  Equazione del calore.  11/ 24

hx=1/m ; % DELTA=hx

matorder=m−1;

ht=time_step_factor ∗(hx ˆ 2) / 2 ; % STEP TEMPORALE.
x =(0: hx : 1 ) ’ ;
x_mid=x ( 2 : l e n g t h ( x ) −1 ,:) ;  % PUNTI INTERNI .
u0=f e v a l ( f , x_mid ) ;  % VETTORE SOLUZIONE AL
TEMPO ”0”.
t =(0: ht : tfin ) ’ ;  % TEMPI DA ANALIZZARE .
t_mid=t ( 2 : l e n g t h ( t ) , 1 ) ;  % TEMPI t > 0 .
% COSTRUZIONE MATRICE LAMBDA.
submat=[ z e r o s ( 1 , matorder ) ; eye ( matorder −1) ( z e r o s ( 1 ,

    matorder −1) ) ’ ] ;
supmat=submat ’ ;
lambda_matrix =(1/hx ˆ 2) ∗( d i a g (−2∗ones ( m −1 ,1) )+submat+
    supmat ) ;

Alvise Sommariva  Equazione del calore.  12/ 24

i f theta < 1
     gt_prev=g ( 0 , x_mid , hx , d0 , d1 , G ) ;

end

V_old=u0 ;
V_hist=[V_old ’ ] ;
err_hist = [ ] ;

f o r index =2: l e n g t h ( t )
% CALCOLO u ’ ( t c u r r )=LAMBDA∗u ( t c u r r )+g ( t c u r r ) ,

t_curr=t ( index ) ;
gt_curr=g ( t_curr , x_mid , hx , d0 , d1 , G ) ;
switch theta
case 1
     A=eye ( s i z e ( lambda_matrix ) )−ht∗lambda_matrix ;
     b=V_old+ht∗gt_curr ; V_new=A\b ;
case 0
     V_new=V_old+ht ∗( lambda_matrix ∗V_old+gt_prev ) ;
     gt_prev=gt_curr ;

otherwise

     A=eye ( s i z e ( lambda_matrix ) )−(ht∗theta ) ∗ . . .
         lambda_matrix ;
     b=V_old+(ht∗(1−theta ) ) ∗lambda_matrix ∗V_old + . . .
         ( ht∗(1−theta ) ) ∗gt_prev+(ht∗theta ) ∗gt_curr ;
     V_new=A\b ; gt_prev=gt_curr ;

end

V_hist=[V_hist ; V_new ’ ] ; V_old=V_new ;

end

Alvise Sommariva  Equazione del calore.  14/ 24

1. Il parametro m determina il passo spaziale h_{x} = 1/m;

2. se theta=0 allora si utilizza il metodo di Eulero esplicito, mentre se theta=1 Eulero implicito;

3. la variabile tfin determina l’istante finale;

4. il parametro timestepfactor determina il passo temporale a partire da quello spaziale; se h_{tmax} = h_{x}^{2}/2 per h_{x} = 1/m il passo usato da Eulero esplicito è timestepfactor*htmax;

5. in seguito la demo valuta uno dei metodi per m = 2^{k} con k = 2, 3, 4, calcolando le ratio e_{2h_{x}}/e_{h_{x}} dove si è posto

e_{h_{x}} = max_{i}|v_{i}(tfin) − u_{x_{i},tfin})|

in cui v_{i}(t) := v^{(h_{x})}(x_{i}, tfin) è la soluzione ottenuta dal metodo scegliendo il parametro temporale uguale a h_{x} = 1/m;

6. il parametro mvect all’interno dello switch iniziale, determina

gli m da analizzare.

Alvise Sommariva  Equazione del calore.  15/ 24

% demoexample : 1 o 2 . t h e t a 0 : EE . ; 0 . 5 : CN. 1 : EI .

demoexample =2; theta =1; mvect =[4 8 16 32 6 4 ] ;
time_step_factor_vect=ones ( s i z e ( mvect ) ) ;
switch demoexample
case 1

tfin =10;  % TEMPO FINALE .
G=inline ( ’ (−0.1+ p i ˆ 2) ∗( exp (( −0.1) ∗t ) . ∗ s i n ( p i ∗x ) ) ’ , …
’ x ’ , ’ t ’ ) ;
d0=inline ( ’ z e r o s ( s i z e ( t ) ) ’ , ’ t ’ ) ;
d1=inline ( ’ exp (( −0.1) ∗t ) . ∗ s i n ( p i ) ’ , ’ t ’ ) ;
f=inline ( ’ s i n ( p i ∗x ) ’ , ’ x ’ ) ;
solution=inline ( ’ exp (( −0.1) ∗t ) . ∗ s i n ( p i ∗x ) ’ , ’ x ’ , ’ t ’ ) ;

case 2

tfin =0.1;  % TEMPO FINALE .
G=inline ( ’ z e r o s ( s i z e ( x ) ) ’ , ’ x ’ , ’ t ’ ) ;
d0=inline ( ’ z e r o s ( s i z e ( t ) ) ’ , ’ t ’ ) ;
d1=inline ( ’ z e r o s ( s i z e ( t ) ) ’ , ’ t ’ ) ;
f=inline ( ’ s i n ( p i ∗x ) ’ , ’ x ’ ) ;
solution=inline ( ’ exp ((− p i ˆ 2) ∗t ) . ∗ s i n ( p i ∗x ) ’ , ’ x ’ , ’ t ’ ) ;

end

Alvise Sommariva  Equazione del calore.  16/ 24

err_hist_prev_m = [ ] ;

f p r i n t f ( ’ \n \ t [THETA] : % 3 . 3 f [ TFIN ] : % 3 . 3 f ’ , theta , tfin ) ;
f o r mindex =1: l e n g t h ( mvect )
err_hist = [ ] ; m=mvect ( mindex ) ; % m IN UNA LISTA .

time_step_factor=time_step_factor_vect ( mindex ) ;
[ V_hist , x_mid , t_mid , hx , ht]=cnheat ( theta , tfin , m , f , d0 , d1
    , G , time_step_factor ) ;

[ X , Y]= meshgrid ( x_mid , t_mid ) ; % VALUTAZIONE ERRORI .
U=f e v a l ( solution , X , Y ) ; % SOL . ESATTA NELLA GRIGLIA .
% ERRORE COMPIUTO PER ”m” FISSATO .
err=norm ( U ( s i z e ( U , 1 ) , : )−V_hist ( s i z e ( V_hist , 1 ) , : ) , inf ) ;

err_hist=[err_hist ; err ] ;

f p r i n t f ( ’ \n \ t [m] : %3.0 f [ERROR ] : %2.2 e [ hx ] : %2.2 e [
ht ] : %2.2 e ’ , m , err , hx , ht ) ;

i f l e n g t h ( err_hist_prev_m ) > 0

f p r i n t f ( ’ [ RATIO ] : %2.2 f ’ , err_hist_prev_m ( s i z e (

     err_hist_prev_m , 1 ) ) /err ) ;

end

err_hist_prev_m=err_hist ;

end

Alvise Sommariva  Equazione del calore.  17/ 24

Esperimento numerico: Eulero esplicito.

Per motivi di stabilità tipici di Eulero esplicito, il passo temporale h_{t} deve essere inferiore o uguale a h_{x}^{2}/2. Vediamo su vari esempi cosa succede numericamente. Dopo aver settato in demoheat il parametro theta=0 scegliamo per esempio timestepfactor=1.5. Quindi dalla shell di Matlab/Octave digitiamo quanto segue

>> demoheatcn
[ THETA ] : 0.000 [ TFIN ] : 5.000 [ TIME STEP FACTOR ] : 1.50 e

+000

[ m ] :  4 [ ERROR ] : 1.40 e+004 [ hx ] : 2.50 e−001 [ ht ] : 4.69 e

−002

[ m ] :  8 [ ERROR ] : 7.32 e+100 [ hx ] : 1.25 e−001 [ ht ] : 1.17 e
−002 [ RATIO ] : 0.00
[ m ] : 16 [ ERROR ] : NaN [ hx ] : 6.25 e−002 [ ht ] : 2.93 e−003 [

    RATIO ] : NaN

>>

Esperimento numerico: Eulero esplicito.

Evidentemente bisogna scegliere un parametro timestepfactor più piccolo. Proviamo ad esempio timestepfactor=1.1.

>> demoheatcn

[ THETA ] : 0 . 0 0 0 [ TFIN ] : 5 . 0 0 0 [ TIME STEP FACTOR ] : 1 . 1 0 e+000
[ m ] : 4 [ ERROR ] : 3 . 2 5 e−002 [ hx ] : 2 . 5 0 e−001 [ ht ] : 3 . 4 4 e−002
[ m ] : 8 [ ERROR ] : 2 . 8 9 e+011 [ hx ] : 1 . 2 5 e−001 [ ht ] : 8 . 5 9 e−003
[ RATIO ] : 0 . 0 0
[ m ] : 16 [ ERROR ] : 3 . 7 7 e+149 [ hx ] : 6 . 2 5 e−002 [ ht ] : 2 . 1 5 e−003
[ RATIO ] : 0 . 0 0

>>

Esperimento numerico: Eulero esplicito.

Il metodo non fornisce evidentemente risultati apprezzabili. Scegliamo ora timestepfactor=1.0: il metodo di Eulero esplicito finalmente converge.

>> demoheatcn
[ THETA ] : 0.000 [ TFIN ] : 5.000 [ TIME STEP FACTOR ] : 1.00 e

+000

[ m ] :  4 [ ERROR ] : 3.25 e−002 [ hx ] : 2.50 e−001 [ ht ] : 3.13 e

−002

[ m ] :  8 [ ERROR ] : 7.93 e−003 [ hx ] : 1.25 e−001 [ ht ] : 7.81 e
−003 [ RATIO ] : 4.10
[ m ] : 16 [ ERROR ] : 1.97 e−003 [ hx ] : 6.25 e−002 [ ht ] : 1.95 e
−003 [ RATIO ] : 4.02

>>

In definitiva affinchè il metodo di Eulero esplicito converga, il passo temporale dev’essere scelto dell’ordine di h_{x}^{2}/2, che in molti casi risulta essere troppo piccolo e rende il metodo non competitivo dal punto di vista computazionale.

Alvise Sommariva  Equazione del calore.  20/ 24

Esperimento numerico: Eulero implicito.

Vediamo in questa sezione il comportamento di Eulero implicito. Osserviamo che a differenza di Eulero esplicito richiede la soluzione di sistemi lineari tridiagonali, ma ciò non è un problema dal punto di vista computazionale (il costo è di 5m per ogni t_{i}). Proviamo il comportamento per timestepfactor=1.5, dopo aver posto theta=1. Il metodo di Eulero implicito, a differenza di Eulero esplicito, converge. Infatti

>> demoheatcn

[ THETA ] : 1 . 0 0 0 [ TFIN ] : 5 . 0 0 0 [ TIME STEP FACTOR ] : 1 . 5 0 e+000
[ m ] : 4 [ ERROR ] : 3 . 2 6 e−002 [ hx ] : 2 . 5 0 e−001 [ ht ] : 4 . 6 9 e−002
[ m ] : 8 [ ERROR ] : 7 . 9 5 e−003 [ hx ] : 1 . 2 5 e−001 [ ht ] : 1 . 1 7 e−002
[ RATIO ] : 4 . 1 1
[ m ] : 1 6 [ ERROR ] : 1 . 9 7 e−003 [ hx ] : 6 . 2 5 e−002 [ ht ] : 2 . 9 3 e−003
[ RATIO ] : 4 . 0 3

>>

Esperimento numerico: Eulero implicito.

Per curiosità proviamo per timestepfactor=10, quindi con un passo temporale h_{t} relativamente grande, ottenendo

>> demoheatcn

[ THETA ] : 1 . 0 0 0 [ TFIN ] : 5 . 0 0 0 [ TIME STEP FACTOR ] : 1 . 0 0 e+001
[ m ] : 4 [ ERROR ] : 3 . 2 6 e−002 [ hx ] : 2 . 5 0 e−001 [ ht ] : 3 . 1 3 e−001
[ m ] : 8 [ ERROR ] : 7 . 9 6 e−003 [ hx ] : 1 . 2 5 e−001 [ ht ] : 7 . 8 1 e−002
[ RATIO ] : 4.10
[ m ] : 1 6 [ ERROR ] : 1 . 9 8 e−003 [ hx ] : 6 . 2 5 e−002 [ ht ] : 1 . 9 5 e−002
[ RATIO ] : 4.02

>>

Esperimento numerico: Eulero implicito.

Si può mostrare che in effetti, il metodo è A-stabile, e quindi non richiede alcun vincolo sullo step temporale. Ciò significa che per ogni valore di h_{x} e h_{t} la propagazione dell’errore avanzando nel tempo è sotto controllo o come si dice il metodo è incondizionatamente stabile. Inoltre se tanto la soluzione quanto d_{0}, d_{1}, g ed f sono sufficientemente regolari allora con ovvia notazione

max_{i,k}|u(x_{i}, t_{k}) − u_{i,k}| = O(h_{t} + h_{x}^{2}).

Esercizio.

Eseguire gli stessi test eseguiti per Eulero esplicito col θ-metodo, per θ = 0.25, θ = 0.5, θ = 0.75. Com’e’ il comportamento del metodo? Simile ad Eulero esplicito o ad Eulero esplicito?