1

9a lezione di laboratorio

Laurea Specialistica in

Ingegneria MatematicaIngegneria dei Sistemi Energetici

Laurea Specialistica in

Ingegneria MatematicaIngegneria dei Sistemi Energetici

a.a. 2007-2008

2

Esercizio 1Dato il problema:

a) Dopo avere stabilito l’esistenza e l’unicità della soluzione, applicare il metodo delle differenze finite suddividendo l’intervallo in N=5 e N=10 sottointervalli;

0,2

0 1

12 2

y t y t t t

y

y

3

Esercizio 1: quesiti b, c

b) si calcoli la maggiorazione teorica dell’errore assoluto e si faccia il confronto con i risultati numerici ottenuti sapendo che la soluzione del problema è:

c) si applichi il metodo dello shooting e si analizzino i risultati che si ottengono facendo un confronto con i risultati ottenuti al punto a).

( ) cos sin ;y t t t t

4

PROBLEMA ai LIMITI: generalità

0 max

0 0

max max

, , ,

0* 1*

0* 1*

y t f t y t y t t t t

A a y t - a y t =

b y t + b y t =

0 max

0 0

max max

* * ,

0* ( ) - 1* ( )

0* ( ) 1* ( )

y t p t y t q t y t r t t t t

B a y t a y t

b y t b y t

PROBLEMA ai LIMITI LINEARE

5

Condizioni generali sufficienti per esistenza ed unicità della soluzione

0, , : ,f t y y C S a b R R

0', ,y yf f C S 1 20 , in y yf L f L S

0* 1 0 0* 1 0 0 0 0a a b b a b

TEOREMA A. IPOTESI:

TESI:

ESISTE UNICA la SOLUZIONE del PROBLEMA (A)

6

ESISTENZA, UNICITA’: caso di PROBLEMA LINEARE

Teorema dell’alternativa: Dato un problemalineare non omogeneo, la soluzione esiste unica se il problema omogeneo associato ammette un unica soluzione: quella nulla. In caso contrario il problema non omogeneo può non avere soluzione, qualora la soluzione esista, essa non è certamente unica.

Nel caso il problema sia lineare, si può, oltre al teorema precedente, utilizzare il seguente:

7

a) Esistenza ed unicità della soluzione del problema assegnato

Si può applicare il Principio dell’alternativa. Il problema omogeneo associato è:

0

0 0 02

z z

z z

la cui soluzione generale è:

1 2cos sinz t c t c t

da cui, per le condizioni ai limiti: ! 0.z t

La soluzione del problema non omogeneo è:

cos siny t t t t

8

Metodo alle differenze finite

Applicando il metodo alle differenze finite si approssima:

Il problema si trasforma allora in un sistema lineare attraverso le relazioni seguenti: Ax b

1 12

2i i ii

y y yy t

h

1 12

0

21,2, , 1

1 12

i i ii i

N

y y yy t i N

h

y y

1 12

0

21,2, , 1

1 12

i i ii i

N

y y yy t i N

h

y y

9

Considerazioni sul sistema lineare

T è matrice tridiagonale simmetrica con: T(i,i)=2 ed elementi codiagonali = -1

2

2

2

2

2 1 0 0

1 2 1 0

0 1

1

0 1 2

h

h

A h I T

h

2 2 2 21 2 2 11 1

2

T

N Nb h t h t h t h t

10

Costruzione della matrice A

%% Attenzione! La matrice dei coefficienti ha dimensione %% N-1xN-1t0=0;tmax=pi/2;det_A=[];for N=[5 10] %numero dei sottointervalli h=(tmax-t0)/N; v = ones(N-2,1); A =(-2+h^2)*eye(N-1)+diag(v,1)+diag(v,-1) detA=det(A)end

La soluzione è unica per entrambi i valori di N, essendo la matrice A non singolare.

detA = 3.2230 detA = -6.3859

11

a) Calcolo della soluzione: istruzioni

