Equazione del calore.

Alvise Sommariva

Università degli Studi di Padova

Dipartimento di Matematica

26 maggio 2014

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/ 36

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.

Equazione del calore.

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

u(t) := [u_{1}(t), …, u_{m−1}(t)]^{T}
u_{0} := [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}

A = \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) è equiv. al sistema di eq. differenziali (lineari)

u′(t) = Au(t) + g(t), u(0) = u_{0}  (5)
Alvise Sommariva  Equazione del calore.  4/ 36

Equazione del calore.

Nota.

Osserviamo che la matrice a predominanza diagonale

A = \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} ]  (6)

è ▶ A è simm. definita negativa per i teoremi di Gershgorin. ▶ Si mostra che se A ∈ ℝ^{(N−1)×(N−1)}, allora i suoi autov. sono

−\frac{2 − 2cos(\frac{jπ}{N})}{h_{x}^{2}}.

▶ Diagonalizzabile e I − h_{t}A è invertibile (A è def. negativa).

Alvise Sommariva  Equazione del calore.  5/ 36

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}  (7)

citiamo il metodo di Eulero esplicito (posto u_{n+1} = u(t_{n+1}))

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

e quello di Eulero implicito

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

Equazione del calore.

Nel nostro caso

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

e quindi il metodo di Eulero esplicito genera la successione

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

mentre Eulero implicito determina

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

o equivalentemente

\{ \begin{matrix}(I − h_{t}A)v_{n+1} = v_{n} + h_{t}g(t_{n+1}) \\ v_{0} assegnato \end{matrix}  (12)
Alvise Sommariva  Equazione del calore.  7/ 36

Equazione del calore.

Nota.

Osserviamo che a differenza del metodo esplicito, poichè

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

ad ogni iterazione si richiede la soluzione di un’equazione (che nel nostro caso è lineare). Usando i primi due teoremi di Gerschgorin, si vede che (I − h_{t}A) è 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}(Av_{n} + g(t_{n})))
+θ (v_{n} + h_{t}(Av_{n+1} + g(t_{n+1})))  (14)

con v_{0} assegnato. Per θ = 0 si ottiene il metodo di Eulero esplicito mentre per θ = 1 si ottiene il metodo di Eulero implicito.

Alvise Sommariva  Equazione del calore.  8/ 36

Stabilità

Consideriamo l’equazione del calore (con G ≡ 0)

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

Si dimostra che la soluzione esatta tende a 0 e quindi si richiede che la soluzione numerica abbia la stessa proprietà . In tal caso la soluzione si dice assolutamente stabile.

Stabilità

Il problema continuo viene discretizzato come proposto precedentemente, ottenendo il sistema differenziale

u′(t) = Au(t), u(0) = u_{0}  (16)

che è possibile risolvere con un generico θ-metodo

v_{n+1} = (1 − θ) (v_{n} + h_{t}Av_{n}) + θ (v_{n} + h_{t}Av_{n+1})  (17)

con v_{0} assegnato.

Alvise Sommariva  Equazione del calore.  10/ 36

Stabilità Eulero esplicito

Nel caso di Eulero esplicito abbiamo

v_{n+1} = v_{n} + h_{t}Av_{n} = (I + h_{t}A)v_{n}  (18)

Siccome A è diagonalizzabile, abbiamo che per qualche S, A = S^{−1}ΛS, con Λ diagonale. Posto z_{n} = Sv_{n},

z_{n+1} = Sv_{n+1} = S(v_{n} + h_{t}Av_{n}) = z_{n} + h_{t}SAv_{n}
= z_{n} + h_{t}SS^{−1}ΛSv_{n}
= z_{n} + h_{t}Λz_{n} = (I + h_{t}Λ)z_{n}  (19)

Si noti che z_{n} = Sv_{n} → 0 implica v_{n} → 0 in quanto

∥v_{n}∥ = ∥S^{−1}z_{n} ≤ ∥S^{−1}∥∥z_{n}∥ → 0.

Alvise Sommariva  Equazione del calore.  11/ 36

Stabilità Eulero esplicito

Se z_{n} = (u_{n,i})_{i}, Λ_{i,i} = λ_{i}, essendo Λ diagonale, da

z_{n+1} = (I + h_{t}Λ)z_{n}  (20)

la i-sima equazione diventa

z_{n+1,i} = (1 + h_{t}λ_{i})z_{n,i}.

che come nel caso di Eulero esplicito scalare, converge asintoticamente a 0 se e solo se

|1 + h_{t}λ_{i}| < 1.

Alvise Sommariva  Equazione del calore.  12/ 36

Stabilità Eulero esplicito

La condizione |1 + h_{t}λ_{i}| < 1 pone i vincoli sul passo

h_{t} ≤ \frac{2}{|λ_{i}|}, i = 1, …, N − 1.

Essendo

