LUISS Corso di Metodi matematici per economia e

LUISS Corso di Metodi matematici per economia e finanza. Prof. F. Gozzi a. a. 2009/10 Integrazione numerica di equazioni differenziali - Prof. ssa G. Rotundo 1

Elenco degli argomenti Formule di integrazione numerica di EDO basati sulla formula di Taylor A. Metodo di Eulero 1. Analisi della qualità del metodo 2. Stima dell’errore globale di approssimazione a. Illustrazione tramite un esempio b. Stima (con dimostrazione) B. Alri metodi di ordine superiore: 1. Metodo di Eulero modificato 2. Metodo di Heun 3. Metodo di Runge-Kutta Formule per il calcolo di integrali Metodo dei rettangoli Metodo dei trapezi Metodo di Simpson 1/3 Corrispondenza tra metodi di calcolo numerico di EDO ed integrali Corrispondenza tra il metodo di Eulero (modificato) ed il metodo dei rettangoli Corrispondenza tra il metodo di Heun ed il metodo dei trapezi Modello preda-predatore Descrizione al modello Formalizzazione tramite un sistema di equazioni differenziali Linearizzazione Integrale primo, grafici di curve di livello e traiettorie. 2

Testi di riferimento MATLAB – manuale di riferimento R. A. Adams, Calcolo Differenziale 2, Casa Editrice Ambrosiana, ISBN 978 – 88 – 408 – 1390 -5 I. Capuzzo Dolcetta, M. Falcone, L’analisi al calcolatore, Zanichelli, ISBN 88 – 03904 - 8 J. Stoer, Introduzione all’analisi numerica, Zanichelli ed. , 1974 3

16 nov 09 Il problema di Cauchy Sotto opportune ipotesi è possibile dimostrare il teorema di esistenza ed unicità delle soluzioni. N. B. : il teorema di esistenza ed unicità fornisce un risultato soltanto in merito ad esistenza ed unicità, ma non illustra alcun metodo per trovare la soluzione. Domanda: come trovare le soluzioni? Anche quando la soluzione esiste ed è unica può essere molto complesso, se non impossibile, determinare la sua espressione analitica. Diventa quindi estremamente importante conoscere 4 alcuni metodi numerici per approssimare la soluzione.

Metodi numerici di integrazione di EDO Basati sulla formula di Taylor A. Metodo di Eulero 1. Analisi della qualità del metodo 2. Stima dell’errore globale di approssimazione a. Illustrazione tramite un esempio b. Stima (con dimostrazione) B. Alri metodi di ordine superiore: 1. Metodo di Heun 2. Metodo di Eulero modificato 3. Metodo di Runge-Kutta 5

Il metodo di Eulero Leonhard Euler (1707 -1783) 6

Il metodo di Eulero E’ il metodo più semplice. Per illustrarlo si parte dalla definizione di derivata: Per h abbastanza piccolo il rapporto incrementale è quasi uguale alla derivata: Sostituendo si ha che da cui 7

La relazione suggerisce di definire una soluzione approssimata a partire dalla successione definita per ricorrenza da Idea: fisso h e valuto l’espressione sull’insieme di punti La quantità h viene chiamata lunghezza del passo oppure passo dello schema di approssimazione 8

In corrispondenza a tk si può definire una successione di xk : 9

La funzione approssimante, lineare a tratti, può essere descritta, t: 0 t T, dalla espressione: x 2 x 1 x 0 t 1 t 2 Interpretazione geometrica del metodo di Eulero: il punto successivo è ottenuto dal precedente spostandosi con pendenza fissata per un intervallo di ampiezza h. La pendenza è data dalla funzione f calcolata nel primo estremo di ciascun 10 intervallo.

Esempio 1 Considero la seguente funzione: Si verifica che soddisfa una equazione differenziale, infatti: Quindi è soluzione della equazione differenziale: con 11

Grafico della soluzione x(t)=t 2 e della sua approssimazione tramite il metodo di Eulero Considero la condizione iniziale t 0=0. 5 e x(0. 5)=0. 25; passo h=0. 2: la generica iterazione col metodo di Eulero è: Problema: la successione approssimante (spezzata in nero con i punti evidenziati) si allontana ben presto dalla soluzione. E’ necessario dare una stima della bontà della approssimazione: Ricordando che: x’(t)=f(t, x(t)) 12

Confronto con lo sviluppo di Taylor arrestato al I ordine Osservo che le due condizioni x’(t)=f(t, x(t)) sono descritte per t=0 dallo sviluppo di Taylor: da cui ricavo una stima dell’errore dell’approssimazione di x(h) tramite x(0)+hx’(0) : (*) Ammesso che esista 13

Esempio 2 Considero la seguente funzione: Si verifica che soddisfa un problema di Cauchy: 14

Costruzione della approssimazione numerica della soluzione di Osservo che la variabile t non appare da sola, ma soltanto tramite x, quindi l’approssimazione di Eulero diventa: Osservo che la funzione f(. ) è la funzione identica (restituisce l’argomento, immutato). Quindi l’approssimazione di Eulero diventa: Scrivo il primo punto che si ottiene tramite l’approssimazione con il metodo di Eulero: Scrivo i punti successivi che si ottengono applicando il metodo: 15