clear all; clct0=0;tmax=pi/2; alpha=1; beta=(pi/2)-1;N=[5 10];Yvera=inline('cos(x)-sin(x)+x');% metodo alle differenzep='0'; q='-1'; r='t';m=[1 0; 1 0];for n=N [T,Y]=Linear_ODE(t0,tmax,alpha,beta,n,m,p,q,r); Err=abs(Yvera(T)-Y); tab=[T Y Err]; fprintf(' Tempo soluzione \t Err\n') fprintf('%7.4f %15.8f %15.3e \n',tab')end

12

Metodo alle differenze: risultati per N = 5

Tempo soluzione Err1 0.0000 1.00000000 0.000e+000 0.3142 0.95657213 3.733e-004 0.6283 0.84974065 1.904e-004 0.9425 0.72105568 1.904e-004 1.2566 0.61422420 3.733e-004 1.5708 0.57079633 1.110e-016

13

Metodo alle differenze: risultati per N = 10

Tempo soluzione Err2 0.0000 1.00000000 0.000e+000 0.1571 0.98840160 6.809e-005 0.3142 0.95629116 9.237e-005 0.4712 0.90833674 8.182e-005 0.6283 0.84959737 4.710e-005 0.7854 0.78539816 1.221e-015 0.9425 0.72119896 4.710e-005 1.0996 0.66245958 8.182e-005 1.2566 0.61450517 9.237e-005 1.4137 0.58239472 6.809e-005 1.5708 0.57079633 1.110e-016

14

b) Limite teorico dell’errore assoluto

Vale la seguente relazione:

Poiché:

2

21

8 12IV

h

b ae h y

sin cos 1 0 in 0,2

y t t t t

sin cos 0 per , cos sin 04

y t t t t y t t t

minimo per 0 in 0, è decresc.4 2

t y y y

15

Quindi: Inoltre:

quindi:

sin cos 0 0, è decresc.2

V IVy t t t t y t

2

max 0 2 2 212.57*10

8 12h

t te h h

2

max 0 2 2 212.57*10

8 12h

t te h h

0 1.y y

max 0 , 12

IV IV IVy t y y

Allora:

16

Confronto maggiorazione teorica e numerica dell’errore assoluto

Norma1=norm(Err1,inf)Norma1 = 3.7334e-004

Norma2=norm(Err2,inf)Norma2 = 9.2370e-005

Limite

teorico

2

2 2 2

2 2 3

2 2 4

tmax t0 12.57*10

8 12

5 0.3142 2.57*10 2.5*10

10 0.1571 2.57*10 6.3*10

h

n

n

e h h

N h e h

N h e h

Limite

numerico

17

c) Istruzioni per il metodo shooting: RK4+secanti