λ_{i} = −\frac{2 − 2cos(\frac{jπ}{N})}{h_{x}^{2}}, j = 1, …, N − 1

abbiamo che λ_{i} < 0 per i = 1, …, N − 1 e

| min_{i}(λ_{i})| = |λ_{N−1}| = | − \frac{2 − 2cos(\frac{(N−1)π}{N})}{h_{x}^{2}}| ≈ \frac{| − 4|}{h_{x}^{2}} = \frac{4}{h_{x}^{2}}

da cui la condizione sulla stabilità asintotica

h_{t} ≤ \frac{h_{x}^{2}}{2} .

Alvise Sommariva  Equazione del calore.  13/ 36

Stabilità Eulero implicito

Nel caso di Eulero implicito abbiamo

v_{n+1} = v_{n} + h_{t}Av_{n+1}  (21)

da cui immediatamente

v_{n+1} = (I − h_{t}A)^{−1}v_{n}  (22)

Siccome A è diagonalizzabile, abbiamo che per qualche S, A = S^{−1}ΛS, con Λ diagonale. Posto z_{n} = Sv_{n},

z_{n+1} = Sv_{n+1} = S(I − h_{t}A)^{−1}v_{n} = S(S^{−1}S − h_{t}S^{−1}ΛS)^{−1}v_{n}
= S(S^{−1}(I − h_{t}Λ)S)^{−1}v_{n} = SS^{−1}(I − h_{t}Λ)^{−1}Sv_{n}
= (I − h_{t}Λ)^{−1}z_{n}  (23)

Se z_{n} = (z_{n,i})_{i}, Λ_{i,i} = λ_{i}, la i-sima equazione diventa

z_{n+1,i} = (1 − h_{t}λ_{i})^{−1}z_{n,i}.

che quindi converge converge asintoticamente a 0 se e solo se

|1/(1 − h_{t}λ_{i})| < 1.

Alvise Sommariva  Equazione del calore.  14/ 36

Stabilità Eulero implicito

Siccoma A è definita negativa e simile a Λ, abbiamo che λ_{i} < 0 e quindi

|1/(1 − h_{t}λ_{i})| < 1.

è verificata per qualsiasi scelta di h_{t}.

Questo significa che qualsiasi sia h_{t}, la successione del metodo di Eulero implicito converge asintoticamente a 0, come la soluzione dell’equazione del calore.

Stabilità Crank-Nicolson

Nel caso di Crank-Nicolson abbiamo

v_{n+1} = (1/2)(v_{n} + h_{t}Av_{n+1}) + (1/2)(v_{n} + h_{t}Av_{n})  (24)

da cui immediatamente

v_{n+1} = v_{n} + (1/2)h_{t}Av_{n+1} + (1/2)h_{t}Av_{n}  (25)

cioè

(1 − \frac{h_{t}}{2} A)v_{n+1} = (1 + \frac{h_{t}}{2} A)v_{n}.  (26)

Siccome A è diagonalizzabile, abbiamo che per qualche S, A = S^{−1}ΛS, con Λ diagonale. Posto z_{n} = Sv_{n},

S(1 − \frac{h_{t}}{2} A)v_{n+1} = S(1 − \frac{h_{t}}{2} A)S^{−1}Sv_{n+1} = (1 − \frac{h_{t}}{2} Λ)z_{n+1} (27)

e pure

S(1 + \frac{h_{t}}{2} A)v_{n} = S(1 + \frac{h_{t}}{2} A)S^{−1}Sv_{n} = (1 + \frac{h_{t}}{2} Λ)z_{n}  (28)
Alvise Sommariva  Equazione del calore.  16/ 36

Stabilità Crank-Nicolson

Da

(1 − \frac{h_{t}}{2} A)v_{n+1} = (1 + \frac{h_{t}}{2} A)v_{n},  (29)
S(1 − \frac{h_{t}}{2} A)v_{n+1} = S(1 − \frac{h_{t}}{2} A)S^{−1}Sv_{n+1} = (1 − \frac{h_{t}}{2} Λ)z_{n+1} (30)

e

S(1 + \frac{h_{t}}{2} A)v_{n} = S(1 + \frac{h_{t}}{2} A)S^{−1}Sv_{n} = (1 + \frac{h_{t}}{2} Λ)z_{n}  (31)

moltiplicando a sinistra ambo i membri di (29) per S, abbiamo

(1 − \frac{h_{t}}{2} Λ)z_{n+1} = (1 + \frac{h_{t}}{2} Λ)z_{n}

da cui

z_{n+1} = (1 − \frac{h_{t}}{2} Λ)^{−1}(1 + \frac{h_{t}}{2} Λ)z_{n}
Alvise Sommariva  Equazione del calore.  17/ 36

Stabilità Crank-Nicolson

La i-sima equazione diventa