Quindi la funzione soluzione del problema è: e le iterate approssimanti sono: k t x(t) xk 0 0. 00 1 0. 40 1. 49 1. 40 2 0. 80 2. 23 1. 96 3 1. 20 3. 32 2. 74 4 1. 60 4. 95 3. 84 Osservo che anche in questo caso la distanza (barre verticali in nero) della soluzione da quella 16 approssimata aumenta al crescere di k.

Programma matlab eulero. m h=0. 4; T=h*4; t=[0: h: T]; %t assume 5 valori: 0, 0. 3, 0. 6, 0. 9, 0. 12 %gli indici slittano di una unità per via della memorizzazione in MATLAB for i=1: 1: length(t) xk(i)=(1+h)^(i-1); end x=exp(t); plot(t, x, 'k', t, xk, 'r+-'); xlabel('t'); ylabel('x(t)'); for i=1: 1: length(t) fprintf(' %2 d | %4. 2 f n', i, i, t(i), xk(i)); end 17

Nota Da questo momento in poi consideriamo, per semplicità, il problema autonomo 18

Osservazione Il metodo di Eulero viene detto metodo ad un passo perché il calcolo del punto successivo dipende unicamente dal punto precedente. Un generico metodo ad un passo è definito da una successione per ricorrenza del tipo Il metodo di Eulero corrisponde ad una particolare scelta di : 19

Analisi della qualità del metodo L’errore di approssimazione è piccolo se h è piccolo? Voglio determinare 1. se l’errore di approssimazione è piccolo se h è piccolo (procedimento: (a) definisco una misura per l’errore e (b) dimostro che il limite è 0 per h 0) 2. l’errore globale di approssimazione (devo esaminare cosa succede se itero l’approssimazione N volte ed N ∞) 20

1. (a)Misura della qualità di un metodo ad un passo Definisco l’errore locale di discretizzazione in x(0) : con 21

1. (b) Dimostro che Svolgimento: utilizzo (*) 22

Osservazioni Poiché l’errore locale di discretizzazione è dello stesso ordine di h, per h 0, si dice che il metodo di Eulero è un metodo del primo ordine. La dimostrazione della convergenza non garantisce che la soluzione approssimata sia buona, vediamo un ulteriore esempio nel paragrafo che segue. 23

2. Stima dell’errore globale di discretizzazione (a)Illustrazione del problema tramite un esempio: dipendenza della qualità della approssimazione dal numero di iterazioni. (b) Stima dell’errore globale di discretizzazione: dipendenza dall’ampiezza del passo. 24

2. (a) Illustrazione del problema tramite un esempio Considero nuovamente il problema di Cauchy La funzione soluzione del problema è: e le iterate della approssimazione di Eulero sono: Voglio stimare l’errore, cioè la differenza tra la soluzione esatta e quella approssimata. Pongo t=hk, così stimo l’errore solo nei punti della successione tk. Devo semplificare questa espressione in modo da capirne facilmente l’andamento quando il numero k di passi aumenta (all’infinito). 25

Promemoria: Teorema di Lagrange Sia f continua in [a, b], derivabile in (a, b). Allora esiste un punto c in (a, b) tale che La dimostrazione fa uso del teorema di Rolle, che a sua volta fa uso del teorema di Weierstrass e sono riportati sul libro di matematica generale. 26

Usando il teorema di Lagrange Con a=klog(1+h) e b=kh Perché >0, quindi e >1 Si può dimostrare anche che: 27

Tramite l’esempio ho dimostrato che la qualità dell’approssimazione peggiora (la distanza tra soluzione esatta e soluzione approssimata va all’infinito) se il numero delle iterazioni aumenta. Un altro aspetto del problema riguarda la dipendenza della qualità della approssimazione dal passo h. Nel prossimo paragrafo dimostramo, in generale, che se il passo h diminuisce la qualità della approssimazione migliora. 28

2. b Errore globale di discretizzazione Misura della qualità della approssimazione a tempi molto maggiori di quello in cui è data la condizione iniziale Definizione: l’errore globale della approssimazione è dato da: 29

23 nov 09 Stima di E(h) per il metodo di Eulero Ipotesi ulteriori: f ed f ‘ limitate in R Fissato T e scelto h=T/N, N intero positivo fissato, vogliamo dimostrare che Tesi Osservazione 1: Questa stima dimostra che si può approssimare la soluzione esatta con la precisione desiderata su un intervallo [0, T] di ampiezza qualsiasi. Osservazione 2: con le ipotesi fatte la soluzione esiste in [0, +∞) Osservazione 3: la stima di questa disuguaglianza permette di dimostrare la convergenza a zero di |E(h)| per h 0 Osservazione 4: se T è grande, h dovrà essere scelto molto piccolo se si vuole avere una buona approssimazione. 30

Dimostrazione Osservazione: se x è soluzione del PC-EDO, allora Caso particolare: Ponendo si ha dunque 31

Applicando il teorema di Lagrange e le ipotesi di limitatezza di f ed f ‘ si ha ora sommo e sottraggo la stessa quantità Vd pagina successiva |dk-1| 32

Osservazione: per t=(k-1)h+ e s=(k-1)h l’equazione diventa Ipotesi di limitatezza Pertanto Iterando questa relazione ‘all’indietro’ fino a k=0 e ricordando che d 0=x(0)-y 0=0 si trova: 33