clear all; clct0=0;tmax=pi/2; alpha=1; beta=(pi/2)-1;N=[5 10]; Yvera=inline('cos(x)-sin(x)+x');% metodo shootings0=-1;s1=1; fs=strvcat('y(2)','-y(1)+t');m=[1 0; 1 0];c0=0;c1=-1;for n=N [T,Y]=Shooting(t0,tmax,alpha,beta,n,m,c0,c1,fs,s0,s1); Err=abs(Yvera(T)-Y(:,1)); tab=[T Y(:,1) Err]; fprintf(' Tempo soluzione \t Err\n') fprintf('%7.4f %15.8f %15.3e \n',tab')end

18

Tempo soluzione Err 0.0000 1.00000000 0.000e+000 0.3142 0.95617735 2.143e-005 0.6283 0.84952440 2.587e-005 0.9425 0.72122733 1.872e-005 1.2566 0.61459013 7.408e-006 1.5708 0.57079633 3.331e-016

Metodo shooting RK4+secanti: risultati per N = 5

19

Metodo shooting RK4+secanti: risultati per N = 10

Tempo soluzione Err2 0.0000 1.00000000 0.000e+000 0.1571 0.98833293 5.811e-007 0.3142 0.95619788 9.052e-007 0.4712 0.90825392 1.007e-006 0.6283 0.84954934 9.329e-007 0.7854 0.78539743 7.354e-007 0.9425 0.72124558 4.737e-007 1.0996 0.66254119 2.096e-007 1.2566 0.61459753 4.633e-009 1.4137 0.58246290 8.281e-008 1.5708 0.57079633 1.110e-016

20

Metodo shooting: RK4+ interpolazione

clear all; clct0=0; tmax=pi/2; n=10;f=strvcat('y(2)','-y(1)+t');s=-0.25:0.25:2;ncase=length(s);b=zeros(1,ncase);l=[];for i=1:ncase y0=[1,s(i)]; [T,Y]=Rungekutta4(t0,tmax,n,y0,f); hold on; plot(T,Y(:,1)); grid on; title('Soluzioni ottenute con i vari valori di s'); xlabel('Tempo'); ylabel('Soluzione'); l=num2str(s(i)); legend(l); gtext(l); [m,nn]=size(Y); b(1,i)=Y(m,1);endp=polyfit(b,s,ncase-1);so=polyval(p,pi/2-1);fprintf('Valore ottimale di s: %20.15f \n',so);y0=[1,so];[T,Y]=Rungekutta4(t0,tmax,n,y0,f);fprintf('Valori trovati nei nodi \n');fprintf('%6.3f %20.12f \n',[T,Y(:,1)]');plot(T,Y(:,1),'r');lo=num2str(so);gtext(lo);legend(lo);

21

Shooting + Interpolazione

22

SOLUZIONE per s = so

23

RISULT. RK4 + Interpolazione

Nodi Valori approssimati 0.000 1.000000000000 0.157 0.988332927165 0.314 0.956197882110 0.471 0.908253915178 0.628 0.849549339920 0.785 0.785397428040 0.942 0.721245580280 1.100 0.662541194674 1.257 0.614597534883 1.414 0.582462901368 1.571 0.570796326795

Nodi Valori approssimati 0.000 1.000000000000 0.157 0.988332927165 0.314 0.956197882110 0.471 0.908253915178 0.628 0.849549339920 0.785 0.785397428040 0.942 0.721245580280 1.100 0.662541194674 1.257 0.614597534883 1.414 0.582462901368 1.571 0.570796326795

1 0.57079632679490!!!2

Valore di so: -8.939137842389755e-006 pertanto le soluzioni corrispondenti a s=0 e s= so, coincidono nel grafico .

24

File ”limiti.m ” - parte 1

clear all; clct0=0;tmax=pi/2;alpha=1;beta=(pi/2)-1;N=input('N =');% METODO linear_odep='0';q='-1';r='t';m=[1 0; 1 0];% METODO shooting (RK4+secanti)s0=-1;s1=1; m=[1 0; 1 0]; c0=0;c1=-1;fs=strvcat('y(2)','-y(1)+t');% Soluzione veraYvera=inline('cos(x)-sin(x)+x');

clear all; clct0=0;tmax=pi/2;alpha=1;beta=(pi/2)-1;N=input('N =');% METODO linear_odep='0';q='-1';r='t';m=[1 0; 1 0];% METODO shooting (RK4+secanti)s0=-1;s1=1; m=[1 0; 1 0]; c0=0;c1=-1;fs=strvcat('y(2)','-y(1)+t');% Soluzione veraYvera=inline('cos(x)-sin(x)+x');

25

File ”limiti.m ” - parte 2

% Implementazione dei metodi

[T,YL]=Linear_ODE(t0,tmax,alpha,beta,N,m,p,q,r);[T,Ys]=Shooting(t0,tmax,alpha,beta,N,m,c0,c1,fs,s0,s1);% Confronto con la soluzione veraErr_L=abs(Yvera(T)-YL);Err_s=abs(Yvera(T)-Ys(:,1));% Stampa dei risultati linear_odetab=[T YL Ys(:,1) Err_L Err_s ];fprintf(' Tempo soluzione_L soluzione_s\t Err_L Err_s\n')fprintf('%7.4f %15.8f %15.8f %15.2e %15.2e \n',tab')

26

Tabella riassuntivaN=10

Tempo soluzione_L soluzione_s Err_L Err_s

0.0000 1.00000000 1.00000000 0.00e+000 0.00e+000 0.1571 0.98840160 0.98833293 6.81e-005 5.81e-007 0.3142 0.95629116 0.95619788 9.24e-005 9.05e-007 0.4712 0.90833674 0.90825392 8.18e-005 1.01e-006 0.6283 0.84959737 0.84954934 4.71e-005 9.33e-007 0.7854 0.78539816 0.78539743 1.22e-015 7.35e-007 0.9425 0.72119896 0.72124558 4.71e-005 4.74e-007 1.0996 0.66245958 0.66254119 8.18e-005 2.10e-007 1.2566 0.61450517 0.61459753 9.24e-005 4.63e-009 1.4137 0.58239472 0.58246290 6.81e-005 8.28e-008 1.5708 0.57079633 0.57079633 1.11e-016 1.11e-016

Tempo soluzione_L soluzione_s Err_L Err_s

0.0000 1.00000000 1.00000000 0.00e+000 0.00e+000 0.1571 0.98840160 0.98833293 6.81e-005 5.81e-007 0.3142 0.95629116 0.95619788 9.24e-005 9.05e-007 0.4712 0.90833674 0.90825392 8.18e-005 1.01e-006 0.6283 0.84959737 0.84954934 4.71e-005 9.33e-007 0.7854 0.78539816 0.78539743 1.22e-015 7.35e-007 0.9425 0.72119896 0.72124558 4.71e-005 4.74e-007 1.0996 0.66245958 0.66254119 8.18e-005 2.10e-007 1.2566 0.61450517 0.61459753 9.24e-005 4.63e-009 1.4137 0.58239472 0.58246290 6.81e-005 8.28e-008 1.5708 0.57079633 0.57079633 1.11e-016 1.11e-016

27

Esercizio 2

Si consideri il seguente problema ai limiti:

a) Si verifichi che:

exp exp

exp 1 exp 1

t ty t t

è soluzione del problema (1);

b) Si applichi il metodo alle differenze finite con partizioni di N=10 e N=40 sottointervalli.

0,1(1)

0 0 1 0

y t y t t t

y y

28

Esercizio 2: quesiti c, dc) si costruisca una tabella che fornisca:

Intestazione: tempo soluzione errore Stampa dei risultati, ogni 5, con formato: 3 cifre decimali e virgola fissa per i tempi;8 cifre decimali e formato esponenziale per la soluzione;2 cifre decimali e formato esponenziale per l’errore assoluto.

d) Sulla stessa figura, si disegnino i grafici, per N=10, della soluzione vera e di quella approssimata.

29

a) Verifica della soluzione t

ee

eety

tt

11

1 11,

t te ey t

e e

11

ee

eety

tt

e sostituendo nel problema (1): tt

ee

ee

ee

ee tttt

1111

verifica le condiz. ai limiti

0011

00

ee

eey 011

11

11

ee

eey

Calcolando per le derivate :

OK!

OK!

30

b) metodo alle differenze finite

% file diff.mclear all; clct0=0;tmax=1;alpha=0;beta=0;N=input(' N= ');% function linear_odep='0';q='1';r='t';m=[1 0; 1 0];YV=inline('(exp(t)-exp(-t))/(exp(1)-exp(-1))-t');[T,Y]=Linear_ODE(t0,tmax,alpha,beta,N,m,p,q,r);

31

c) Istruzioni per la stampa dei risultati

% Calcolo dell’errore assolutoerr=abs(YV(T)-Y);% Stampa dei risultatitab=[T Y err ];disp('tempo soluzione errore')fprintf('%6.3f %17.8e %15.2e \n',tab(1:5:end,:)')

32

Metodo differenze: risultati N = 10

tempo soluzione errore

0.000 0.00000000e+000 0.00e+000 0.500 -5.65479225e-002 4.26e-005 1.000 0.00000000e+000 0.00e+000

tempo soluzione errore

0.000 0.00000000e+000 0.00e+000 0.500 -5.65479225e-002 4.26e-005 1.000 0.00000000e+000 0.00e+000

33

Metodo differenze: risultati N = 40

tempo soluzione errore

0.000 0.00000000e+000 0.00e+000 0.125 -1.83571711e-002 8.55e-007 0.250 -3.50459642e-002 1.64e-006 0.375 -4.83719003e-002 2.27e-006 0.500 -5.65878902e-002 2.67e-006 0.625 -5.78669016e-002 2.75e-006 0.750 -5.02733767e-002 2.41e-006 0.875 -3.17329508e-002 1.54e-006 1.000 0.00000000e+000 0.00e+000

tempo soluzione errore

0.000 0.00000000e+000 0.00e+000 0.125 -1.83571711e-002 8.55e-007 0.250 -3.50459642e-002 1.64e-006 0.375 -4.83719003e-002 2.27e-006 0.500 -5.65878902e-002 2.67e-006 0.625 -5.78669016e-002 2.75e-006 0.750 -5.02733767e-002 2.41e-006 0.875 -3.17329508e-002 1.54e-006 1.000 0.00000000e+000 0.00e+000

34

d) Istruzioni e grafico per N=10

plot(T,YV(T),T,Y,'*')xlabel('T')

35

Esercizio 3Si consideri il trasferimento di calore in una

barra isolata con dimensione dominante data

dalla lunghezza, con sezione trasversale di area

A, lunghezza L, conduttività k, coefficiente di

convezione , temperatura ambiente pari a Il

problema che si genera è: .T

0

0

ˆ0 0x L

d dTkA f x L

dx dx

dTT T kA T T q

dx

36

