Calcolo Numerico 2

Z

y(t ) = y(t ) + h f (t + a h, y(t + a h)) d a (2.3.3)

i i−1 i i i i i

0

Ora possiamo utilizzare una formula di quadratura per approssimare l’integrale sui nodi a con

i

dei pesi c scelti a seconda del metodo che si sta utilizzando.

i

Si ottiene quindi una approssimazione del valore y della funzione:

1

X

Y = Y + h c f (t + a h, Y ) (2.3.4)

n+1 n i i i i

i

dove Y sono valori approssimati di y(t + a h), ignoti a priori, da determinare anche essi attra-

i i i

verso formule di quadratura con pesi appositi.

Definizione 2.25. (Metodi di Runge-Kutta)

I metodi di Runge-Kutta sono una famiglia di metodi iterativi impliciti ed espliciti. In generale

un metodo di Runge-Kutta viene caratterizzato da tre parametri; un vettore ~a

, una matrice B

quadrata, anche chiamata matrice di Butcher, e un vettore ~c

.

a b b . . . b

1 1,1 1,2 1,p

a b b . . . b

2 2,1 2,2 2,p ~a B

.. .. ..

.. .. =

. .

. . . ~c

b b . . . b

a p,1 p,2 p,p

p c c ... c

1 2 p

Il metodo generico, chiamato a s stadi è nella forma p

s

X X

n n

Y = Y + h c K , dove K = f (t + a h, Y + h b K ) (2.3.5)

n+1 n i n i n i,j j

i i

i=1 j

I metodi di Runge Kutta possono essere sia espliciti, sia impliciti. Se la matrice B è triangolare

inferiore stretta, allora il metodo è esplicito. Se la matrice A è triangolare, allora il metodo si

dice semi-esplicito, altrimenti implicito. 18

Esempio 2.26. (Metodi gia visti)

Qui di seguito vi è un elenco dei metodi visti in precedenza reinterpretati come metodi di Rutte-

Kutta.

