aa
9 10
lezione di laboratorio
Laurea Ingegneria CIVILE
Lauree Specialistiche in Ingegneria
CHIMICA, ELETTRONICA,
AMBIENTE
a.a. 2007-2008
1
Generalità su un problema del 2° ordine
di tipo iperbolico
 utt  v 2 uxx x   x0 , x N  t  t 0

u  x , t 0     x  x   x0 , x N 

(1) 
 ut  x , t0     x  x   x0 , x N 
u x , t  g t , u x , t  g t t  t
 N  2 
1 
0
  0 
Si fa notare che le funzioni al bordo sono non nulle.
Si può assumere per semplicità:
x0  0, t0  0.
2
Approssimazioni utilizzate
Si approssimano le derivate parziali seconde con le
differenze finite:
ui , j 1  2ui , j  ui , j 1
 2u

,
2
2
t  x ,t 
k
i j
ui 1, j  2ui , j  ui 1, j
 2u

2
x  x ,t 
h2
i j
con errore locale nell’approssimazione delle derivate:
  , O h 
O k
2
2
La condizione iniziale sulla derivata si approssima con
una differenza centrale:
ui ,1  ui ,1
u

   xi   errore locale: O k 2
t  xi ,t0 
2k
 
3

vk
,
h
Problema discreto
Si indica
si considera dapprima j = 0, si
utilizzano le condizioni iniziali e si ottiene la soluzione
approssimata per t  t1 . Per j  1 ci serviremo della
seconda equazione del sistema (2) e delle condizioni
iniziali/al contorno:
1 2

2
u






1


 i  k  i i  1, , N  1
 i ,1 2  i -1 i 1 

u   2u  2 1   2 u   2u  u
, N 1 j  1
i 1, j
i, j
i 1, j
i , j 1 i  1,
 2   i , j 1
u  
i  0,1, , N
i ,0
i

 u0, j  g1 j uN , j  g2 j j  0




Si è posto:  i    xi  ;  i    xi  ; g1 j  g1  t j  ; g2 j  g2  t j 
4
Costruzione della forma vettoriale
Dalla seconda equazione dello schema assumendo ad
esempio N=5 e j=1, si ha:
Vj(1)
u1, 2

u 2, 2

u 3, 2

u4, 2

 2 1  
 2 1  
 2 1  

u
u
u
Uj-1
  2 u0,1  2 1   2 u1,1   2 u 2,1  u1,0
  2 u1,1
  2 u 2,1
  2 u 3,1
2
2
2
2


u 3,1  u 2,0
2,1
2


u4,1  u 3,0
3,1
2


u5,1  u4,0
4,1
Vj (N-1)
il vettore a primo membro è la soluzione approssimata nei
punti al livello 2.
5
Costruzione della forma vettoriale
Poniamo quindi:
T
U1   ui ,1 
i  1,
, N 1
Per j  1 :
V j   u0, j , 0,
U 0   1 ,  2 ,  ,  N 1 
T
T
1 2

2
    i -1   i 1   1    i  k  i 
2
 i 1,

U j   u1, j , u2, j ,
 

T

, uN 1, j 
, 0, uN , j    g1 t j , 0,
T
, N 1
 
, 0, g2 t j 
N.B. Il vettore V j ha N-1 componenti di cui
N–3 uguali a zero!!!
T
6
Forma vettoriale
U j 1  AU j  U j 1   V j
2
Si ottiene allora:
con la matrice A tridiagonale e simmetrica:
2 1   2  2

0


 2

2
A

2
I


T

2



2
2 

0

2
1




I matrice identità, T matrice tridiagonale con 2 sulla
diagonale principale e –1 sulle due codiagonali



Si ricorda che il metodo converge se

 1