f è la sorgente distribuita di calore (spesso = 0);

la seconda condizione al limite tiene conto del trasferimento di calore convettivo attraverso l’estremo x=L e è il flusso di calore nel suddetto punto. Supponiamo:

a =kA=86Wcm/°C, L=10cm,

Il problema diventa:

q̂

20

ˆ10W/cm , 30 , 100 , 0, 0.T C T C f q

86 0 0 10

(1)

0 100 86 10 30 0x L

d dTx

dx dx

dTT T

dx

37

Caso soluzione regolare

Trasferimento di calore in una barra isolata

x

Temp. ambiente T∞

T= T0

38

Casi di possibile non regolarità della soluzione

L-L1L1

A2A1

(b)

T∞

(Ambiente)T=T0

L-L1L1

k2k1

(a) (a)barra di materiali diversi (compositi)

(b) Barra di sezioni diverse

39

Soluzione analitica del problema (1)

Si integrano due volte i due membri

dell’equazione differenziale; si ottiene:

Si impongono le condizioni:

Si ricava T(10) dalla soluzione generale, quindi

la soluzione vera è:

1286

cT x x c

2 1100 , 10 10 300 0c c T

700100

186T x x 700

100186

T x x

40

Soluzione numerica con shooting e differenze finite

Il problema viene risolto anche

numericamente utilizzando i metodi: shooting+secanti e differenze finite

mediante lo script memorizzato col nome

asta e riportato nella diapositiva seguente.

41

clc;clear allt0=0;tmax=10;alpha=100;beta=300;m=[1 0;10 86];c0=0;c1=-1;N=100;s0=1;s1=1.2;p='0';q='0';r='0';f=strvcat('y(2)','0');[T,Y,s]=Shooting(t0,tmax,alpha,beta,N,m,c0,c1,f,s0,s1);[T1,Y1]=Linear_ODE(t0,tmax,alpha,beta,N,m,p,q,r);xvera1=-700*T./186+100;xvera2=-700/186+T.*0; xvera=[xvera1,xvera2]; err=abs(xvera-Y); err1=abs(xvera1-Y1);tab=[T Y err];tabr=tab(1:10:end,:);ss='---------------------------------------------------------';disp(ss)fprintf('X Tapp. dTapp errT errdT\n')disp(ss)fprintf('%9.4f %15.7e %15.7e %8.2e %8.2e\n',tabr')fprintf('\n\n\n\n\n\n')tab1=[T1 Y1 err1];tabr1=tab1(1:10:end,:);disp(ss)fprintf('X Tapp. errT \n')disp(ss)fprintf('%9.4f %15.7e %8.2e \n',tabr1')

42

Risultati con shooting+secanti

X Tapp. dTapp errT errdT

0.0000 1.0000000e+002 -3.7634409e+000 0.00e+000 2.89e-014 1.0000 9.6236559e+001 -3.7634409e+000 1.42e-014 2.89e-014 2.0000 9.2473118e+001 -3.7634409e+000 2.84e-014 2.89e-014 3.0000 8.8709677e+001 -3.7634409e+000 4.26e-014 2.89e-014 4.0000 8.4946237e+001 -3.7634409e+000 5.68e-014 2.89e-014 5.0000 8.1182796e+001 -3.7634409e+000 5.68e-014 2.89e-014 6.0000 7.7419355e+001 -3.7634409e+000 7.11e-014 2.89e-014 7.0000 7.3655914e+001 -3.7634409e+000 8.53e-014 2.89e-014 8.0000 6.9892473e+001 -3.7634409e+000 9.95e-014 2.89e-014 9.0000 6.6129032e+001 -3.7634409e+000 1.14e-013 2.89e-014 10.0000 6.2365591e+001 -3.7634409e+000 1.21e-013 2.89e-014

X Tapp. dTapp errT errdT

0.0000 1.0000000e+002 -3.7634409e+000 0.00e+000 2.89e-014 1.0000 9.6236559e+001 -3.7634409e+000 1.42e-014 2.89e-014 2.0000 9.2473118e+001 -3.7634409e+000 2.84e-014 2.89e-014 3.0000 8.8709677e+001 -3.7634409e+000 4.26e-014 2.89e-014 4.0000 8.4946237e+001 -3.7634409e+000 5.68e-014 2.89e-014 5.0000 8.1182796e+001 -3.7634409e+000 5.68e-014 2.89e-014 6.0000 7.7419355e+001 -3.7634409e+000 7.11e-014 2.89e-014 7.0000 7.3655914e+001 -3.7634409e+000 8.53e-014 2.89e-014 8.0000 6.9892473e+001 -3.7634409e+000 9.95e-014 2.89e-014 9.0000 6.6129032e+001 -3.7634409e+000 1.14e-013 2.89e-014 10.0000 6.2365591e+001 -3.7634409e+000 1.21e-013 2.89e-014