Iterando questa relazione ‘all’indietro’ fino a k=0 e ricordando che d 0=x(0)-y 0=0 si trova: Scegliendo ora k=N e ricordando che h=T/N, osservando che (1+LT/N)N<e. LT q. e. d 34

Metodi di ordine superiore basati sullo sviluppo di Taylor Vogliamo ora studiare dei metodi di approssimazione che, a parità di passo h, diano luogo a errori locali e globali di approssimazione più piccoli rispetto a quelli ottenuti per il metodo di Eulero. In particolare, ci occuperemo di vari metodi ad un passo di ordine p, intendendo con questo che Ovvero che l’errore locale tende a zero con la stessa rapidità di hp per h 0 In tal caso si dice che l’errore locale è un infinitesimo di ordine p 35 e(t, h)=O(hp)

Utilizzando lo sviluppo di Taylor è possibile costruire facilmente dei metodi ad un passo di ordine comunque elevato, a patto che f sia derivabile un numero sufficiente di volte. Supponiamo dunque che f sia dotata di (p-1) derivate continue. Il suo sviluppo di Taylor in un intorno di t è Da cui si ricava 36

Poiché x è soluzione del PC-EDO E così via. Possiamo cioè esprimere tutte le derivate di y in funzione di f e delle sue derivate calcolate nel punto y(t). Di conseguenza per avere un metodo di ordine 3, ad esempio, basterà definire E la successione come 37

Calcolo l’errore locale di discretizzazione in x(0) : Ricordo anche l’espansione in serie di Taylor Da cui quindi Sostituendo l’espressione di (x 0, h) : 38

Osservazioni 1. Questi metodi possono essere usati per definire metodi di approssimazione ad un passo di ordine comunque elevato. 2. I metodi ottenuti in questo modo richiedono la conoscenza di una formula esplicita di ciascuna delle derivate di f che compaiono in . 3. Siccome le derivate possono non essere facili, in alcuni casi si possono sostituire con i corrispondenti rapporti incrementali, ma questo appesantisce il calcolo. 39

Altri metodi di ordine superiore Obbiettivo: Ottenere metodi che abbiano ordine di convergenza elevato senza complicare troppo la trattazione. Descriviamo la costruzione di alcuni metodi ad un passo di ordine 2 aventi queste caratteristiche. 40

Ipotesi: derivabile due volte Per ottenere un metodo di ordine due basterà allora scegliere in modo tale che per ogni x si abbia =0 41

Esempio Con A, B, C parametri non negativi da determinare. La condizione diventa Che è verificata se A, B, C soddisfano il sistema di equazioni 42

Possibili scelte per A, B, C Evidentemente le equazioni non determinano univocamente i parametri. A differenti soluzioni corrispondono differenti metodi ad un passo di ordine due. Metodo di Heun: A=1/2, B=1/2, C=1 Corrisponde ad utilizzare come pendenza la media delle pendenze all’inizio ed alla fine di ciascun intervallo della partizione su t. Metodo di Eulero modificato: A=0, B=1, C=1/2 Corrisponde ad utilizzare come pendenza la pendenza calcolata nel punto intermedio di ciascun intervallo della partizione. 43

Metodo di Runge-Kutta Si ottiene definendo Con Si può dimostrare che è un metodo ad un passo di ordine 4. 44

Osservazioni 1. Questi metodi hanno vantaggi rispetto a quelli basati sulla formula di Taylor. Per esempio si chiede solo che f abbia la derivata prima continua e sarà sufficiente calcolare f in un minor numero di punti. 2. Il risparmio del calcolo di f anche in un solo punto può portare ad un grande risparmio nel tempo di calcolo complessivo in applicazioni che necessitano di un elevato numero di iterazioni 3. Sotto l’ipotesi che f sia dotata di derivate continue e limitate si puo’ dimostrare che E(h)=O(h 2) per i metodi di Heun ed Eulero modificato E(h)=O(h 4) per il metodo di Runge-Kutta 45

Nota sui sistemi del primo ordine I metodi visti finora si estendono ad un sistema di N equazioni del primo ordine. Che in forma vettoriale diventa Lo spazio euclideo RN dove è assegnata la condizione iniziale e dove evolve la traiettoria x(t) prende il nome di spazio delle fasi 46

Metodi numerici per il calcolo approssimato di integrali 30 nov 09 • Metodo dei rettangoli • Metodo dei trapezi • Metodo di Simpson 47

Osservazione Nel caso in cui f(t, x(t)) dipende solo da t il problema di Cauchy corrispondente si riduce al calcolo di una particolare primitiva di f(t). 48

Obbiettivo: calcolare l’integrale di una funzione reale su [a, b] Strumenti noti da matematica generale: il teorema fondamentale del calcolo integrale ( teorema di Newton-Leibniz), che risolve il problema nel caso in cui f è continua ed è nota una primitiva di f, cioè una funzione derivabile tale che F’(t)=f(t) per ogni x: a<x<b. In tal caso vale la formula: 49