Eulero Esplicito (EE) (

Y = Y + h [K ]

0 0 n+1 n 1

1 K = f (t , Y )

1 n n

Eulero Implicito (EI) (

Y = Y + h [K ]

1 1 n+1 n 1

1 K = f (t + h, Y + h[K ]) = f (t , Y )

1 n n 1 n+1 n+1

Heun (H) h

 [K + K ]

Y = Y + 1 2

n+1 n 2



0 0 0 

 K = f (t , Y )

 1 n n

1 0

1 ∗

K = f (t + h, Y + K ) = f (t , Y )

2 n 1 n+1

1/2 1/2 n+1





 ∗

Y = Y + f (t , Y )

 n n n

n+1

Eulero Modificato - Punto Medio (EM)  Y = Y + h[K ]

0 0 0 n+1 n 2





 K = f (t , Y )

1/2 1/2 0 1 n n

0 1 1 1

 K = f t + h, Y + K = f t , Y

 1 1

2 1 n+ n+

 2 2 2 2

Trapezi (T)  h

Y = Y + [K + K ]

0 0 0 n+1 n 1 2

2





1 1/2 1/2 K = f (t , Y )

1 n n

1/2 1/2 h



K = f t + h, Y + [K + K ] = f (t , Y )

 2 n n 1 2 n+1 n+1

2

Consistenza

Supponiamo che K venga valutato nella soluzione esatta.

i 2

0 00 0 0

h

n 2

P P

− − − − −

hσ (h) = y y h c K = y + hy + y + o(h ) y h c (y + o(h)) = hy (1

n n+1 n i n n i

n n n n

i

i 2

2

P c ) + o(h )

i P

⇔

Si ha quindi che il metodo è consistente c = 1 e per essere di almeno primo ordine è

i

P ∀i

sufficiente avere a = b = 1 . . . p

i,j i

j 1

P

Si può dimostrare che l’ordine è invece almeno 2 se c b =

i i 2

Vale inoltre anche la seguente relazione tra l’ordinde di convergenza e il numero di stadi minimo

per i metodi di Runge-Kutta espliciti, riportato nella seguente tabella:

|1

ordine 2 3 4 5 6 7 8 . . .

|1

s 2 3 4 6 7 8 9 . . .

min

In particolare vale il seguente teorema

Teorema 2.27. ≤

Un metodo Runge-Kutta a s stadi, esplicito, ha un ordine di convergenza p s se s < 5,

≥

altrimenti ha un ordine di convergenza p < s per s 5

19

2.3.1 Codici d’implementazione

La seguente function è adatta per utilizzare i metodi di Runge-Kutta espliciti, cioè quelli in cui

la matrice A è triangolare inferiore stretta. n

Alla riga 19, viene applicata la formula dei metodi di Runge-Kutta: Y = Y + hcK (riscritta

n+1 n

in forma matriciale), è però necessario creare allora prima una matrice K che contiene le valuta-

zioni della f ad ogni passo in maniera ricorsiva.

1 function [T , Y ]= rkesp (f ,T , y0 ,a ,B , c )

% [T , Y ]= RKesp (f ,T , y0 ,a ,B , c )

3 % Input : f = stringa funzione di due variabili (t , y )

% T = vettore riga dei tempi

5 % y0 = dato iniziale

% a ,B , c parti di una matrice di Butcher , B e ’ la matrice , c il

vettore la cui somma da 1

7 % Output : T = vettore colonna dei tempi

% Y = vettore colonna soluzione approssimata agli istanti contenuti

in T

9 h = T (2) -T (1) ;

11 Y = y0 ’;

d = length ( Y ) ;

13 ns = length ( a ) ;

K = zeros (d , ns ) ;

15 for n =1: length ( T ) -1

for r =1: ns

17 K (: , r ) = feval (f , T ( n ) + a ( r ) *h , Y (: , n ) + h * K * B (r ,:) ’) ;

end

19 Y (: , n +1) = Y (: , n ) + h * K *c ’;

end

21 Y =Y ’; T =T ’;

end

Questa function invece è in grado di trattare i metodi di Runge-Kutta semi-impliciti, il metodo

usato per risolvere l’equazione implicita è quello di Newton.

function [T ,Y , nit ]= rksimp (f ,T , y0 , dfy ,a ,B ,c , toll , nitmax )

2 % Runge Kutta semi - implicito

% [T , Y ]= RKesp (f ,T , y0 ,a ,B , c )

4 % Input : f = stringa funzione di due variabili (t , y )

% T = vettore riga dei tempi ;

6 % y0 = dato iniziale ( vettore riga , ogni componente inizializza un ’ eqz

del sist ) ;

% a ,B , c parti di una matrice di Butcher , B e ’ la matrice , c il

8 % vettore la cui somma da 1

% Output : T = vettore colonna dei tempi ;

10 % Y = vettore colonna soluzione approssimata agli istanti contenuti

in T .

if nargin ==7

12 toll =1 e -12;

nitmax =100;

14 elseif nargin ==8

nitmax =100;

16 end 20

dim = length ( y0 ) ;

18 h = T (2) -T (1) ;

Y = y0 ’;

20 ns = length ( a ) ;

K = zeros ( dim , ns ) ;

22 bd = diag ( B ) ;

B =B - diag ( bd ) ;

24 for n =1: length ( T ) -1

for r =1: ns

26 it =0;

err = toll +1;

28 u = K (: , r ) ;

tnr = T ( n ) + a ( r ) * h ;

30 aux = Y (: , n ) + h * K * B (r ,:) ’;

while it < nitmax & err > toll

32 J = eye ( dim ) -h * bd ( r ) * feval ( dfy , tnr , aux + h * bd ( r ) * u ) ;

tn = feval (f , tnr , aux + h * bd ( r ) * u ) -u ;

34 delta = J \ tn ;

err = norm ( delta , inf ) ;

36 it = it +1;

u = u + delta ;

38 end

K (: , r ) = u ;

40 nit (n , r ) = it ;

if it == nitmax & err > toll

42 disp ( ’ non converge nel numero massimo di iterazioni assegnato ’)

disp ([ ’ per K numero ’ , num2str ( r ) ])

44 disp ([ ’ per h = ’ , num2str ( h ) , ’ in t = ’ , num2str ( T ( n ) ) ])

return

46 end

end

48 Y (: , n +1) = Y (: , n ) + h * K *c ’;

end

50 Y =Y ’; T =T ’;

Questa function è un particolare metodo di Runge-Kutta, in cui la matrice di Butcher è quella

relativa per il metodo dei trapezi.

function [T , Y ]= rk4 (f ,T , y0 )

2 % [T , Y ]= RKesp (f ,T , y0 )

% Input : f = stringa funzione di due variabili (t , y ) che restituisce un

vettore colonna ;

4 % T = vettore riga dei tempi ;

% y0 = dato iniziale ( vettore riga , ogni componente inizializza un ’ eqz

del sist ) ;

6 % Output : T = vettore colonna dei tempi ;

% Y = vettore colonna soluzione approssimata agli istanti contenuti

in T .

8 a =[0 1/2 1/2 1] ’;

10 B =[0 0 0 0; 1/2 0 0 0; 0 1/2 0 0; 0 0 1 0];

c =[1/6 1/3 1/3 1/6];

12 h = T (2) -T (1) ;

14 Y = y0 ’; 21

d = length ( Y ) ;

16 ns = length ( a ) ;

K = zeros (d , ns ) ;

18 for n =1: length ( T ) -1

for r =1: ns

20 K (: , r ) = feval (f , T ( n ) + a ( r ) *h , Y (: , n ) + h * K * B (r ,:) ’) ;

end

22 Y (: , n +1) = Y (: , n ) + h * K *c ’;

end

24 Y =Y ’;

T =T ’;

26 end

Questa function disegna le regioni di stabilità per i metodi di Runge-Kutta.

1 function rkregstab (B ,c , x1 , x2 , y1 , y2 , punti )

% disegna gli insiemi di assoluta stabilita ’

3 % x1x2y1y2 sono gli estremi del rettangolo su cui disegniamo , punti e ’ il

numero di punti per disegnare l ’ insieme

x = linspace ( x1 , x2 , punti ) ;

5 y = linspace ( y1 , y2 , punti ) ;

[X , Y ]= meshgrid (x , y ) ; % crea la ’ tabella ’ su cui valutare la funzione

7 Z = X + i * Y ; % cosi ho i numeri complessi

for l =1: punti

9 for m =1: punti

R (l , m ) = feval ( ’ funstab ’ ,Z (l , m ) ,B , c ) ;

11 end

end

13 figure ()

contour (X ,Y , abs ( R ) ,[1 ,1] , ’k ’) ;

15 axis ([ x1 x2 y1 y2 ]) ;

grid on ;

17 figure

spy ( abs ( R ) <=1)

19 end

21 function f = funstab (q ,B , c )

% valuta la funzione che serve per calcolare la stabilita

23 ns = length ( c ) ;

U = ones ( ns ,1) ;

25 f =1+ q * c * inv ( eye ( ns ) -q * B ) * U ;

end 22

2.4 Metodi di Estrapolazione

Definizione 2.28. Metodi di estrapolazione

⊆

Sia I la soluzione esatta del problema (ad esempio un integrale, oppure un equazione diffe-

R

renziale), sia I la soluzione approssimata (ottenuta ad esempio con una formula di quadratura)

n α

|I − ≤

Supponiamo di essere in grado di controllare l’errore; che valga I| ch , e in modo par-

n

α α α

− ∈

ticolare che I I = ch + o(h ) e che o(h ) sia sviluppabile nel seguente modo fino a m N

n

fisstao α 2α mα mα+1

− · · ·

I I = c h + c h + + c h + o(h )

n 1 2 m

α 2 m m

Se si può fare tutto questo, ponendo H = h , si ottiene: I = I +c H +c H +. . . c H +o(H ),

k 1 2 m

che è un polinomio di grado m del quale a noi interessa il termina noto.

Utilizzando quindi l’interpolazione polinomiale nel seguente modo possiamo approssimare il ter-

mine noto I:

1. Si scelgono m + 1 valori distinti per H

2. Si costruisce con questi valori il polinomio interpolante

3. Una volta ottenuto il polinomio, si ha anche il suo termine noto

Definizione 2.29. Algoritmo di Neville

Si tratta di una algoritomo simile alle differenze divise di Newton, ma permette di calolcolare un

polinomio valutato in un punto. È quindi più adatto ai nostri scopi perchè di tutto il polinomio

ci interessa un solo punto.

L’algoritmo funziona nel seguente modo, supponendo che (x , y ) siano i nostri nodi d’interpola-

i i

6 6

zione con la condizione x = x per ogni i = j e ponendo

i j · − · −

P (x) (x x ) + P (x) (x x)

i,i+1,...,j−1 j i+1,...,j i

P = y , P (x) =

i i i,i+1,...,j −

x x

i j

Da cui si ha ∀k

P (x ) = f = i, i + 1, . . . , j

i,i+1...,j k k

−

y (1) (2) ... (l 1) l

x

x y P P ... P P

0 0 0 0,1 0,...,l−1 0,...,l