Corso di EQUAZIONI DIFFERENZIALI ORDINARIEPrima esercitazione di laboratorio MATLAB1

Clelia Marchionna(clelia.marchionna@polimi.it)

Simone Pezzuto(simone.pezzuto@mail.polimi.it)

Versione del: 18 ottobre 2012

1 Introduzione

Matlab e un software per il calcolo scientifico ampiamente utilizzato da scienziati applica-ti e ingegneri in tutto il mondo. Sviluppato dalla MathWorks (http://www.mathworks.com),e un software protetto da copyright e non disponibile gratuitamente. Esistono tuttavia pro-grammi OpenSource che ne imitano (con successo) il funzionamento, in particolare Octa-ve (http://www.octave.org), che utilizza comandi di fatto identici per nome e sintassi, eSciLab (http://www.scilab.org), che seppure direttamente ispirato a Matlab utilizza unasintassi leggermente diversa. Entrambi questi software riproducono solo le funzionalita basedi Matlab. In particolare, per poter utilizzare le funzioni grafiche di Octave e necessariodisporre anche del software grafico GnuPlot, che e anche esso disponibile senza restrizioni(http://www.gnuplot.info).

Lo scopo di questa dispensa e quello di fornire una prima introduzione agli strumenti presentiin Matlab per la soluzione di equazioni e sistemi di equazioni differenziali ordinarie. Si presup-pone che il lettore gia possegga una certa familiarita con l’uso di questo software; per una primaintroduzione si rimanda alle dispense del corso EAMAG e ai Tutorials disponibili in rete. Nell’am-bito di queste esercitazioni si utilizzeranno infine numerosi scripts didattici sviluppati da JohnPolking presso la Rice University (USA) e reperibili sul sito http://math.rice.edu/˜dfield.

In questa dispensa si utilizzeranno comandi e grafica relativi alla versione R2012a. Si noti,quindi, che nel caso di versioni differenti gli output e la sintassi potrebbero essere leggermentedifferenti.

2 Richiami sul Symbolic Math Toolbox di MATLAB

Matlab viene generalmente utilizzato come software per la soluzione numerica di problemimatematici, ed e uno dei migliori strumenti attualmente disponibili a livello commerciale. Sesi dispone del pacchetto software aggiuntivo Symbolic Math Toolbox, Matlab pero puo essereanche utilizzato per eseguire calcoli simbolici, sebbene le sue funzionalita simboliche siano ingenerale meno efficaci di quelle di altri software scientifici (vedi Mathematica o Maple).

Per lavorare in modalita simbolica si devono dichiarare dei simboli. Per esempio, se inMatlab scriviamo

>> a = 1/3

1Contributi: questa dispensa e un’edizione che utilizza il materiale delle esercitazioni degli anni passati.Ringraziamo LORIS BIONDO, LUCA BONAVENTURA e NAOMI VALENTINA SLAVAZZA per le versioniprecedenti.

Prima esercitazione di laboratorio MATLAB Richiami sul Symbolic Math Toolbox di MATLAB

il software creera una variabile di nome a contenente il valore in virgola mobile (floating–point)piu vicino ad 1/3 (ma non esattamente). Per dichiarare invece il simbolo 1/3, si puo scrivere:

>> a = sym( ’1/3 ’ )

sym e un comando specifico del toolbox simbolico, e crea una variabile di nome a contenente ilsimbolo 1/3, interpretato come tale, e non in virgola mobile. E chiaramente possibile usare ilcomando anche per espressioni piu complesse, come:

>> b = sym( ’1+ sqrt (5) ’ )

che corrisponde a 1 +√

5. Bisogna fare attenzione che il comando sym converte la stringa fornitain argomento in un simbolo, ma e anche in grado di trovare il simbolo corrispondente ad unaespressione numerica (notare l’assenza delle virgolette):

>> c = sym(1/3)

Ancora una volta, c contiene il simbolo 1/3. Pero se proviamo con:>> d = sym(1+sqrt ( 5 ) )

il toolbox prova a calcolare il numero razionale piu vicino, e se ci riesce tratta quello comesimbolo. Se invece non e in grado di far cio, memorizza il numero nella forma floating–point2.

Per definire stringhe piu complesse, la notazione e sempre la stessa:>> f = sym( ’xˆ2 + sin(x)’ )

In questo caso sym interpreta la stringa come la somma di due addendi, dove il primo e il quadratodell’oggetto x, e il secondo e la funzione sin applicata ad x. Notare invece che il comando:

>> f = sym( xˆ2 + sin ( x ) )

fornisce un errore, in quanto la variabile x non e stata definita. Pero si puo scrivere:>> x = 1 ;>> f = sym( xˆ2 + sin ( x ) )

e in questo caso f contiene il valore (eventualmente approssimato) di 1 + sin(1).Ora definiamo due simboli:

>> f = sym( ’xˆ2 ’ ) ;>> g = sym( ’xˆ3 ’ ) ;

e proviamo a scrivere:>> h = f ∗g

La variabile h contiene il simbolo x5, ossia x2 · x3. Quindi si possono eseguire molte operazionialgebriche ai simboli definiti, e automaticamente verra creato il simbolo contenente quell’espres-sione. Questo suggerisce un’idea: supponiamo di definire il simbolo x, e che per comodita vengamemorizzato nella variabile di nome x

>> x = sym( ’x’ ) ;

Ora proviamo nuovamente il comando>> f = xˆ2 + sin ( x )

e questa volta l’interpretazione e quella corretta.Attenzione: variabile e oggetto corrispondente sono due cose differenti! Come specificato,

chiamare la variabile contenente il simbolo x con x e una pura comodita, non una necessita.Infatti:

>> x = sym( ’topolino ’ ) ;>> f = xˆ2 + sin ( x )

2Versioni diverse di Matlab potrebbero comportarsi in modo diverso da questo punto di vista.

– 2 –

Prima esercitazione di laboratorio MATLAB Richiami sul Symbolic Math Toolbox di MATLAB

restituisce l’espressione sin(topolino) + topolino2.Se l’espressione contiene molte variabili, puo diventare tedioso scrivere ogni volta il comando

sym:>> x = sym( ’x’ ) ;>> y = sym( ’y’ ) ;>> z = sym( ’z’ ) ;>> w = sym( ’w’ ) ;

Per questo il toolbox simbolico mette a disposizione il comando syms:>> syms x y z w;

che e l’equivalente compatto di quanto scritto prima.Con i simboli si possono fare un certo numero di operazioni algebriche, come espandere

un prodotto, semplificare l’espressione, calcolare derivate ed integrali, trasformate di Laplace oFourier, o risolvere equazioni differenziali. E inoltre possibile usare molti comandi “numerici” diMatlab, come calcolare determinante, traccia, autovalori. In quest’ultimo caso e utile ricordareche e possibile definire anche vettori o matrici di oggetti simbolici, in modo del tutto naturale:

>> syms x ;>> A = [ x , 0 ; 1−x , x ] ;

Per esempio il comando det(A) fornira il simbolo x2, mentre eig(A) fornisce il vettore [ xx ].

Proviamo con un esempio piu elaborato: data una matrice generica, vogliamo calcolare gli au-tovalori usando la definizione. Definite le variabili simboliche, costruiamo la matrice e calcoliamoil polinomio caratteristico:

>> syms a b c d lambda ;>> A = [ a b ; c d ] ;>> pol inomio = det (A − lambda∗eye ( 2 ) ) ;

Ora risolviamo il polinomio rispetto a λ, con il comando:>> e i g s = s o l v e ( pol inomio , lambda ) ;

Il vettore eig contiene gli autovalori, come si puo verificare confrontando il risultato con:>> eig (A) ;

Esercizio 1

Utilizzando la guida di Matlab (comandi help e doc), provare ad utilizzare int, collect,simplify, expand, pretty, ....

Per valutare una funzione numericamente (per esempio sostituire ad x un valore numerico),si utilizza il comando subs nel seguente modo:

>> syms x ;>> y = x ∗ sin ( x ) ;

che puo essere valutata numericamente con il comando subs nel seguente modo:>> subs (y , x , pi /2)ans =

1.5708

Se si sostituisce alla variabile indipendente un vettore, il programma effettua il calcolodell’espressione per ciascuna componente, restituendo un vettore di pari dimensioni:

– 3 –

Prima esercitazione di laboratorio MATLAB Risoluzione simbolica di EDO

>> subs (y , x , 0 : 3 )ans =

0 0.8415 1 .8186 0 .4234

Per visualizzare il grafico di una funzione definita come espressione simbolica, si puo utilizzareil comando ezplot, che fornisce direttamente il grafico della funzione simbolica per l’intervallodi valori inserito sotto forma di vettore nel secondo argomento.

>> e z p l o t (y , 4∗pi ∗[−1 , 1 ] ) ;

3 Risoluzione simbolica di EDO

Se il comando fondamentale per la soluzione simbolica di equazioni e solve, quello per la solu-zione di equazioni differenziali ordinarie e dsolve. Tale comando fornisce l’integrale generale diun’equazione differenziale.

Consideriamo l’equazione differenziale:

dydt = y2.

Vi sono due casi: y(t) ≡ 0 e y(t) 6= 0. In quest’ultimo, procedendo per separazione di variabili siha:

dyy2 = dt ⇒ −1

y+ C = t ⇒ y(t) = − 1

t− C

Per risolvere la stessa equazione con il Toolbox simbolico di Matlab si puo eseguire il seguentecomando:

>> dso lve ( ’Dy = yˆ2 ’ )ans =

0−1/(C3 + t )

Si osservi che, per default, la variabile indipendente viene indicata con la lettera t. Se necessario ilnome della variabile indipendente puo essere modificato (purche sia differente dal nome utilizzatoper indicare la variabile dipendente):

>> dso lve ( ’Dy = yˆ2 ’ , ’x’ )ans =

0−1/(C7 + x )

Si noti, inoltre, come il solutore simbolico sia in grado di individuare anche la soluzioneparticolare y ≡ 0.

Nota: nelle prime versioni di Matlab (prima della R2009b) il toolbox simbolico non era ingrado di trovare entrambe le soluzioni. Ad esempio nella versione R2007a l’output che si ottienee il seguente:

>> dso lve ( ’Dy = yˆ2 ’ )ans =

−1/(C1 + t )

e quindi necessario controllare quale versione del programma si usa e prestare la dovuta atten-zione.

– 4 –

Prima esercitazione di laboratorio MATLAB Risoluzione simbolica di EDO

Per risolvere invece un problema di Cauchy (ossia un’equazione differenziale ordinaria conassegnato un valore iniziale), come ad esempio:{

y′ = y2

y(0) = 1

si puo procedere come segue:>> dso lve ( ’Dy = yˆ2 ’ , ’y (0)=1 ’ )ans =

−1/( t − 1)

>> dso lve ( ’Dy = yˆ2 ’ , ’y (0)=0 ’ )ans =

0

I risultati forniti dal comando dsolve sono a loro volta delle espressioni simboliche che pos-sono essere assegnate come valore a nuove variabili (non necessariamente gia dichiarate) esuccessivamente plottate:

>> y = dso lve ( ’Dy = yˆ2 ’ , ’y (0)=1 ’ )y =

−1/( t − 1)

>> e z p l o t (y , [ 0 , 1 ] ) ;

Con il comando dsolve si possono risolvere anche equazioni differenziali di ordine superioreal primo. Se y′ si traduce in Dy, y′′ sara D2y, y′′′ diviene D3y, e cosı via.

Consideriamo ad esempio l’equazione lineare:

y′′′ + y′′ + y′ + y = 0.

Il suo integrale generale si trova eseguendo:>> y = dso lve ( ’D3y + D2y + Dy + y = 0’ )

y =

C18∗exp(−t ) + C19∗ sin ( t ) + C20∗cos ( t )

Per risolvere invece il problema di Cauchy{y′′′ + y′′ + y′ + y = 0y(0) = 0, y′(0) = 1, y′′(0) = 2

si puo utilizzare il comando>> y = dso lve ( ’D3y + D2y + Dy + y = 0’ , ’y(0) = 0’ , ’Dy (0) = 1’ , ’D2y (0) = 2’ )

y =

exp(−t ) + 2∗ sin ( t ) − cos ( t )

– 5 –

Prima esercitazione di laboratorio MATLAB Risoluzione simbolica di EDO

Esercizio 2

Calcolare con Matlab l’integrale generale delle seguenti equazioni e tracciare il grafico dellesoluzioni per alcuni valori a scelta del parametro C1:

1 y′ + y2 = 0 2 y′ + 2ty = 0 3 y′ + 2ty2 = 0

4 (1 + t)y′ − 4y = 0 5 ty′ + 2y − 3t = 0 6 ty′ + y − t2 = 0

7 y′ + 2ty − t = 0 8 2ty′ − 3y − t3 = 0 9 y′ + 3y − 2te−3t = 0

Esercizio 3

Risolvere con Matlab i seguenti problemi di Cauchy e rappresentare graficamente lesoluzioni:

1{y′ + 3t2y = 0y(0) = 3

2{y′ − y2/t = 0y(1) = 2

3{y′ + 2ty = 0y(0) = 1/2

4{y′ − 1− y2 = 0y(0) = 0

5{y′ − 3 + 2y/t = 0y(2) = 3

6{y′ + y = 0y(0) = 2

7{ty′ + 2y − 3t = 0y(2) = 5

8

y′′ + y′ − sin(t) = 0y(0) = 0y′(0) = 1

9

y′′ + y′ − t− 1 = 0y(0) = 0y′(0) = −1

10

y′′′ + 2y′′ + y′ − t = 0y(0) = 0y′(0) = 1y′′(0) = 2

Esercizio 4

Risolvere con Matlab il problema di Cauchy parametrico{y′ = ay + b cos(ωt)y(0) = c

– 6 –

Prima esercitazione di laboratorio MATLAB Risoluzione simbolica di Sistemi di EDO

4 Risoluzione simbolica di Sistemi di EDO

Il comando dsolve puo essere utilizzato anche per la risoluzione di sistemi di equazioni differen-ziali ordinarie. Si consideri ad esempio il sistema lineare del I ordine{

x′ = y

y′ = −xche e equivalente a d

dt

[xy

]=[

0 1−1 0

] [xy

]la cui soluzione e x = x(t), y = y(t).

In Matlab si puo risolvere nel seguente modo:>> [ x , y ] = dso lve ( ’Dx = y’ , ’Dy = -x’ )

x =

C21∗ sin ( t ) + C22∗cos ( t )

y =

C22∗cos ( t ) − C21∗ sin ( t )

Si noti che in questo caso e fondamentale l’assegnazione della soluzione in una coppia di variabiliper il funzionamento del comando.

Se vengono invece assegnate ad esempio le condizioni iniziali x(0) = 0, y(0) = 1 si ha:>> [ x , y ] = dso lve ( ’Dx = y’ , ’Dy = -x’ , ’x(0) = 0’ , ’y(0) = 1’ )

x =

sin ( t )

y =

cos ( t )

Esercizio 5

Calcolare usando Matlab le soluzioni dei seguenti sistemi di equazioni differenziali ordinarie,rappresentando poi graficamente le soluzioni. Scrivere inoltre le equazioni in forma matriciale.

1

x′ = 3x− 2yy′ = 2x+ y

x(0) = 1y(0) = 0

2

x′ = 2x+ 4y + 3et

y′ = 5x− y − t2

x(0) = 0y(0) = −1

3

x′ = 3x− 2yy′ = −x+ 3y − 2zz′ = −y + 3zx(0) = 0y(0) = 0z(0) = 4

– 7 –

Prima esercitazione di laboratorio MATLAB Limiti dei solutori simbolici di EDO

Esercizio 6

Si trasformi il problema di Cauchy parametrico:mx′′ + γx′ + kx = a cos(ωt)x(0) = α

x′(0) = β

in un sistema di due equazioni e se ne calcoli la soluzione con Matlab.

5 Limiti dei solutori simbolici di EDO

Il fatto che sia possibile calcolare la soluzione esatta di alcune (semplici) equazioni differenzialinon deve indurre a pensare che per tale tipo di problema possa sempre essere calcolata lasoluzione mediante una formula esplicita.

In generale, e vero esattamente il contrario. Per la maggior parte delle equazioni differen-ziali ordinarie e dei sistemi, anche qualora l’esistenza e l’unicita della soluzione siano garantitedall’applicabilita dei relativi teoremi, non esistono soluzioni facilmente calcolabili ed esprimibilimediante formule chiuse.

Un esempio tipico di questa situazione e il seguente problema di Cauchy:{y′ = e−t2

, t ∈ (0, 3]y(0) = 1

L’equazione differenziale associata a questo problema e concettualmente banale, dato che equivalead una semplice integrazione definita e la soluzione non e altro che

y(t) = 1 +∫ t

0e−s2

ds, t ∈ [0, 3].

D’altra parte, dato che l’integrale definito della funzione gaussiana non puo essere espressomediante una semplice formula chiusa, tale formula non e altro che una rappresentazione dellasoluzione che deve poi essere approssimata numericamente qualora si vogliano effettivamentecalcolare i valori della soluzione al variare di t (cio non deve spaventare, anche le funzioni sin, cos,log devono essere approssimate numericamente dal calcolatore, in quando non sono esprimibiliattraverso le operazioni elementari di somma e prodotto).

In effetti in Matlab otteniamo:>> dso lve ( ’Dy = exp(-tˆ2) ’ , ’y (0)=1 ’ )

ans =

( pi ˆ(1/2)∗ erf ( t ) )/2 + 1

La funzione erf e un oggetto matematico ben noto in statistica (non a caso si trova anche inquasi tutte le calcolatrici tascabili).

Dato che questo problema non deriva semplicemente dalla complessita del calcolo da eseguire,ma da un ostacolo concettuale, non sara possibile in casi come questo (che, ripetiamo ancora,costituiscono la regola se si eccettuano particolari classi di equazioni come quelle lineari ed

– 8 –

Prima esercitazione di laboratorio MATLAB Introduzione alla risoluzione numerica di EDO

omogenee a coefficienti costanti di primo e secondo ordine) non sara possibile ottenere unasoluzione neanche con Matlab. Si otterranno invece dei messaggi d’errore di svariato tipo, odelle soluzioni formali in termini di cosidette funzioni speciali, che debbono poi a loro volta esserecalcolate e approssimate numericamente.

Nel caso di equazioni di ordine superiore al primo, pur restringendosi al caso di equazio-ni lineari, a tali difficolta si deve aggiungere quella legata alla determinazione delle radici delpolinomio caratteristico associato alla equazione data.

La determinazione di tali radici e possibile attraverso la modalita simbolica solo per alcunipolinomi piu facilmente trattabili tramite semplici trasformazioni algebriche. In generale, peruna arbitraria equazione di ordine abbastanza alto (ad esempio oltre il quarto), Matlab nonsara in grado di determinare una soluzione lavorando in modalita simbolica. Considerazioni deltutto analoghe possono essere fatte per i sistemi, per cui non e in generale conveniente cercare dideterminare la soluzione di un sistema di piu di due equazioni differenziali lavorando in modalitasimbolica.

If no closed-form (explicit) solution is found, an implicit solution is attempted.When an implicit solution is returned, a warning is given.

If neither an explicit nor implicit solution can be computed, then a warning is givenand the empty sym is returned.

In some cases involving nonlinear equations, the output will be an equivalent lowerorder differential equation or an integral3.

6 Introduzione alla risoluzione numerica di EDO

Come si e visto nella precedente sezione, in generale non e quasi mai possibile calcolare esatta-mente la soluzione di una equazione differenziale ordinaria. Esistono pero degli algoritmi che neforniscono delle soluzioni approssimate e che possono essere implementati (anche) in Matlab.Lo studio approfondito di tali algoritmi viene svolto nell’ambito dei corsi di Calcolo Numerico eAnalisi Numerica.

In questa esercitazione ci si limitera ad introdurre una semplice implementazione del metododi Eulero, mentre nella seguente esercitazione si utilizzeranno metodi numerici piu sofisticati epiu accurati implementati in funzioni predefinite in Matlab.

Consideriamo il generico problema di Cauchy (scalare):{y′(t) = f(t, y(t)), t ∈ (0, T ]y(0) = y0,

con y : [0, T ] → R, f : (0, T ] × R → R e y0 ∈ R. L’equazione ci dice che geometricamente latangente di y(t) e esattamente f(t, y(t)), per ogni t ∈ (0, T ]. Quindi, la cosa piu ovvia che sipuo fare, dato un valore noto di y(t) (per esempio il dato iniziale y(0)), e procedere passo perpasso lungo la tangente (che fornisce una direzione). Piu precisamente, se supponiamo di avere adisposizione al tempo tk il valore approssimato uk, possiamo determinare uk+1 facendo un passodi dimensione hk sommando ad uk il valore hk · f(tk, uk).

3Stralcio dell’help di Matlab

– 9 –

Prima esercitazione di laboratorio MATLAB Introduzione alla risoluzione numerica di EDO

In formule, il metodo di Eulero esplicito (o in avanti) si scrive come:{uk+1 = uk + hk · f(tk, uk), tk ∈ (0, T ], k = 1, 2, . . .u0 = y0.

Il vettore{uk

}k=kmax

k=0 sara l’approssimazione numerica di y(t), ossia uk approssimera y(tk)con un certo margine di errore. Se hk = h costante, allora tk = h · k e k variera tra 0 e l’interopiu vicino a T/h.

Vogliamo ora implementare l’algoritmo in Matlab: prima pero apriamo una parentesi sullefunzioni definite da utente, che torneranno molto utili anche nella prossima esercitazione.

6.1 Script e Funzioni in MATLAB

Fino ad ora abbiamo utilizzato Matlab in modalita interattiva, ossia immettendo i comandiuno alla volta direttamente nell’interfaccia grafica (il prompt dei comandi) e aspettando, di voltain volta, l’output. Matlab ci dice quando e in attesa di comandi con il simbolo >> nel prompt.

E buona norma pero utilizzare script e funzioni per operazioni piu strutturate e complesse(come, ad esempio, risolvere un’equazione differenziale ordinaria) o anche solo per fare una seriedi comandi ripetitivi. Uno script non e nient’altro che un file di testo che contiene una sequenzadi operazioni che altrimenti potreste scrivere uno alla volta al prompt dei comandi. Una funzioneinvece e un particolare tipo di script con un’intestazione ben precisa (la firma della funzione)necessaria per essere richiamata.

Script e funzioni possono essere scritti direttamente con l’editor di Matlab:>> e d i t

oppure con un qualsiasi editor di testo (come Notepad su sistemi Windows).Apriamo dunque l’editor di script e scriviamo al suo interno:

1 function dy = funz ione ( t , y )2

3 % Questa funz ione c a l c o l a s i n ˆ2(y ) + t4

5 dy = sin ( y ) . ˆ 2 + t ;

Listing 1: Funzione che valuta f(t, y) = sin2(y) + t.

Salviamo il file con il nome funzione.m (il nome del file che contiene la funzione, suggerito daMatlab stesso, deve essere identico al nome della funzione).

Ora dal prompt dei comandi digitiamo:>> funz ione (2 , pi /2)ans =

3

>> help funz ione

Questa funz ione c a l c o l a s i n ˆ2(y ) + t

Operativamente Matlab richiama la funzione passando i valori t = 2 e y = pi/2. Le linee dicommento (non obbligatorie ma utili) sono stampate a video come guida della funzione attraversoil comando help. Come vedremo, molte funzioni di Matlab prendono in ingresso a sua voltauna funzione, la quale pero deve avere una firma ben precisa.

– 10 –

Prima esercitazione di laboratorio MATLAB Introduzione alla risoluzione numerica di EDO

Nota: Nelle ultime versioni di Matlab (superiori alla R2008a), se un parametro deve esserepresente ma non e utilizzato, e buona norma sostituirne il nome con il simbolo4 ∼. Ad esempio(nel caso della versione R2007a questo non funzionera, quindi si inserisca t come nell’esempioprecedente):

1 function dy = fun (∼ , y )2 dy = y.∗(1−y ) ;

Listing 2: Funzione che valuta f(t, y) = y · (1− y).

ossia nient’altro che il membro a destra dell’equazione della logistica: y′ = y · (1− y). Si ottiene:>> fun ( pi , 2)ans =

−2

solo che questa volta MATLAB non passa π alla funzione, guadagnando qualcosa in efficienza.

6.2 Implementazione del metodo di Eulero

Iniziamo introducendo il tipo di problema che vogliamo affrontare. Inplementiamo la funzionefunction dy = fun ( t , y )dy = y.∗(1−y ) ;

ossia nient’altro che il membro a destra dell’equazione della logistica: y′ = y · (1− y). Si ottiene:>> fun ( pi , 2)ans =

−2

Il metodo di Eulero dipende dai parametri dell’equazione e dai parametri di discretizzazione,come il passo temporale h e il tempo totale T . Fissiamoli dunque una volta per tutti:

1 %% Dati i n i z i a l i2 t0 = 0 ;3 y0 = 0 . 1 ;4 %% Parametri d i d i s c r e t i z z a z i o n e5 T = 1 0 . 0 ; % Tempo t o t a l e s imu laz ione6 h = 0 . 5 ; % Passo temporale

Abbiamo poi bisogno di un vettore ove memorizzare la soluzione e i passi temporali:7 %% Memorizzazione s o l u z i o n e8 t = t0 : h :T;9 N = s ize ( t , 2 ) ; % Numero d i p a s s i t empora l i che e f f e t tue r emo .

10 % R i c o r d a r s i che s i z e ( t ) r e s t i t u i s c e11 % l e dimens ion i d i t in formato v e t t o r e12 % ( [ 1 ,N] in questo caso ) e che s i z e ( t , 2)13 % f o r n i s c e l a seconda componente ( o s s i a N) .14 u = zeros (1 , N) ; % Vettore contenente l a s o l u z i o n e : y ( t ( k ) ) ∼ u( k )

Procediamo ora con l’algoritmo vero e proprio:15 %% I n i z i a l i z z a z i o n e16 u (1) = y0 ; % Assegno i l dato i n i z i a l e17 %% Pass i tempora l i

4Per digitare il simbolo “∼” su tastiere italiane si deve tener premuto il tasto Alt e digitare il numero 126 sultastierino numerico. Il simbolo comparira una volta rilasciato Alt . Ricordarsi di attivare Num . Su sistemiUNIX-like, come GNU/Linux, basta tener premulto AltGr e premere il tasto della “ı” (“i” accentata, vicino a7−→ ).

– 11 –

Prima esercitazione di laboratorio MATLAB Introduzione alla risoluzione numerica di EDO

18 for k=1:N−119 % Suppongo che fun s i a d e f i n i t a20 u( k+1) = u( k ) + h ∗ fun ( t ( k ) , u ( k ) ) ;21 end

Per convincerci del risultato plottiamo il risultato sovrapponendolo alla soluzione esatta:22 %% Soluz ione e s a t t a23 % Scr ivo l e c o n d i z i o n i i n i z i a l i assemblando l a s t r i n g a co r r i spondente24 i n i t c ond = [ ’y(’ ,num2str( t0 ) , ’)= ’ ,num2str( y0 ) ] ; % Cost ru i s co y ( t0)=y025 % Riso lvo s imbol icamente con s v a r i a b i l e ind ipendente26 y = dso lve ( ’Dy = y * (1-y)’ , in i t cond , ’s’ ) ;27 % S o s t i t u i s c o s s imbo l i ca con i l v e t t o r e s , e p l o t t o in t r a t t o cont inuo nero28 s = linspace ( t0 , T, 1000 ) ; % l i n s p a c e (a , b ,N) d iv id e i l segmento29 % [ a , b ] in N e lement i30 c l f ;31 plot ( s , subs (y , ’s’ , s ) , ’-k’ ) ;32 %% Soluz ione appross imata33 hold on ; % Mantengo l a f i n e s t r a g r a f i c a co r r en t e aperta34 plot ( t , u , ’*-r’ ) ;35 grid on ;36 legend ( ’Sol. esatta ’ , ’Sol. numerica ’ ) ;37 t i t l e ( ’Metodo di Eulero in avanti ’ ) ;38 xlabel ( ’t’ ) ;

0 1 2 3 4 5 6 7 8 9 100

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

t

y

Metodo di Eulero in avanti

Sol. esatta

Sol. numerica

Figura 1: Soluzione esatta ed approssimata sovrapposte.

Si osserva la classica “sigmoide” o curva logistica, fondamentale in dinamica delle popolazioni.

Esercizio 7

Studiare la soluzione del problema seguente:{y′ = y(y − 1)(0.5− y)y(0) = y0

al variare di y0. Questa e detta equazione bistabile, in quanto vi sono due equilibri stabili (cioeverso i quali le traiettorie convergono) e uno instabile a meta. E un’equazione molto importantein chimica e biochimica (nonche in fisiologia).

– 12 –

Prima esercitazione di laboratorio MATLAB Introduzione alla risoluzione numerica di EDO

Esercizio 8

Risolvere il seguente problema di Cauchy:{y′ = 100y(y − 1)(0.5− y)y(0) = 0.6

per T = 5.0 al variare di h tra 0.01 e 0.1. Cercare di descrivere come evolve la situazione.

6.3 Funzioni anonime

Puo essere utile a volte utilizzare le funzioni anonime, ad esempio nel caso in cui un solutore(come il nostro programma di Eulero esplicito) richiede sempre una funzione di due argomentima noi ne abbiamo scritta una che dipende anche da parametri aggiuntivi.

Modifichiamo la funzione fun nel seguente modo:function dy = funpar ( t , y , a , b )%FUNPAR Funzione l o g i s t i c a d i parametr i a e b%% FUNPAR( t , y , a , b ) c a l c o l a a∗y∗(b−y )dy = a∗y . ∗ ( b−y ) ;

Se provassimo ad eseguire nuovamente il solutore di Eulero riceveremmo un errore da parte diMatlab. In effetti basterebbe modificare la riga 16 mettendo funpar (t(k), u(k), a, b) alposto di fun (t(k), u(k)). Ma non e cio che vogliamo fare (si immagini la situazione in cui loscript e fornito ma non modificabile).

Il trucco che possiamo usare sfrutta l’idea seguente:>> a = 1 . 0 ;>> b = 1 . 0 ;>> fun = @( t , y ) funpar ( t , y , a , b ) ;>> funpar ( 0 . 0 , 0 . 5 , a , b )

ans =

0.2500

>> fun ( 0 . 0 , 0 . 5 )

ans =

0.2500

All’atto pratico, quando chiamo funpar, Matlab chiama fun sostituendo automaticamente ilvalore d’ambiente dei parametri a e b. In questo modo possiamo usare lo script di Eulero senzamodificare nulla: basta ridefinire fun.

6.4 Accuratezza del metodo

Una volta deciso il metodo, il passo successivo e capire quanto questo sia “buono”. Scelto lo sche-ma rimangono i parametri di discretizzazione da selezionare. Come si sceglie, dunque, il passotemporale? C’e quale vincolo? Quanto migliora la soluzione al diminuire di h? Mi avvicino davve-ro indefinitivamente alla soluzione vera? Queste domande portano a tre importantissimi concetti

– 13 –

Prima esercitazione di laboratorio MATLAB Introduzione alla risoluzione numerica di EDO

0 1 2 3 4 5 6 7 8 9 100

0.2

0.4

0.6

0.8

1

1.2

1.4Metodo di Eulero in avanti

t

y

Sol. esatta

Sol. numerica

(a) h = 2.0

0 1 2 3 4 5 6 7 8 9 100

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

t

y

Metodo di Eulero in avanti

Sol. esatta

Sol. numerica

(b) h = 1.0

0 1 2 3 4 5 6 7 8 9 100

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

t

y

Metodo di Eulero in avanti

Sol. esatta

Sol. numerica

(c) h = 0.5

0 1 2 3 4 5 6 7 8 9 100

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

t

y

Metodo di Eulero in avanti

Sol. esatta

Sol. numerica

(d) h = 0.1

Figura 2: Approssimazione numerica al variare del passo temporale h.

in Analisi Numerica: stabilita, convergenza e consistenza che saranno trattati profusamente inesami successivi e che qui abbiamo intenzione solo di sfiorare.

In Figura 2 possiamo vedere come varia l’approssimazione calcolata dal metodo di Eulero alvariare di h: pare che vi sia convergenza (si puo mostrare che scala linearmente con h).

Gli errori possono essere calcolati relativamente all’istante finale T = tN o relativamente atutto l’intervallo [0, T ]. L’errore assoluto all’istante finale puo essere definito come

errabs(T ) = |y(tN )− uN |.

Tenendo conto pero che i valori numerici della soluzione esatta potrebbero essere molto grandi omolto piccoli, ci si rende conto che e piu opportuno rapportare la dimensione dell’errore a quelladella soluzione stessa, introducendo quindi il cosidetto errore relativo

errrel(T ) = |y(tN )− uN ||y(tN )| .

Se si vuole considerare invece l’errore compiuto su tutto l’intervallo [0, T ] si deve tenereconto anche degli istanti precedenti all’ultimo. Una possibile definizione e quella che consiste nel

– 14 –

Prima esercitazione di laboratorio MATLAB Introduzione alla risoluzione numerica di EDO

prendere il massimo errore commesso su tutto l’intevallo (detto errore in norma `∞):

errabs∞ (0, T ) = max

k=0,...,N|y(tk)− uk|.

Il corrispondente errore relativo e invece

errrel∞ (0, T ) = maxk=0,...,N |y(tk)− uk|

maxk=0,...,N |y(tk)| .

Tale misura di errore considera il massimo scarto rilevato tra la soluzione esatta e quella numericain tutti gli istanti in cui si e proceduto a calcolare l’approssimazione uk.

Un’altra possibile definizione e invece quella del cosidetto errore quadratico medio, analogoa quello che si vede in statistica (detto anche errore in norma `2):

errabs2 (0, T ) =

√√√√ N∑k=0

h · |y(tk)− uk|2.

Il corrispondente errore relativo e invece:

errrel2 (0, T ) =

√√√√∑Nk=0 h · |y(tk)− uk|2∑N

k=0 h · |y(tk)|2.

Si noti che nella definizione dell’errore relativo il termine∑N

k=0 |y(tk)|2h costituisce di fatto unaapprossimazione dell’integrale

∫ T

0 |y(s)|2ds.

Esercizio 9

Calcolare gli errori assoluti e relativi compiuti nell’approssimazione della soluzione delproblema di Cauchy {

y′ + y = 2 cos(3t)y(0) = 1

sull’intervallo [0, 6] utilizzando il metodo di Eulero con un passo h = 0.1, h = 0.01, h = 0.001.

– 15 –