z_{n+1,i} = (1 − \frac{h_{t}}{2} λ_{i})^{−1}(1 + \frac{h_{t}}{2} λ_{i})z_{n,i}

e come nel caso scalare converge a 0 se e solo se

\frac{|1 + \frac{h_{t}}{2} λ_{i}|}{|1 − \frac{h_{t}}{2} λ_{i}|} < 1

che è sempre verificata poichè λ_{i} < 0.

Questo significa che qualsiasi sia h_{t}, la successione del metodo di Crank-Nicolson converge asintoticamente a 0, come la soluzione dell’equazione del calore.

Nota.

Si dimostra che le matrici dei sistemi lineari dei metodi di Eulero implicito e di Crank-Nicolson sono molto meglio condizionate della matrice A. Infatti si dimostra che

λmin(A) ≈ −4/h_{x}^{2}

ed essendo 1 − cos(x) ≈ x^{2}/2 per x ≈ 0, Nh_{x} = 1,

λmax(A) = −\frac{2(1 − cos(π/N))}{h_{x}^{2}} ≈ −\frac{2(π/N)^{2}}{2h_{x}^{2}} = −π^{2}

da cui cond_{2}(A) = \frac{4}{π^{2}h_{x}^{2}} mentre

cond_{2}(1 − h_{t}A) = \frac{1 − h_{t}λmin(A)}{1 − h_{t}λmax(A)} ≈ \frac{1 + 4h_{t}/h_{x}^{2}}{1 + h_{t}π^{2}} ≈ 1 + 4h_{t}/h_{x}^{2}.

Alvise Sommariva  Equazione del calore.  19/ 36

Nota.

Quanto visto, presuppone che nell’equazione del calore sia g ≡ 0, ma una analisi simile può essere effettuata anche per il caso in cui g ≠ 0.

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}  (32)

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}  (33)

avente quale soluzione

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

Alvise Sommariva  Equazione del calore.  21/ 36

Esperimento numerico.

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

Alvise Sommariva  Equazione del calore.  22/ 36

Esperimento numerico.

 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 I L 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 : I L 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.
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)) ’ ] ) ’ ;

lambda_matrix =(1/hx ˆ 2) ∗( d i a g (−2∗ones ( m −1,1) )+submat+supmat ) ;
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

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 ) ;

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.  25/ 36

% 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

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

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 . 5 0 e+000
[ m ] :  4 [ ERROR ] : 1 . 4 0 e+004 [ hx ] : 2 . 5 0 e−001 [ ht ] : 4 . 6 9 e−002
[ m ] :  8 [ ERROR ] : 7 . 3 2 e+100 [ hx ] : 1 . 2 5 e−001 [ ht ] : 1 . 1 7 e−002 [ RATIO ] : 0 . 0 0
[ m ] : 16 [ ERROR ] : NaN [ hx ] : 6 . 2 5 e−002 [ ht ] : 2 . 9 3 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 . 0 0 e+000
[ m ] :  4 [ ERROR ] : 3 . 2 5 e−002 [ hx ] : 2 . 5 0 e−001 [ ht ] : 3 . 1 3 e−002
[ m ] :  8 [ ERROR ] : 7 . 9 3 e−003 [ hx ] : 1 . 2 5 e−001 [ ht ] : 7 . 8 1 e−003 [ RATIO ] : 4 . 1 0
[ m ] : 16 [ ERROR ] : 1 . 9 7 e−003 [ hx ] : 6 . 2 5 e−002 [ ht ] : 1 . 9 5 e−003 [ RATIO ] : 4 . 0 2

>>

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.

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 . 1 0
[ m ] : 1 6 [ ERROR ] : 1 . 9 8 e−003 [ hx ] : 6 . 2 5 e−002 [ ht ] : 1 . 9 5 e−002
[ RATIO ] : 4 . 0 2

>>

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?

Esercizio (dal Quarteroni-Saleri).

Consideriamo lequazione del calore \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

con

G(x, t) = − sin(x) sin(t) + sin(x) cos(t) u(x, 0) = sin(x)

e condizioni al bordo

u(0, t) = 0,

u(1, t) = sin(1) cos(t).

Alvise Sommariva  Equazione del calore.  35/ 36

Esercizio.

In questo caso la soluzione

u(x, t) = sin(x) cos(t).

Si confrontino via plot semi-logaritmico il comportamento degli errori

max_{i=0,...,N} |u(x_{i}, 1) − u_{i}^{M}|

rispetto al passo temporale su una griglia uniforme lungo la direzione spaziale con h_{x} = 0.002, h_{x} = 0.004, h_{x} = 0.008 e su una griglia temporale uniforme con h_{t} = 10^{−3} + 2k · 10^{−3}, k = 0, …, 10, per un θ-metodo con θ = 0, θ = 0.5, θ = 1.

I valori {u_{i}^{M}} rappresentano la soluzione alle differenze finite calcolata al tempo t^{M} = 1.