43

Risultati con differenze finite

X Tapp. errT

0.0000 1.0000000e+002 0.00e+000 1.0000 9.6236559e+001 6.68e-013 2.0000 9.2473118e+001 1.31e-012 3.0000 8.8709677e+001 1.58e-012 4.0000 8.4946237e+001 1.59e-012 5.0000 8.1182796e+001 1.51e-012 6.0000 7.7419355e+001 1.46e-012 7.0000 7.3655914e+001 1.65e-012 8.0000 6.9892473e+001 1.75e-012 9.0000 6.6129032e+001 1.73e-012 10.0000 6.2365591e+001 1.57e-012

Commentare e motivare la differenza dei risultati ottenuti con i due metodi.

44

Si determini lo spostamento verticale di una fune lunga 1 u.m., soggetta ad un carico trasversale di intensità w per unità di lunghezza. T è la tensione lungo l’asse x e k è il coefficiente di elasticità della fune. Se u(x) è lo spostamento verticale, si ha il seguente problema :

4 - Compito d’esame (07/01/2003)

0,1

0 0 1 0

w xku x u x x

T Tu u

Si assume 1 sin 4 , 1 0.1w x x T k

45

1) Si dica motivando la risposta, se sono verificate le condizioni affinché esista unica, la soluzione del problema.

2. Si costruisca un file MATLAB: Cognome_studente_matricola.m, che una volta avviato:

a) faccia visualizzare una schermata con i dati personali ed una breve presentazione del problema;

Quesiti 1, 2a

46

Quesiti 2b, 2cb) Risolva il problema col metodo alle diff. finite con i passi

.40,20,10,1

jj

h j

c) costruisca una tabella che riporti l’intestazione:

nodi soluzione1 soluzione2 soluzione3;

su ogni riga, il valore di xi per i nodi coincidenti

nelle varie partizioni, il valore della soluzione

calcolata nei nodi, con i formati di stampa:

3 cifre decimali e formato virgola fissa per xi;

12 cifre decimali e formato virgola fissa per i valori della

soluzione nei 3 casi considerati.

47

1) Verifica dell’esistenza ed unicità della soluzione

0

0 1 0 1 0 1 0 1

0 0

, , 0,1 min ,max

0.1 0 0,1 , 0,

1 0 0, 1 0 0

2 0

u u

w xkf x u u u C u u

T Tk

f x fT

a a a a b b b b

a b

Poiché :

! la soluzione! OK Teorema A

N.B. Il problema dato è lineare si potrebbe utilizzare anche il Teorema dell’alternativa.

48