Perché cambiare metodo? L’applicabilità pratica di questo metodo analitico è limitata, infatti: 1. Esistono funzioni continue che si incontrano frequentemente nelle applicazioni dell’analisi matematica la cui primitiva non può essere espressa in termini di funzioni elementari ( per esempio sen(t)/t, 1/((1 -t 2)(1 -k 2 t 2))1/2, exp(-t 2); 2. Il calcolo può essere complicato ( p. es. nel caso di integrazione di funzioni razionali fratte, che si basano sulla conoscenza di radici di polinomi). 50

Obbiettivo: calcolare UN VALORE APPROSSIMATO dell’integrale Strumenti: metodi di integrazione numerica facilmente implementabili: Rettangoli, Trapezi, Simpson [. . . ] Idea: integrare, invece di f, una sua approssimazione g che si sappia integrare in maniera esatta (a meno, ovviamente, di errori di arrotondamento). 51

Punto fondamentale Calcolare a priori l’errore commesso con la specifica procedura adottata. Serve per garantire l’affidabilità del metodo e bisognerà pertanto stimare In ciascun metodo la costruzione di g è fatta definendo una partizione dell’intervallo [a, b]. Precisamente, si fissa un intero positivo N e si divide l’intervallo di integrazione [a, b[ in parti uguali con i punti Il numero h=(b-a)/N è il passo della discretizzazione 52

METODO DEI RETTANGOLI a b Devo calcolare l’area sottesa dal grafico della funzione. 53

METODO DEI RETTANGOLI Idea: costruisco la partizione di [a, b] Obbiettivo: approssimare l’area con rettangoli di base (tk+1 -tk) ed altezza uguale all’altezza della funzione nel punto intermedio di ciascun intervallo (tk+tk+1)/2 Preparo gli elementi che servono per questa approssimazione a=t 0 t 1 t 2 b=t. N 54

Considero i punti intermedi in ciascun intervallo a=t 0 t 1 t 2 b=t. N 55

Considero il valore della funzione in ciascun punto intermedio in ciascun intervallo a=t 0 t 1 t 2 b=t. N 56

Considero la funzione costante a tratti, avente, su ciascun intervallo, il valore costante uguale all’altezza della funzione calcolata nel punto intermedio dell’intervallo a=t 0 t 1 t 2 b=t. N 57

Considero la funzione costante a tratti, avente, su ciascun intervallo, il valore costante uguale all’altezza della funzione calcolata nel punto intermedio dell’intervallo Ovviamente g coincide con f in tutti i punti intermedi di ciascun intervallo. a=t 0 t 1 t 2 b=t. N 58

Approssimo l’area da calcolare con la somma delle aree (base x altezza) di questi rettangoli. Questo calcolo è immediato: • tutte le basi hanno ampiezza costante (b-a)/N • le altezze sono date dal valore della funzione nel punto intermedio dell’intervallo a=t 0 t 1 t 3 b=t. N 59

Osservazioni • Il metodo dei rettangoli è ispirato dalla definizione di integrale definito. • I è una particolare somma integrale. • Il numero I dipende dalla partizione scelta e quindi dal passo h (h=(b-a)/N). 60

Stima dell’errore 61

MEMO: da matematica generale Proprietà dell’integrale definito L’area si può calcolare dividendola in due parti: (a, c) e (c, b) y f(x) 0 x a c b 62

Applico la proprietà dell’integrale ripetutamente: Posso ripetere il procedimento per g e quindi 63

MEMO In generale Il valore assoluto di una somma è minore od uguale alla somma dei valori assoluti 64

Scrivo l’espressione di g Ho applicato la proprietà del valore assoluto MEMO: teorema di Lagrange: sotto opportune ipotesi f(b)-f(a)=f ‘ (c ) (b-a), a<c<b Applico il teorema considerando b=t e a= (tk+tk+1)/2, c= k , t< k < (tk+tk+1)/2 65

Osservazione t tk (tk+tk+1)/2 tk+1 Un punto t dista dal centro dell’intervallo meno di metà della lunghezza dell’intervallo. 66

Utilizzo questa proprietà per svolgere il passaggio 67

Conclusione: ho dimostrato che 68

METODO DEI TRAPEZI (o di Crank-Nicolson) Mentre nel metodo dei rettangoli si approssima f con g costante a tratti, nel metodo dei trapezi si approssima f con g costruita partendo dalle rette secanti in ciascun intervallo a = t 0 t 1 t 2 b=t. N Quindi la funzione g è lineare in ogni intervallo e coincide con f negli estremi di ciascun intervallo. Il metodo si chiama ‘dei trapezi’ perché approssimo l’area tramite trapezi (eventualmente degeneri) invece che tramite rettangoli. 69

L’approssimazione corrisponde al calcolo dell’area evidenziata in rosso. Osservazione: alcuni trapezi possono essere degeneri, per esempio composti da triangoli o rettangoli Obbiettivi: • calcolo dell’area • stima dell’errore 70

Osservazione: rispetto alla normale visualizzazione, i trapezi sono ruotati di 90 gradi: Base minore altezza Base maggiore MEMO: area di un trapezio (base minore + base maggiore)*altezza /2 71

MEMO: da matematica generale Equazione di una retta che passa per i punti y y 2 (x 1, y 1) e (x 2, y 2) y y 1 da cui x 1 x x 2 da cui 72

Calcolo della retta secante nel singolo intervallo Applico la formula f(tk+1) f(t) Considerando le particolari coordinate dei punti: f(tk) … tk t Osservo che la funzione g coincide con la funzione f negli estremi dei singoli intervalli ed è lineare in ogni intervallo. tk+1 … 73 t

Area approssimata = somma delle aree dei singoli trapezi Base minore+base maggiore altezza Perché gli intervalli hanno tutti la stessa ampiezza 74

Stima dell’errore Ipotesi in più: f derivabile due volte in (a, b) con derivata f’’ limitata. Calcolo l’errore nel singolo intervallo, poi sommo per ottenere l’errore totale. Inizio il calcolo partendo dall’espressione di g 75

Per il teorema di Lagrange Inoltre g ed f coincidono sui punti della partizione Sottraggo Da cui Ora cerco di semplificare questa quantità 76

Considero quindi Lavoro dapprima su Applico la formula di Lagrange ed ottengo Sostituisco ed ottengo: 77

Applico il teorema di Lagrange alla funzione derivata prima: Sostituisco ed ottengo che in ciascun intervallo 78

Errore totale: sommo sugli N intervalli Quindi la stima cercata è: 79

METODO DI SIMPSON I passi effettuati finora per il calcolo approssimato dell’integrale hanno portato da una prima approssimazione mediante una funzione costante a tratti ad una approssimazione tramite una funzione rettilinea a tratti, in cui l’errore va a zero più velocemente quando h 0. Il passo successivo riguarda l’approssimazione mediante una funzione quadratica. Passi ulteriori possibili riguardano l’approssimazione della funzione assegnata mediante polinomi di ordine maggiore. 80

Interpolata quadratrica Considero una funzione in [-1, 1] e la sua interpolata quadratica definita da: Dove le tre costanti sono determinate dalle tre condizioni che chiedono che le funzioni coincidano su quei tre punti. Osservo che l’unica soluzione del sistema lineare è: 81

Interpolata quadratrica Quindi Questo calcolo si generalizza facilmente ad un arbitrario intervallo [c, d] tramite un cambio di variabile Che trasforma l’intervallo [-1, 1] in [c, d]. L’integrale diventa quindi 82

Su un intervallo qualsiasi [c, d] g(x) f(x) con semiampiezza dell’intervallo h h=(m-c)=(d-m) c h m h d Il nome con cui il metodo è più conosciuto è “Simpson 1/3”, per distinguerlo da un altro metodo di Simpson in cui la costante è diversa (“Simpson 3/8”) 83

Su più intervalli consecutivi Applico la formula considerando i punti della partizione a gruppi di tre: iniziale ( c ), finale (d) ed intermedio (m) e sommando tutti i risultati t 0 t 1 t 2 Quindi la funzione g è lineare in ogni intervallo e coincide con f negli estremi di ciascun intervallo. b=t. N 84

Su più intervalli consecutivi Applico la formula considerando i punti della partizione a gruppi di tre: iniziale ( c ), finale (d) ed intermedio (m) e sommando tutti i risultati. (h=(b-a)/N) h t 0 h t 1 h t 2 h h h b=t. N 85

Osservazione Le formule ottenute per la quadratura sono casi particolari di formule di Newton-Cotes in quanto la partizione sull'asse delle ascisse è costituita da punti equidistanti. Le formule si dicono di tipo chiuso quando x 0=a, xn=b (per esempio metodo dei trapezi e metodo di Simpson); si dicono di tipo aperto quando x 0=a+q, b=b-q (per esempio nel metodo dei rettangoli). 86

Osservazione Nel caso in cui f(t, x(t)) dipende solo da t il problema di Cauchy corrispondente si riduce al calcolo di una particolare primitiva di f(t). In questo caso si può dimostrare che il metodo di Heun si riduce al metodo dei trapezi per il calcolo di Si può anche dimostrare che il metodo di Eulero modificato si riduce a quello dei rettangoli, il metodo di Runge-Kutta a quello di Simpson. Nella prossima sezione mostriamo la corrispondenza tra – il metodo di Eulero (modificato) ed il metodo dei rettangoli; – il metodo di Heun ed il metodo dei trapezi. 87

14 dic 09 Confronto tra metodi per la risoluzione di equazioni differenziali e metodi per il calcolo di integrali Mostriamo la corrispondenza tra • il metodo di Eulero (modificato) ed il metodo dei rettangoli; • il metodo di Heun ed il metodo dei trapezi. 88

Confronto tra metodi per la risoluzione di EDO e metodi per il calcolo di integrali Formule di integrazione numerica di EDO basati sulla formula di Taylor: A. Metodo di Eulero 1. Analisi della qualità del metodo 2. Stima dell’errore globale di approssimazione a. Illustrazione tramite un esempio b. Stima (con dimostrazione) B. Alri metodi di ordine superiore: 1. Metodo di Eulero modificato 2. Metodo di Heun 3. Metodo di Runge-Kutta Nel caso in cui f(t, x(t)) dipende solo da t il problema di Cauchy corrispondente si riduce al calcolo di una particolare primitiva di f(t). Formule per il calcolo di integrali Metodo dei rettangoli Metodo dei trapezi Metodo di Simpson 1/3 Mostriamo la corrispondenza tra il metodo di Eulero (modificato) ed il metodo dei rettangoli e tra il metodo di Heun ed il metodo dei trapezi. 89

Corrispondenza tra il metodo di Eulero (modificato) ed il metodo dei rettangoli Procedimento: scrivo le iterazioni di Eulero, raggruppo i termini in maniera diversa, mostro che, modificata, la formula è la stessa dei rettangoli. Ricordo che il metodo di Eulero individua la successione approssimante è definita per ricorrenza da 90

scrivo xk+1 raggruppando i termini in modo diverso (**) quindi 91

Quindi la trasformazione della approssimazione che si ottiene con il metodo di Eulero permette di scrivere Considero, per semplicità, il caso in cui x 0=0 : Ricordo l’approssimazione dell’integrale che si ottiene dal metodo dei rettangoli: Le due formule coincidono se , e questo di fatto significa che stiamo usando il metodo di Eulero modificato. 92

Infatti, la modifica che viene considerata nel metodo di Eulero modificato corrisponde a scrivere il generico passo della iterazione con ponendo A=0, B=1, C=1/2, quindi Reintroducendo esplicitamente la dipendenza dal tempo si ha: Poiché ora stiamo trattando il caso di funzioni che dipendono soltanto dal tempo l’espressione si semplifica e si ha: Questo significa che xk+1 puo’ essere scritto, analogamente a quanto svolto in (**), come Avendo posto, per semplicità, x 0=0 , fissando k=N e sostituendo si ottiene: Ovvero la stessa approssimazione del metodo dei rettangoli, che considera come altezza del rettangolo il valore della funzione nel punto intermedio di ciascun intervallo della partizione. Senza questa modifica, la formula di Eulero considera rettangoli con altezza data dal valore della 93 funzione nel primo estremo di ciascun intervallo.

Corrispondenza tra il metodo di Heun ed il metodo dei trapezi Nel metodo di Heun l’iterazione è definita da: con A=1/2, B=1/2, C=1 e quindi diventa Reintroducendo esplicitamente la dipendenza dal tempo si ha: Poiché ora stiamo trattando il caso di funzioni che dipendono soltanto dal tempo l’espressione si semplifica e si ha: ovvero il valor medio di f rispetto agli estremi del singolo intervallo della partizione. 94

Avendo posto, per semplicità, x 0=0 , fissando k=N e sostituendo si ottiene: Considerare come pendenza la media delle pendenze agli estremi di ciascun intervallo della partizione dovrebbe migliorare l’approssimazione rispetto al metodo di Eulero, in cui si considera la tangente in uno soltanto degli estremi del segmento. t t+h Ricordo il valore che si ottiene con il metodo dei trapezi: Base minore+base maggiore altezza 95 Cioè esattamente ciò che si ottiene con il metodo di Heun

21 dic 09 Esempi numerici (in excel) Mediante il metodo di Eulero determinare i valori approssimati della soluzione del problema ai valori iniziali, nell’intervallo [0, 1], usando (a) 5 passi di lunghezza h=0. 2 (b) 10 passi di lunghezza 0. 1 Calcolare l’errore ad ogni passo, dal momento che il problema coinvolge un’equazione (lineare) risolvibile in modo esatto la cui soluzione è y=x-1+exp(-x) 96

Soluzione (a) In questo caso abbiamo f(x, y)=x-y, x 0=0, y 0=1 e h=0. 2, per cui xn=n/5, yn+1=yn+0, 2*(xn-yn) e l’errore è y-yn, per n=0, 1, 2, 3, 4, 5. I risultati del calcolo si ottengono tramite il foglio elettronico. (b) h=0. 1, n=10. File excel: esercizio. xls Osservazione Si noti che l’errore alla fine del primo passo è circa un quarto dell’errore alla fine del primo passo trovato nella partte (a), ma l’errore finale in x=1 è solo la metà di quello della parte (a). Questo comportamento è caratteristico del metodo di Eulero. Se diminuiamo la lunghezza del passo h occorrono più passi per raggiungere dal punto di partenza il valore particolare cui vogliamo arrivare. 97

Esempio 2 Usare il metodo di Heun con h=0. 2 per trovare i valori approssimati della soluzione del problema ai valori iniziali dell’esempio 1 in [0, 1]. Confrontare gli errori con quelli ottenuti tramite il metodo di Eulero. [excel] Osservazione: gli errori sono notevolmente inferiori ad 1/10 di quelli ottenuti nell’esempio 1(a). Ovviamente sono necessari piu’ calcoli per ogni passo, ma il numero richiesto di valutazioni di f(x, y) è solo doppio di quello richiesto nell’esempio 1(a). Come nel caso dell’integrazione numerica, se f è complicata, è proprio il numero di valutazioni della funzione a costituire gran parte del ‘costo’computazionale del calcolo della soluzione numerica. 98

Esercizio Ricavare i valori in x=0. 4, x=0. 8 e x=1 della soluzione del problema ai valori iniziali mediante i metodi appena descritti ed usando passi di lunghezza h=0. 2, h=0. 1 e h=0. 05 per ciascun metodo. Che cosa indicano i risultati riguardo ai valori della soluzione in questi punti? Confrontare i risultati con la soluzione esatta y=1/(1 -x). 99

Modello preda-predatore (Lotka-Volterra) 22 dic 09 100

Le equazioni di Lotka-Volterra, note anche come equazioni preda-predatore, sono un sistema di equazioni differenziali non lineari del primo ordine. Tali equazioni forniscono un modello matematico in grado di descrivere la dinamica di un ecosistema in cui interagiscono soltanto due specie animali: una delle due come predatore, l'altra come la sua preda. Questa modellizzazione matematica è stata proposta indipendentemente da Alfred J. Lotka nel 1925 e Vito Volterra nel 1926. 101

Supponiamo che nello stesso ambiente convivano due sole specie di animali. Quelli della prima specie sono erbivori, mentre quelli della seconda, i predatori, sono carnivori e si nutrono degli animali appartenenti alla prima specie, le prede. Il nostro interesse è essenzialmente rivolto allo studio dei legami tra le due specie: come varia nel tempo il numero degli animali appartenenti alle due specie? Possono le due specie convivere nello stesso habitati oppure una di esse è destinata ad estinguersi? Che influenza hanno le caratteristiche biologiche delle singole specie nell’equilibrio finale dell’ecosistema? 102

Se centriamo la nostra analisi su questi problemi possiamo, almeno in prima approssimazione, trascurare alcuni aspetti biologici che sono certamente connessi allo studio della evoluzione delle due specie. Supponiamo, ad esempio, che il cibo della prima specie sia disponibile sul territorio senza limitazioni e ignoriamo ogni differenziazione di età e di sesso tra gli animali. E’ evidente che queste semplicifazioni portano a considerare un ecosistema astratto. Nella costruzione del modello viene quindi isolato l’aspetto che piu’ interessa, in questo caso l’interazione preda-predatore. 103

La validità del modello matematico come modello di previsione dipende pertanto da quanto i fattori trascurati siano ininfluenti rispetto allo studio delle interazioni predatore. Supponiamo che sia le prede che i predatori si riproducano con tassi di natalità costanti e si estinguano con tassi di mortalità costanti. 104

Quindi, indicati con x(t) ed y(t) rispettivamente il numero delle prede e dei predatori presenti sul territorio al tempo t, in assenza di predazione l’evoluzione nel tempo delle due specie è descritta dalle due equazioni differenziali: x’(t)= a x(t) y’(t)=-b y(t) con a>0, b>0. 105

Si tratta di due equazioni, ciascuna a variabili separabili, non collegate fra di loro. La soluzione è x’(t)= x(0)*exp(a t) y’(t)= y(0)* exp(-b t) Il che vuol dire che il numero delle prede tende a crescere indefinitamente, coerentemente con l’ipotesi che il loro cibo è disponibile senza limitazioni, mentre il numero dei predatori decresce verso zero perche’ non si è tenuto conto della attività di predazione, che è l’unica forma di sostentamento della seconda specie. 106

Per descrivere l’attività di predazione supponiamo che il numero delle prede diminuisca di un fattore proporzionale a x(t)*y(t) Questo numero si interpreta come numero di ‘incontri’ tra le due specie. Supponiamo che il numero dei predatori cresca in maniera analoga. 107

Si ottiene il seguente sistema: La determinazione esplicita della soluzione di questo sistema non è facile, a causa del legame dato dal termine non lineare x(t)y(t). 108

E’ possibile effettuare un’analisi qualitativa che fornisca in maniera semplice qualche informazione sull’evoluzione delle due specie. Il primo passo consiste nel determinare i punti di equilibrio, cioè quelli in cui la derivata prima si annulla. Omettiamo, per semplicità di scrittura, la dipendenza da t e risolviamo il sistema x’=0, y’=0. 109

Partiamo dal sistema: x’(t)= a x(t) – c x(t)*y(t) y’(t)= -b y(t) + d x(t)*y(t) Imporre x’=0 ed y’=0 vuol dire chiedere a x – c x*y=0 -b y + d x*y=0 Metto in evidenza: x( a – c y) = 0 y(-b + d x) = 0 Verifico che E 1=(0, 0) è soluzione e che E 2=(b/d, a/c) è soluzione. Ho quindi due punti di equilibrio. 110

Questo significa che se il numero iniziale delle prede e dei predatori è dato dalle coordinate di E 1 (oppure dalle coordinate di E 2), le due popolazioni (preda e predatore) conservano nel tempo esattamente lo stesso numero di membri. Naturalmente questa situazione è assai improbabile e ci interessa piuttosto sapere se, partendo da una arbitraria situazione iniziale (x(0), y(0)), l’evoluzione converga o meno ad uno degli equilibri. Analizziamo ora il comportamento qualitativo del sistema per una generica scelta dei dati iniziali x(0), y(0) strettamente positivi. 111

La situazione è descritta in figura, che rappresenta il segno delle componenti nel campo vettoriale f 1(x, y)=ax-cxy ; f 2(x, y)=-by+dxy nel piano delle fasi. Se si parte da un punto sull’asse delle x le prede tenderanno a moltiplicarsi perché non ci sono predatori sul territorio. Se invece la condizione iniziale è sull’asse delle ordinate, allora i predatori tenderanno ad estinguersi in assenza di prede e quindi la soluzione del sistema convergerà all’equilibrio E 1 112

E’ interessante osservare che non appena x o y sono uguali a zero anche le corrispondenti derivate x’ e y’ si annullano; : se il numero dei membri di una delle popolazioni è nullo ad un certo tempo, rimarrà nullo per tutti i tempi successivi. Partendo invece da un punto con componenti strettamente positive, il segno delle derivate x’ ed y’ varierà a seconda che ci si trovi nelle regioni I, III, o IV. 113

Quali tecniche conosciamo per studiare un sistema di equazioni non lineare? 1. Si linearizza e si studia il sistema linearizzato. 2. Tecniche alternative. 114

1. Dal sistema linearizzato La stabilità dei punti di equilibrio può essere determinata linearizzando il sistema utilizzando le derivate parziali. La matrice jacobiana del modello preda-predatore è 115

La matrice jacobiana calcolata nel primo punto di equilibrio è: Gli autovalori sono uguali, rispettivamente, ad a e –b. Siccome a>0 e b>0, si ha che gli autovalori i segni dei due autovalori sono sempre diversi. Dunque il punto di equilibrio nell'origine è un punto di sella. La stabilità di questo punto è importante: se fosse stabile, valori di popolazione diversi da zero potrebbero essere attratti da esso, e perciò la dinamica del sistema porterebbe all'estinzione di entrambe le specie per molti valori iniziali delle popolazioni. Visto però che il punto è di sella, l'equilibrio è instabile, e l'estinzione di entrambe le specie è quindi difficile (in effetti, può accadere solo se le prede vengono estinte completamente in modo artificiale, provocando la morte dei predatori a causa della mancanza di cibo. Se invece sono i predatori ad essere estinti, la popolazione delle prede cresce senza limite in questo 116 semplice modello).

La matrice jacobiana calcolata nel primo punto di equilibrio è: Gli autovalori si ottengono annullando il determinante: Dato che gli unici valori che annullano il determinante sono complessi, il punto è di equilibrio stabile. La parte reale è uguale a zero in entrambi i casi, quindi il punto di equilibrio è un centro. Ciò significa che il livello dei predatori e delle prede 117 è ciclico, ed oscilla attorno a questo punto fisso.

2. Tecniche alternative 118

119

120

y Proiezione sul piano (x, y) delle curve di livello H(x, y) 121 x

122

Per analizzare il comportamento del sistema fuori dai punti di equilibrio si può dimostrare che esso ammette la costante del moto H(x, y). Dunque tutte le traiettorie del sistema nello spazio x y giacciono sulle curve di livello della funzione H. La funzione ha un minimo corrispondente al punto di equilibrio stabile ed è convessa sul primo quadrante, da questo si può dedurre che in una situazione generica con due popolazioni iniziali x e y il sistema ha un comportamento oscillante che torna periodicamente nello stato iniziale, con oscillazioni anche molto grandi. Le soluzioni non possiedono una espressione semplice in termini delle funzioni trigonometriche. Se si approssimano le soluzioni che si trovano in un intorno del punto di equilibrio stabile linearizzando il sistema si ottiene un moto armonico semplice in cui la popolazione dei predatori precede quella delle prede con uno sfasamento uguale a π / 2. 123

124

Domande per l'esame - Elementi propedeutici di matematica generale o enunciato del teorema di Bolzano-Weierstrass o teorema di Rolle o teorema di Lagrange o formula di Taylor; dimostrazione della convergenza del resto di ordine uno o definizione di max, sup: quale relazione li lega? o equazione di una retta per due punti - Metodi per la risoluzione di equazioni differenziali o Illustra il metodo di Eulero o Mostra il modo in cui viene ricavato il metodo a partire dalla formula di Taylor o Calcola l'errore ad un passo o Calcola l'errore globale o Dimostra l'ordine di grandezza dell'errore rispetto al passo. o Scrivi il generico passo ricavato con il metodo di Eulero per la risoluzione dei problemi di Caucy: A. x'(t)=t^2; x(0)=0; B. x'(t)=x(t); x(0)=1; o Illustra il modo in cui si ricavano metodi di ordine superiore a partire dalla formula di Taylor o Illustra il modo in cui si ricava il metodo di Heun a partire dalla formula di Taylor 125 o Descrivi il generico passo della iterazione del metodo di Runge-Kutta ed il suo ordine di grandezza (solo enunciato).

- Metodi per la risoluzione di integrali o Illustra il metodo dei rettangoli o Calcola l'errore del metodo dei rettangoli o Dimostra l'ordine di grandezza dell'errore o Illustra il metodo dei trapezi o Calcola l'errore del metodo dei trapezi o Dimostra l'ordine di grandezza dell'errore o Illustra il metodo di Simpson - Legame tra metodi per la risoluzione di equazioni differenziali e metodi per la risoluzione di integrali o Illustra la problematica (in generale) o Mostra come il metodo di Eulero si riduce al metodo dei rettangoli nel caso di problemi di Cauchy in cui l'equazione differenziale è del tipo x'(t)=f(t) o Mostra come il metodo di Heun si riduce al metodo dei trapezi nel caso di problemi di Cauchy in cui l'equazione differenziale è del tipo x'(t)=f(t) - Modello preda-predatore o E' un sistema di equazioni lineari? o Quali tecniche sono disponibili per la sua risoluzione? o Cos'e' un integrale primo? o Quali particolarità hanno le traiettorie? o Come si riduce il sistema in assenza di attività predatoria? 126