7
Function PDE_iperboliche:parte 1
function [x,t,sol]=PDE_iperboliche(t0,M,x0,xN,h,k,v,r,f,l,
g1,g2);
alfa=v*k/h;x=(x0:h:xN)';x(end)=xN;N=length(x)-1;
f=eval(f).*ones(size(x)); %condizione iniziale U(x,t0)
%La presenza del vettore è necessaria per il caso f=cost.
U0=f(2:N);
vv=eval(l).*ones(size(x)); %condizione iniziale Ut(x,t0)
% Calcolo della soluzione al passo 1
t=t0;x=x(2:end-1); tnoto=eval(r).*ones(size(x));
U1=alfa^2/2*(f(1:N-1)+f(1+2:N+1))+(1alfa^2)*f(2:N)+k*vv(2:N) + 0.5*k^2*tnoto;
tM=k*M+t0;t=linspace(t0,tM,M+1)';
v1=eval(g1).*ones(size(t)); %condizione al contorno U(xo,t)
v2=eval(g2).*ones(size(t)); %condizione al contorno U(xN,t)
Vj=zeros(N-1,1);
sol=[f'; [v1(2) U1' v2(2)]];Uj=U1;Uj_1=U0;
8
Function PDE_iperboliche:parte 2
T=2*eye(N-1)-diag(ones(1,N-2),-1)-diag(ones(1,N-2),+1);
A=2*eye(N-1)-alfa^2*T;
t=t(2);
for j=2:M
Vj([1,N-1])=[v1(j),v2(j)];
tnoto=eval(r).*ones(size(x));
Uj1=A*Uj-Uj_1+alfa^2*Vj + k^2*tnoto;
sol=[sol;[v1(j+1); Uj1 ;v2(j+1)]'];
Uj_1=Uj;
Uj=Uj1;
t=t+k;
end
t=linspace(t0,tM,M+1)';
x=[x0;x;xN];
9
Esercizio 1
 utt  v 2 uxx x   0,10  t  0 v  1

2
x

5

x

5
 
u  x, 0  f  x  
e
x   0,10

0.76

 ut  x , 0   l  x   0 x   0,10

 u  0, t   g1  t   0 u  10, t   g2  t   0
Per le funzioni che forniscono le condizioni, si sono
utilizzati gli stessi nomi della function
PDE_iperboliche
10
Quesiti a, b
a) calcolare la soluzione approssimata del
problema con il metodo esplicito utilizzando i
seguenti valori per i passi k=0.25, e h=0.30 ed un
numero di intervalli temporali pari a M=40;
b) graficare la superficie che si ottiene e le sue
curve di livello in due figure distinte, utilizzando
i comandi surf, per la figura 1, e contour per
la figura 2. Analizzare bene e commentare le
figure.
11
Quesiti c, d
c) Cosa accade per la teoria se si modifica la
condizione sull’asse t in g1( t ) = 1 ?
d)
Risolvere ancora il problema, eseguire i grafici
come al punto b) ed evidenziarne le differenze col
caso precedente. Riportare inoltre sulla figura 2
(contour) anche la retta di equazione t = x/v; cosa
rappresenta?
Lo studente risolva ed analizzi il problema
assumendo M = 80, ed utilizzi il comando
subplot per le figure di cui al punto b).
12
Calcolo della soluzione: file onda.m
clear all;close all;clc
% Dati del problema
t0=0; M=40;
x0=0; xN=10;
k=.25;h=.30;
%cond.iniziale su t=0
f='(x-5).*exp(-(x-5).^2 )/.76 ';
l='0'; %altra cond.iniz.su t=0
g1='0'; %cond. al contorno su x=x0
g2='0'; %cond.al contorno su x=xN
r='0';
v=1;alfa=v*k/h;
% implementazione del metodo
[x,t,sol]=PDE_iperboliche(t0,M,x0,xN,h,k,v,r,f,l,g1,g2);
13
%
%
%
%
N.B. Non si costruisce la matrice dei nodi con
meshgrid perché non si valuta la soluzione vera;
l’assegnazione particolare di H1 è stata fatta per
avere una grafica più significativa
figure(1)
H1=surf(t,x,sol'); % vedi N.B.
xlabel('t','FontWeight','bold','Fontsize',12)
ylabel('x','FontWeight','bold','Fontsize',12)
title(['Soluzione con alfa= '
num2str(alfa)],'FontWeight','bold','Fontsize',12)
set(gca,'FontWeight','bold','Fontsize',12)
figure(2)
[C,H]=contour(x,t,sol,20); axis square
xlabel('x','FontWeight','bold','Fontsize',12)
ylabel('t','FontWeight','bold','Fontsize',12)
title('Curve di
livello','FontWeight','bold','Fontsize',12)
set(H,'LineWidth',2) %spessore della linea
set(gca,'FontWeight','bold','Fontsize',12)
14
Grafico della soluzione numerica
Figura 1
15
Altro grafico di sol
H2=surf(x,t,sol);% osservare bene
% l’istruzione!!!
16
Curve di livello
Figura 2
17
Grafico della condizione iniziale
f='(x-5).*exp(-(x-5).^2)/.76';
fplot(f, [0 10]), xlabel('x');ylabel('y');
title(['Grafico di f(x)=', f]);
18
Formula di D’Alembert
La soluzione analitica del problema di Cauchy
associato al problema (1), con gli stessi dati iniziali
è data, in base alla formula di D’Alembert, da:
1
1 x  vt

u( x , t )   f ( x  vt )  f ( x  vt )   l ( z )dz 
2
v x vt

con i dati assegnati su una linea che non è
caratteristica.
19
Soluzione analitica dell’esercizio 1
x  5  ( x  5)2 l ( x )  0,
e
,
Tenendo conto che: f ( x ) 
0.76
la soluzione analitica del problema di Cauchy
ottenuto dall’esercizio 1 non considerando le
condizioni al contorno, è:
1 
 ( x  t  5)2
 ( x  t  5)2 
u( x , t ) 
( x  t  5)e
 ( x  t  5)e

2  0.76 
20
Soluzione del problema di Cauchy:
formula di D’Alembert
Figura 3
[X,T]=meshgrid(x,t); sv=inline(f);
solvera=1/2*(sv(X+T)+sv(X-T));surf(t,x,solvera')
21
Curve di livello della soluzione di
D’Alembert
Osservare
le differenze
tra le figure 1
e 3 e tra le
figure 2 e 4;
giustificare
teoricamente
le differenze.
Figura 4
figure(4)
[C,H]=contour(x,t,solvera,20); axis square
22
Soluzione numerica: M=3*40
N.B. 1- Il numero degli intervalli temporali
è stato triplicato per evidenziare il fenomeno delle onde;
23
2- la figura è stata opportunamente ruotata
c) Caso con g1 ( t ) = 1
Imporre la condizione:
u  0, t   g1  t   1 t  0
equivale ad introdurre una discontinuità nel punto
(0,0), discontinuità che si propaga lungo la
caratteristica passante per tale punto di equazione
tx
. Dal punto di vista del calcolo, è
v
sufficiente modificare nel file onda.m solo
l’istruzione:
% condizione al contorno in (0,tj)
g1='0'sostituendola con g1='1';
24
Grafico della soluzione con g1 ( t )=1
[x,t,sol]=PDE_iperboliche(t0,M,x0,xN,h,k,v,r,f,l,g1,g2);
figure(1)
... % rimangono le stesse istruzioni
figure(2)
[C,H]=contour(x,t,sol,20);
... % rimangono le stesse istruzioni
set(gca,'FontWeight','bold','Fontsize',12)
% grafico della retta di equazione t=x/v
t1=x/v; % caratteristica passante per (0,0)
hold on
H3=plot(x,t1,'k');
set(H3,'LineWidth',3) %spessore della linea
set(gca,'FontWeight','bold','Fontsize',12)
25
Grafico soluzione numerica in presenza
di singolarità
Figura 1
26
Altro grafico di sol in presenza di
singolarità
H2=surf(x,t,sol);
27
Curve di livello senza la retta t = x/v
Figura 2
28
Curve di livello e retta caratteristica per
(0,0)
29
Generalità su un problema del 2°
ordine di tipo parabolico
 ut  c uxx x   x0 , x N  t  t 0

 u  x , t 0     x  x   x0 , x N 

t  t0
 u  x0 , t   g  t 
u x , t  l t
t  t0
  N   
N.B. Condizioni al contorno di Dirichlet
Si può assumere per semplicità: x0  0, t0  0.
30
Metodo di Crank-Nicolson
Si colloca l’equazione differenziale in

uxx xi , t j 1 2

 x ,t .
j 1 2
i
è approssimata con la media delle differenze
centrali relative ai livelli j e j+1:


uxx xi , t j 1 2 
ui 1, j  2ui , j  ui 1, j  ui 1, j 1  2ui , j 1  ui 1, j 1
2 h2
La ut  xi , t j 1 2  è approssimata con una
differenza centrale di passo k/2:
j+1
j 1 2
ut  xi , t j 1 2  
ui , j 1  ui , j
k
j
i-1
i
i+1
31
Schema di Crank-Nicolson
c
ui 1, j  2ui , j  ui 1, j  ui 1, j 1  2ui , j 1  ui 1, j 1
2h
2

ui , j 1  ui , j
k
k
Posto   c 2 , si ottiene:
h
  u i 1, j 1  2 1    u i , j 1   u i 1, j 1   u i 1, j  2 1    u i , j   u i 1, j

i  1, 2 , N  1, j  0


i  0,1, , N
 ui ,0   i

j0
 u0, j  g j
u  l
j0
j
 N,j
32
Costruzione della forma vettoriale
Se si assume N=5 , dalla prima equazione dello
schema per j=0, si ottiene il sistema:
Vj+1(1)
  u0,1  2  1    u1,1   u2,1

  u1,1  2  1    u2,1   u3,1

  u2,1  2  1    u3,1   u4,1

  u3,1  2  1    u4,1   u5,1
Vj(1)
  u0,0  2  1    u1,0   u2,0
  u1,0  2  1    u2,0   u3,0
  u2,0  2  1    u3,0   u4,0
  u3,0  2  1    u4,0   u5,0
Vj+1(N-1)
Vj(N-1)
33
Forma vettoriale
U 0   1 ,  2 ,  ,  N 1 
Per j = 0 abbiamo :
T
Per j=1,2,…poniamo:
U j   u1, j , u2, j ,
V j   u0, j , 0,
 
, 0, uN , j    g t j , 0,

T
, uN 1, j 
 
, 0, l t j 

T
T
N-3
N-3
Si ottiene il sistema lineare:
AU j 1  BU j   V j 1  V j  , j  0,1, ...
La quantità nell’ovale è b
34
Caratteristiche delle matrici A, B
A,
N 1 N 1

R
B
sono tridiagonali e simmetriche
 2  1    

0




  2I  T ;
A

 


0
 2  1    

 2 1    



B


 0


  2I  T .
 

2 1    
0

Si ricorda che il metodo converge
 h, k .
35
Function PDE_paraboliche_CN:parte 1
function [x,t,sol]=PDE_paraboliche_CN(t0,M,x0,xN,h,k,c,r,
f,g,l);
alfa=k*c/h^2;x=(x0:h:xN)'; x(end)=xN; N=length(x)-1;
tM=M*k+t0; t=linspace(t0,tM,M+1)';
f=eval(f).*ones(size(x)); %condizione iniziale U(x,0)
v1=eval(g).*ones(size(t)); %condizione al contorno U(xo,t)
v2=eval(l).*ones(size(t)); %condizione al contorno U(xN,t)
U0=f(2:N); sol=f';
Vj=zeros(N-1,1);
Uj=U0;
bb=alfa*ones(N-2,1);
A=-diag(bb,-1)+2*(1+alfa)*eye(N-1)-diag(bb,1);
B= diag(bb,-1)+2*(1-alfa)*eye(N-1)+diag(bb,1);
N.B. Questa function risolve un problema parabolico del tipo:
ut  cuxx  r ( x , t ).
36
Function PDE_paraboliche_CN:parte 2
t=t0+k/2;x=x(2:end-1);
Attenzione all’indice
for j=1:M
tnoto=eval(r).*ones(N-1,1);
Vj([1,N-1])=[v1(j)+v1(j+1),v2(j)+v2(j+1)];
b=B*Uj+alfa*Vj + 2*k*tnoto;
Uj1=A\b;
sol=[sol;[v1(j+1); Uj1 ;v2(j+1)]'];
Uj=Uj1;
t=t+k;
end
t=linspace(t0,tM,M+1)';x=[x0;x;xN];
j!!!
37
Esercizio 2
Si risolva con il metodo di Crank-Nicolson il seguente
problema con le condizioni al contorno miste
Dirichlet-Neumann:
 ut  uxx x   0,1 t  0

u  x, 0  f  x   x

u
0,
t

g
t

0
u
1,
t

l
t

t








x

con passo spaziale h = 0.2, passi temporali k=[0.5,
0.25,0.125] e numero di intervalli temporali M=30.
1- Si rappresenti la soluzione e le sue linee di livello.
2- Si calcoli la soluzione approssimata in x=0.4 e t
variabile, nel caso M fissato e nel caso tM fissato.
38
Condizioni iniziali ed al contorno
Si noti la continuità del problema nel punto (x0,t0)=(0,0):
u  x , 0  x  0  u  0, t  t  0  0.
Per conservare l’ordine di approssimazione, la
condizione ux  x N , t   l  t  viene approssimata con
una differenza centrale:


u x x N , t j 1 
uN 1, j 1  uN 1, j 1
2h
 l j 1
j  0.
39
Approssimazione del problema in esame
k
Posto   h2 , si assume:
quindi:
x0  0, x N  1, si ottiene
  ui 1, j 1  2 1    ui , j 1   ui 1, j 1   ui 1, j  2 1    ui , j   ui 1, j

i  0,1, , N
 ui ,0  xi

 u0, j  0
j0

 uN 1, j 1  uN 1, j 1
 t j 1 j  0
*

2h
Dall’uguaglianza (*) si ricava:
uN 1, j 1  2ht j 1 + uN 1, j 1
j0
40
Costruzione del sistema
Per i=1,2,3,4 e j = 0, lo schema diventa:
Vj(1)
Vj+1(1)
  u0,1  2  1    u1,1   u2,1   u0,0  2  1    u1,0   u2,0

  u1,1  2  1    u2,1   u3,1   u1,0  2  1    u2,0   u3,0

  u2,1  2  1    u3,1   u4,1   u2,0  2  1    u3,0   u4,0

 2 u3,1  2  1    u4,1  2h l1  2 u3,0  2  1    u4,0  2h l0
Vj+1(N)
N.B In questo caso la dimensione di V j e del
sistema è N; occorre infatti calcolare la soluzione
anche in x=xN.
Vj(N)
41
Sistema relativo al problema in esame

AU j 1  BU j   V j 1  V j
Il sistema diventa:
dove:
U j   u1, j , u2, j ,
V j   u0, j , 0,
 2  1    



A


0

, 0, 2 h l j 
T
, uN , j  ,
T
V j 1   u0, j 1 , 0,
, 0, 2h l j 1 

 2 1    


 B 

 


0

2  1    
0
2

T






2  1    
0
2
42
Soluzione del problema
Il sistema è lineare ad ogni livello j; la matrice dei
coefficienti di dimensione NxN ha la forma:
 2  1    



A


0



;
 

2  1    
0
2
A è matrice diagonalmente dominante ( A non singolare ),
quindi la soluzione approssimata nei nodi, per i=
1,2,3,4 ad ogni livello, esiste ed è unica.
43
Function Pde_paraboliche_CN1
...
U0=f(2:N+1);sol=f';Uj=U0;Vj=zeros(N,1);
bb=alfa*ones(N-1,1);
A=-diag(bb,-1)+2*(1+alfa)*eye(N)-diag(bb,1);
A(N,N-1)=-2*alfa;
B= diag(bb,-1)+2*(1-alfa)*eye(N)+diag(bb,1);
B(N,N-1)= 2*alfa;
t=t0+k/2;x=x(2:end);
for j=1:M
tnoto=eval(r).*ones(N,1);
Vj([1,N])=[v1(j)+v1(j+1),2*h*(v2(j)+v2(j+1))];
b=B*Uj+alfa*Vj +2*k*tnoto;
Uj1=A\b;
sol=[sol;[v1(j+1); Uj1]'];
Uj=Uj1;
t=t+k;
end
t=linspace(t0,tM,M+1)';x=[x0;x];
44
Inizializzazione dei dati
clear all;close all;
f='x'; %condizione iniziale nei nodi (xi,0)
g='0'; %condizione al contorno U(xo,t)
l='t'; %condizione al contorno Ux(xN,t)
r='0'; %termine noto
x0=0;xN=1;
t0=0;
h=.2;c=1;
xsol=0.4;
ind=round((xsol-x0)/h)+1;
tab=[];
k=[.5 .25 .125]; % attenzione si richiede che i valori di k siano
% multipli altrimenti la soluzione non può essere confrontata
% nei nodi corrispondenti
45
a - Calcolo della soluzione approssimata
con M fissato
% si può usare il comando sol(:,:,i)solo se le
% matrici hanno le stesse dimensioni
M=30;
tM3=t0+M*k(3);
num=fix((tM3-t0)/k(1));
t_fin=t0+num*k(1); % 3.5 valore limite comune
for i=1:length(k)
[x,t,sol(:,:,i)]=PDE_paraboliche_CN1(t0,M,x0,xN,
h,k(i),c,r,f,g,l);
n=round((t_fin-t0)/k(i))+1;
tab=[tab sol(1:round(k(1)/k(i)):n,ind,i)];
end
tab=[t(1:round(k(1)/k(i)):n) tab];
fprintf('
t
sol1
sol2
sol3\n')
46
fprintf('%7.4f %10.6f %10.6f %10.6f\n',tab')
Soluzione approssimata U(0.4,t)
t
0.0000
0.5000
1.0000
1.5000
2.0000
2.5000
3.0000
3.5000
sol1
0.400000
0.225572
0.212672
0.438797
0.597636
0.815526
1.004790
1.209085
sol2
0.400000
0.187795
0.262851
0.425218
0.613716
0.810050
1.008788
1.208307
sol3
0.400000
0.206611
0.265497
0.424684
0.612849
0.809411
1.008411
1.208120
47
Rappresentazione della soluzione
%____________________________________________________
%
% Queste istruzioni inserite nel ciclo consentono di
% graficare tutte le soluzioni e relative curve di
% livello.
%____________________________________________________
jj=1; alfa=k(i)*c/h^2;
figure(jj)
surf(x,t,sol(:,:,i)); colorbar('vert')
title(['Soluzione con alfa= ' num2str(alfa)])
xlabel('x');ylabel('t')
figure(jj+1)
contour(x,t,sol(:,:,i),20); colorbar('vert')
title(['Curve di livello con alfa= ' num2str(alfa)])
xlabel('x');ylabel('t')
48
jj=jj+2;
Grafico della soluzione approssimata con
M fissato
Soluzione con alfa= 3.125
3
4
2.5
3
2
2
1.5
1
1
0
4
1
2
Figura 1
t
0.5
0.5
0 0
0
x
49
Curve di livello con alfa= 3.125
Figura 2
3.5
3
3
2.5
2.5
2
t
2
1.5
1.5
1
1
0.5
0
0
0.5
0.2
0.4
0.6
0.8
1
x
Si invitano gli studenti a risolvere l’esercizio 2 con il
metodo esplicito, modificando opportunamente la
function relativa a tale metodo.
50
b - Calcolo della soluzione approssimata
con tM fissato
Si riportano ora le istruzioni relative al caso tM fissato.
tM=3.5; tab=[];
for i=1:length(k)
M=round((tM-t0)/k(i));
[x,t,sol]=PDE_paraboliche_CN1(t0,M,x0,xN,h,k(i),
c,r,f,g,l);
tab=[tab sol(1:round(k(1)/k(i)):end,ind)];
end
tab=[t(1:round(k(1)/k(i)):end) tab];
fprintf('
t
sol1
sol2
sol3\n')
fprintf('%7.4f %10.6f %10.6f %10.6f\n',tab')
51
Esercizio 3
Sia dato il seguente problema alle derivate parziali di
tipo parabolico:
1

 ut  2 u xx  3(1  x )
 u( x ,0)  x 3

 u(0, t )  3t

 u( 2, t )  8  3t
con soluzione vera:
0  x  2, t  0
0 x2
t0
t  0.
(1) u( x, t )  x  3t .
3
52
Quesito 1) e 2)
1) Si verifichi che la funzione (1) è soluzione del problema
proposto.
2) Si valuti, per il passo spaziale h=0.1, passo temporale
k=0.02 e M=20, l’errore assoluto massimo che si
commette usando il metodo di Crank-Nicolson ed il
metodo alle differenze esplicito. Si confrontino e
commentino i risultati; soddisfano le aspettative
teoriche?
Lo studente risolva ed analizzi il problema anche con M =
40, facendo le dovute considerazioni sull’errore di
troncamento delle approssimazioni effettuate, e su come
una perturbazione dovuta alla macchina potrebbe
influenzare i risultati.
53
Quesito 3)
3) Si costruiscano due tabelle che riportino
l’intestazione: x sol1 sol2 err1 err2
con le quantità x, sol1, sol2, err1, err2 rappresentanti,
rispettivamente i nodi spaziali, la soluzione numerica e
l’errore ottenuti con i due metodi, da riportare uno ogni
due, valutati in corrispondenza dei valori t=0.1 e t=0.4
utilizzando i seguenti formati di stampa:
3 cifre decimali e formato virgola fissa
per i nodi,
8 cifre decimali e formato esponenziale
per la soluzione,
2 cifre decimali e formato virgola mobile
per l’errore.
54
Istruzioni relative al quesito 2)
clc; clear all
t0=0;M=20;x0=0;xN=2;h=0.1; k=0.02;c=1/2;
r='-3*(1+x)';f='x.^3';g='-3*t';l='8-3*t';
% Implementazione dei metodi
[x,t,sol_CN]=PDE_paraboliche_CN(t0,M,x0,xN,h,k,c,r,f,g,l);
[x,t,sol_esp]=PDE_paraboliche(t0,M,x0,xN,h,k,c,r,f,g,l);
% Confronto con la soluzione vera
[X,T]=meshgrid(x,t);
Uvera=X.^3-3*T;
err_CN=abs(Uvera-sol_CN);
err_esp=abs(Uvera-sol_esp);
errmax_CN=max(max(err_CN)); errmax_esp=max(max(err_esp));
55
Istruzioni relative al quesito 3)
alfa=k*c/h^2;
disp(' alfa
err_CN
err_esp')
fprintf('\n %7.3f %11.2e %11.2e \n\n',
[alfa errmax_CN errmax_esp])
for t_val=[0.1 0.4]
i=round((t_val-t0)/k)+1;
tab=[x sol_CN(i,:)' sol_esp(i,:)', err_CN(i,:)',
err_esp(i,:)'];
tab_rid=tab(1:2:end,:);
fprintf([' \n\n Tabella per t=', num2str(t_val),
' \n\n x \t\t sol_CN \t\t sol_esp \t\t
err_CN \t\t err_esp \n'])
fprintf(' %7.3f
%16.8e %16.8e %10.2e %10.2e \n',
tab_rid')
end
56
Risultati quesito 3):tabella 1
Tabella per t=0.1
sol_esp
x
sol_CN
0.000
0.200
0.400
0.600
0.800
1.000
1.200
1.400
1.600
1.800
2.000
-3.00000000e-001
-2.92000000e-001
-2.36000000e-001
-8.40000000e-002
2.12000000e-001
7.00000000e-001
1.42800000e+000
2.44400000e+000
3.79600000e+000
5.53200000e+000
7.70000000e+000
Errori massimi
nei due metodi:
alfa
1.000
-3.00000000e-001
-2.92000000e-001
-2.36000000e-001
-8.40000000e-002
2.12000000e-001
7.00000000e-001
1.42800000e+000
2.44400000e+000
3.79600000e+000
5.53200000e+000
7.70000000e+000
err_CN
1.20e-014
err_CN
err_esp
0.00e+000
5.55e-017
2.78e-017
1.94e-016
4.16e-016
6.66e-016
8.88e-016
2.66e-015
7.11e-015
7.11e-015
0.00e+000
0.00e+000
1.67e-016
1.75e-015
3.43e-015
4.52e-015
1.73e-014
6.00e-015
3.33e-014
9.99e-014
1.10e-013
0.00e+000
err_esp
6.72e-007
57
Risultati quesito 3): tabella 2
x
0.000
0.200
0.400
0.600
0.800
1.000
1.200
1.400
1.600
1.800
2.000
sol_CN
Tabella per t=0.4
sol_esp
-1.20000000e+000
-1.19200000e+000
-1.13600000e+000
-9.84000000e-001
-6.88000000e-001
-2.00000000e-001
5.28000000e-001
1.54400000e+000
2.89600000e+000
4.63200000e+000
6.80000000e+000
-1.20000000e+000
-1.19200000e+000
-1.13600000e+000
-9.84000021e-001
-6.88000027e-001
-1.99999954e-001
5.28000251e-001
1.54400053e+000
2.89600067e+000
4.63200048e+000
6.80000000e+000
err_CN
err_esp
0.00e+000
6.66e-016
8.88e-016
1.44e-015
2.00e-015
3.39e-015
5.22e-015
7.77e-015
1.15e-014
8.88e-015
0.00e+000
0.00e+000
4.32e-009
1.91e-009
2.14e-008
2.72e-008
4.64e-008
2.51e-007
5.27e-007
6.72e-007
4.85e-007
0.00e+000
N.B. Il metodo di Crank-Nicolson riproduce ancora la soluzione (errore circa della
precisione di macchina). Il metodo esplicito fornisce errore più grande rispetto al
caso precedente; cosa succede per valori di t ancora più grandi?
58
Rappresentazione della soluzione
figure(1)
subplot(211), surf(x,t,Uvera);
title('Soluzione vera')
xlabel('x');ylabel('t')
subplot(212), surf(x,t,sol_esp);
title(['Soluzione metodo esplicito con alfa= '
num2str(alfa)])
xlabel('x');ylabel('t')
figure(2)
surf(x,t,err_esp);colorbar('vert')
title(['Errore metodo esplicito con alfa= ',
num2str(alfa)])
xlabel('x');ylabel('t')
59
Grafici della soluzione vera ed
approssimata
60
Grafico dell’errore col metodo esplicito
Si noti che per piccoli valori di t l’errore è piccolo; ma, essendo il
metodo instabile, all’aumentare del numero delle iterazioni e per
effetto dell’errore di macchina, amplificato dall’instabilità numerica,
61
l’errore diventa sempre più grande.