2a) file XXXX_XXX.mclcdisp('Cognome e nome studente: XXXX_XXX')disp('Numero di matricola: XXXX')disp(' Corso di Laurea: XXXX ')disp(' ')disp(' Questo programma consente di calcolare ')disp(' e visualizzare la soluzione ')disp(' del seguente problema ai limiti ')disp(' / ')disp(' \ -u''''(x)+ k/T u(x)= w/T ')disp(' / u(0)=0 u(1)=0')disp(' \ ')disp(' ')disp(' utilizzando il metodo delle differenze finite. ')

49

2b) Istruzioni clear allt0=0;tmax=1;alpha=0; beta=0;T=1;k=0.1;c=k/T;% METODO linear_odeq=num2str(c,'%25.16e');p='0';r=['(-1-sin(4*pi*t))/',num2str(T,' %25.16e')];m=[1 0; 1 0];j=[10 20 40];N=j;nj=length(N);tab=[];st= '%6.3f ';for i =1:nj [T,Y]=Linear_ODE(t0,tmax,alpha,beta,N(i),m,p,q,r); tab=[tab,Y(1:round(N(i)/N(1)):end)]; st=[st ' %16.12f '];endtab=[T(1:round(N(nj)/N(1)):end) tab];disp('nodi soluzione1 soluzione2 soluzione3')fprintf([st '\n‘],tab')

50

2c) Risultati

nodi soluzione1 soluzione2 soluzione3

0.000 0.000000000000 0.000000000000 0.000000000000 0.100 0.051468557578 0.050815060106 0.050663617021 0.200 0.083478018551 0.083077332036 0.082984534444 0.300 0.099693105020 0.100105126571 0.100200758273 0.400 0.111885737116 0.112549342408 0.112703312762 0.500 0.123700820113 0.123708471034 0.123710384002 0.600 0.125639603929 0.124990690864 0.124840394022 0.700 0.108193462187 0.107794302783 0.107701887017 0.800 0.074977661384 0.075388155824 0.075483405701 0.900 0.037714690765 0.038373711650 0.038526535761 1.000 0.000000000000 0.000000000000 0.000000000000

nodi soluzione1 soluzione2 soluzione3

0.000 0.000000000000 0.000000000000 0.000000000000 0.100 0.051468557578 0.050815060106 0.050663617021 0.200 0.083478018551 0.083077332036 0.082984534444 0.300 0.099693105020 0.100105126571 0.100200758273 0.400 0.111885737116 0.112549342408 0.112703312762 0.500 0.123700820113 0.123708471034 0.123710384002 0.600 0.125639603929 0.124990690864 0.124840394022 0.700 0.108193462187 0.107794302783 0.107701887017 0.800 0.074977661384 0.075388155824 0.075483405701 0.900 0.037714690765 0.038373711650 0.038526535761 1.000 0.000000000000 0.000000000000 0.000000000000

51

Grafico soluzione approssimata

plot(T,Y), grid,title('Soluzione approssimata')

N=40

52

Metodo alle differenze: caso non lineare

0 max

0

max

, , ,y t f t y t y t t t t

y t =

y t =

21 12

2i i ii

y y yy x O h

h

21 1'2

i ii

y yy x O h

h

La funzione f è non lineare. Utilizzando:

53

Metodo alle differenze: caso non lineare

0

1 1 1 12

1

2, , 1, ...,

2i i i i i

i i

n

y =

y y y y yf t y i n

hh

y =

h h hA y f y b

2 2

2 1 0 0

01 2 1 01 1

0 1

01

0 1 2

hA bh h

Sistema di equazioni non lineari

con

54

Esercizio 5 (esame 09-01-2006)

Si consideri il seguente problema ai limiti:

.0)1(,0)0(

)1,0(,0''

yy

xey y

.0)1(,0)0(

)1,0(,0''

yy

xey y

1) Si provi che esso ha un’unica soluzione che è non negativa in [0, 1].

2) Si costruisca un file Matlab: Cognome_Nome_studente.m che, una volta avviato:

55

Quesiti 2a), 2b), 2c) e 2d)

a) faccia visualizzare una schermata con i dati personali ed una breve presentazione del problema; b) permetta di dare in input il numero di sottointervalli spaziali N1=20, N2=40;c) risolva il problema dato nei due casi con il metodo Shooting+secanti e con le differenze finite;d) Visualizzi una tabella che riporti l’intestazione: x sol1s sol2s sol1dif sol2dif con le quantità x, sol1s, sol2s, sol1dif, sol2dif rappresentanti, rispettivamente i nodi coincidenti nelle due partizioni, la soluzione numerica in ognuno dei due casi, calcolata con i due metodi, con i seguenti formati di stampa:

56

Quesiti 2e) e 3) 3 cifre decimali e formato virgola fissa per i nodi, 8 cifre decimali e formato virgola mobile per la soluzione nei due casi;e) costruisca una figura in cui si riporti, mediante subplot con due finestre su una riga, nella prima finestra la soluzione approssimata calcolata con i due metodi nel caso N1=20, nella seconda la soluzione approssimata con i due metodi nel caso N2=40.3) Si commentino i risultati ottenuti specificando quante cifre della soluzione numerica si possono ritenere corrette utilizzando i due metodi.

57

Quesito1:Esistenza ed unicità della soluzione

Si consideri il problema associato

'' , , (0,1)

(0) 0, (1) 0.

y f x y x

y y

'' , , (0,1)

(0) 0, (1) 0.

y f x y x

y y

con se 0,

2 se 0

y

y

e yf x y

e y

58

Proprietà del problema associatoPer il nuovo problema risulta:

se 0, 0.

se 0

y

y yy

e yf f

e y

Sono pertanto soddisfatte le ipotesi del teorema di esistenza ed unicità della soluzione sulla funzione f .Le ipotesi sui coefficienti delle condizioni ai limiti sono banalmente soddisfatte.Il problema associato ammette quindi soluzione unica cheè sempre positiva per le condizioni ai limiti. Si deduce allora che anche il problema assegnato ha soluzione unica e sempre positiva.

59

Istruzioni relative ai quesiti 2b) e 2c)clear all;clc t0=0;tmax=1;alfa=0;beta=0; m=[1 0; 1 0];N=[20 40];% richiesti in inputfun='-exp(-y(i ))';dfi_1='0'; dfi='exp(-y(i ))'; dfip1='0';c0=0; c1=-1;f1='y(2)'; f2='-exp(-y(1))';f=strvcat(f1,f2);s0=-0.5;s1=0.5; % Metoditab=[];st= '%8.3f ';for i=1:length(N) y0=zeros(N(i)+1,1); [T,Y_df]=Nonlin_ODE(t0,tmax,alfa,beta,N(i),fun,dfi_1, dfi,dfip1,y0); [T,Y_s,s]=Shooting(t0,tmax,alfa,beta,N(i),m,c0,c1,f, s0,s1); tab=[tab Y_s(1:round(N(i)/N(1)):end,1) Y_df(1:round(N(i)/N(1)):end)]; st=[st, ' %16.8e %16.8e '];end

60

Istruzioni relative al quesito 2d)

% Formattazione tabella st=[st '\n'];tab=[T(1:round(N(end)/N(1)):end) tab]; tab(:,[3,4])=tab(:,[4,3]); str='________________________________________________';disp(' TABELLA RIASSUNTIVA: SHOOTING-DIFFERENZE')disp(str)fprintf(' x \t\t Y1_s \t\t Y2_s \t\t Y1_df \t Y2_df \n')disp(str)fprintf(st,tab')

61

Istruzioni relative al quesito 2e)

% Grafici (istruzioni all’interno del ciclo for)if i==1 c1= 'r'; c2= 'g o';else c1= 'm-.'; c2= 'g ';endsubplot(1,length(N),i),plot(T,Y_s(:,1),c1,T,Y_df,c2);title(['M. shooting e differenze N=',num2str(N(i))])xlabel('x');ylabel('y');legend('Y_s','Y_df')

62

Risultati TABELLA RIASSUNTIVA: SHOOTING-DIFFERENZE________________________________________________________________________________ x Y1_s Y2_s Y1_df Y2_df ________________________________________________________________________________0.000 0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000 0.050 2.19409754e-002 2.19409770e-002 2.19367745e-002 2.19399259e-002 0.100 4.14356200e-002 4.14356230e-002 4.14277939e-002 4.14336649e-002 0.150 5.85311907e-002 5.85311949e-002 5.85202667e-002 5.85284616e-002 0.200 7.32683762e-002 7.32683814e-002 7.32548417e-002 7.32649950e-002 0.250 8.56816979e-002 8.56817038e-002 8.56660067e-002 8.56777779e-002 0.300 9.57998406e-002 9.57998472e-002 9.57824200e-002 9.57954886e-002 0.350 1.03645923e-001 1.03645931e-001 1.03627179e-001 1.03641241e-001 0.400 1.09237713e-001 1.09237721e-001 1.09218034e-001 1.09232797e-001 0.450 1.12587788e-001 1.12587795e-001 1.12567553e-001 1.12582733e-001 0.500 1.13703648e-001 1.13703656e-001 1.13683229e-001 1.13698547e-001 0.550 1.12587788e-001 1.12587795e-001 1.12567553e-001 1.12582733e-001 0.600 1.09237713e-001 1.09237721e-001 1.09218034e-001 1.09232797e-001 0.650 1.03645923e-001 1.03645931e-001 1.03627179e-001 1.03641241e-001 0.700 9.57998407e-002 9.57998472e-002 9.57824200e-002 9.57954886e-002 0.750 8.56816979e-002 8.56817038e-002 8.56660067e-002 8.56777779e-002 0.800 7.32683763e-002 7.32683814e-002 7.32548417e-002 7.32649950e-002 0.850 5.85311907e-002 5.85311949e-002 5.85202667e-002 5.85284616e-002 0.900 4.14356201e-002 4.14356231e-002 4.14277939e-002 4.14336649e-002 0.950 2.19409754e-002 2.19409770e-002 2.19367745e-002 2.19399259e-002 1.000 7.97972799e-017 1.73472348e-018 0.00000000e+000 0.00000000e+000

63

Grafici