LIBRERIA C++ PER L'INTEGRAZIONE DI SISTEMI DI ODE...

Politecnico di Milano

Facoltà di Ingegneria dei sistemi

Corso di Laurea in Ingegneria matematica

Progetto di Programmazione Avanzata per il Calcolo Scientifico

LIBRERIA C++ PER L'INTEGRAZIONE

DI SISTEMI DI ODE DRH

di

Massimiliano Beccaria, Giorgio Sandrini

Matricola 754742�750097

Anno Accademico 2011�2012

Contents

1 Teoria per sistemi di ODE DRH . . . . . . . . . . . . . . . . . . . . . . 2

1.1 De�nizione del problema . . . . . . . . . . . . . . . . . . . . . . . 2

1.2 Analisi di Filippov . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.3 Condizioni di ordine superiore . . . . . . . . . . . . . . . . . . . . 3

1.4 Ricerca punti di transizione . . . . . . . . . . . . . . . . . . . . . 4

1.5 Proiezione su soglia . . . . . . . . . . . . . . . . . . . . . . . . . 5

2 Struttura del codice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.1 Inclusioni e test . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.2 Libreria per sistemi di ODE . . . . . . . . . . . . . . . . . . . . . 6

2.3 Libreria per sistemi di ODE-DRH . . . . . . . . . . . . . . . . . 11

3 Casi test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.1 Libreria ODE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.2 Libreria ODEDRH-Autonomi . . . . . . . . . . . . . . . . . . . . 20

3.3 Libreria ODEDRH-Non Autonomi . . . . . . . . . . . . . . . . . 23

3.4 Libreria ODEDRH-Ordine superiore . . . . . . . . . . . . . . . . 26

Bibliogra�a 28

1

1 Teoria per sistemi di ODE DRH

1.1 De�nizione del problema

La libreria sviluppata ha come obiettivo un trattamento esaustivo dei sistemi di equazioni

di�erenziali ordinarie con forzanti discontinue (ODE-DRH) in cui le discontinuità si

presentano su una o più super�ci (super�ci di passaggio). In generale il problema modello

è il seguente:

dy

dt=

f1(t,y), se g(t,y) < 0

f2(t,y), se g(t,y) > 0(1)

Allo scopo di trattare la presenza di tali super�ci di discontinuità, de�nibili a livello

matematico come la curva di livello zero di una funzione nelle variabili di stato del

sistema e nel tempo (g(t,y) = 0), ci si baserà sulla teoria di Filippov [1] estesa alle con-

dizioni di ordine superiore da Dieci-Lopez nell'articolo "Sliding motion in discontinuous

di�erential systems: Theory and a computational approach" [2]. Inoltre, a di�erenza

di quanto fatto in tali lavori, sarà considerata anche la possibilità che il sistema, come

evidenziato in (1), possa essere non autonomo sia nelle forzanti che nella funzioni di

soglia.

1.2 Analisi di Filippov

Sia la variabile di stato yi appartenente alla funzione di soglia gi(t,y) = 0 e de�niamo

con n(t,y) la normale a tale super�cie. Essa potrà essere calcolata come:

n(t,y) =∇gi(t,y)‖∇gi(t,y)‖

. (2)

De�niamo la variabile ausiliaria γ1,2(t,y(t)) = ∂g/∂t+∇g · f1,2 dove si è omesso l'indice

i della soglia attiva ma si considera sempre il campo di forzanti sotto (γ1) e sopra (γ2)

tale soglia. I termini γ saranno le proiezioni dei campi di forzanti corrispondenti sulla

normale alla soglia attiva. Sia y appartenente ad una super�cie di discontinuità del

sistema, si de�nisce A come:

A1,2 = γ1,2(t,y(t)). (3)

La analisi di Filippov de�nisce le seguenti condizioni per valutare le interazioni con le

super�ci di discontinuità delle variabili di stato:

1. se (A1) · (A2) > 0, la soluzione lascia la soglia (crossing). Essa entra nella regione

sotto soglia (S1) se (A1) < 0, entra nella regione sopra soglia (S2) se (A2) >

0. Il campo di forzanti con cui si proseguirà sarà corrispondente alla regione di

destinazione.

2. se (A1) · (A2) < 0 e (A1) > 0 e (A2) < 0 si ha una soluzione che rimane sulla soglia

(sliding motion). Inoltre tale soglia risulta essere attrattiva e le variabili di stato

2

saranno soggette ad un campo di forzanti fF che permetterà di proseguire nel suo

percorso rimanendo "ancorata" alla soglia. Tale campo è de�nibile come:

fF (t,y) = (1− α(t,y))f i1(t,y) + α(t,y)f i2(t,y). (4)

dove i è l'indice della variabile coinvolta nella sliding motion e:

α(t,y) =∂gi(t,y)/∂t+ (nT f i1(t,y))

nT (f i1(t,y))− f i2(t,y)). (5)

3. se (A1) · (A2) < 0 e (A1) < 0 e (A2) > 0 si ha uno sliding motion repulsivo.

Tale condizione risulta instabile e l'unicità della soluzione non è garantita perchè

le variabili di stato possono passare spontaneamente in S1 o S2.

Tale teoria non è esaustiva delle possibili condizioni che possono veri�carsi per un

problema di ODE−DRH in quanto non considera i casi in cui i termini A1,2 assumano

valori nulli. Per poter trattare anche tali casi si utilizza la teoria di Dieci-Lopez per

condizioni di ordine superiore sempre estesa al caso generico non autonomo.

1.3 Condizioni di ordine superiore

Sia y un set di variabili di stato che risulta sulla soglia gi(t,y) = 0 o passante da essa

al tempo t e de�niamo i seguenti termini ausiliari di analisi derivanti dallo sviluppo al

secondo ordine delle funzioni γ de�nite in precedenza:

B±1,2 =

∂γ1,2∂t

+∇γ1,2(t,y(t))y′(t)

t±

(6)

C±1,2 =

∂2γ1,2∂t2

+∇γ1,2(t,y(t))y′′(t) + y′T (t)∂γ1,2

∂y2y′(t)

t±

(7)

Ora se A1 6= 0 e A2 6= 0 si applica l'analisi al primo ordine di Filippov. Al contrario,

i casi in cui A1 = 0 e/o A2 = 0 determinano un andamento della soluzione che può

essere descritto in funzione di A,B e C:

1. A1 = 0 e A2 6= 0. Se si arriva da uno sliding motion attrattivo (A2 < 0) si dimostra

che necessariamente si ha B−1 ≤ 0. In tale condizione, si possono veri�care i

seguenti andamenti delle variabili di stato:

• B−1 < 0 allora la soluzione lascia la soglia e va in S1 con il campo corrispon-

dente.

• B−1 = 0 allora si dimostra che C−1 ≥ 0. Inoltre, lo sliding motion continua

con un campo che è corrispondente a quello in S1.

Se la soluzione arriva da sotto soglia si hanno le seguenti possibilità:

3

• A2 < 0, allora si applica la analisi del punto precedente.

• A2 > 0 allora si ha uno sliding motion repulsivo ed instabile se B−1 < 0, una

soluzione di crossing in S2 la soglia se B−1 = 0 e C−1 > 0.

Se si arriva da S2 allora A2 < 0, y′(t−) è uguale al campo sopra soglia e y′(t+) è

pari al campo sotto soglia. In tal caso si hanno le seguenti possibilità:

• B+1 < 0, allora c'è soluzione di crossing in S1 con il campo corrispondente di

forzanti.

• B+1 > 0 allora si soluzione di sliding che rimane con un campo corrispondente

ad S1.

2. A2 = 0 e A1 6= 0. L'analisi è la stessa di prima con gli indici 1 e 2 ribaltati.

3. A1 = 0 e A1 = 0. In tale condizione si può arrivare a soluzione molto complesse e

non del tutto de�nibili. Tuttavia si hanno dei risultati parziali per casi nei quali

si può decidere che fare. In particolare se B−1 = B−2 = 0 lo sliding continuerà, se

B−1 < 0,B−2 = 0 e C−2 < 0 c'è crossing in S1, se B−1 = 0,B−2 > 0 e C−1 > 0 c'è

crossing in S2, se B−1 < 0 e B−2 > 0 c'è una soluzione instabile e non unica.

1.4 Ricerca punti di transizione

Ci si pone nel caso in cui le variabili di stato incontrino una soglia nella loro evoluzione.

Le condizioni andranno testate al tempo e nel punto preciso di transizione. In generale è

molto probabile che l'integrazione numerica fornisca una soluzione che, a distanza di un

passo temporale tra un nodo e l'altro, appartenga a regioni opposte rispetto alla soglia.

E' quindi necessario de�nire un metodo per localizzare con precisione il punto in cui

si cade sulla discontinuità. Inoltre tale localizzazione, se esatta, permette di preservare

l'ordine di convergenza dello schema numerico in presenza di DRH.

Immaginiamo per esempio di avere nella nostra soluzione numerica y(n) che appar-

tiene ad una regione S1 sotto soglia e y(n+1) che appartiene ad una regione S2 sopra

soglia. Si vuole cercare il tempo t∗ a cui si arriva sulla soglia e la corrispondente soluzione

y(t∗). Tale problema può essere formulato matematicamente come la ricerca dello zero

della funzione G(η) = g(η,y(η)). Usando il metodo delle secanti si ottiene il seguente

algoritmo iterativo di risoluzione:

ηk+1 = ηk −ηk − ηk−1

G(ηk)−G(ηk−1))G(ηk) (8)

inizializzato con η0 = t(n) e η1 = t(n+1). Nella libreria per trovare il t∗ si utilizza tale

metodo con un criterio di arresto sulla di�erenza tra un passo ed il successivo. Si noti

inoltre che per ogni ηk provvisorio si deve risolvere il corrispondente problema �no a

tale valore temporale per ottenere G(ηk).

4

Può essere conveniente de�nire il t∗ di transizione delle variabili di stato anche quando

si veri�cano le condizioni di uscita da una soglia. In tal caso, per campi continui e regolari

al di fuori delle super�ci di discontinuità, il vantaggio relativo al calcolo del t∗ esatto

sarà minore in quanto il campo passando da uno stato di sliding attrattivo di�cilmente

avrà valori che si discosteranno molto dalla tangente alla soglia. In ogni caso tale calcolo

di t∗ è implementato nella libreria nel caso in cui le condizioni del primo ordine siano

utilizzabili.

Ipotizziamo per esempio che A1 > 0 e A2 < 0 ad un passo temporale ci si trova

in una condizione di sliding motion attrattivo e che al successivo si abbia A1 > 0 e

A2 > 0 e quindi si debba uscire dalla soglia verso il campo superiore ad essa (S2). Per

ottenere il tempo preciso di uscita si cerca lo zero della funzione A2(t,y(t)). Il procedi-

mento numerico per trovare lo zero di tale funzione è il medesimo del caso precedente.

Ovviamente, nel caso si vada verso S1, si ricercherà lo zero di A1(t,y(t)).

1.5 Proiezione su soglia

Ipotizziamo di essere in condizione di sliding attrattivo su di una soglia attiva. A livello

numerico, seppur utilizzando il campo corretto per rimanere vincolati alla soglia, le

soluzioni ai passi successivi potranno discostarsi da essa. Per questo si implementa un

algoritmo di proiezione sulla soglia attiva.

Sia y un punto in prossimità di una soglia g = 0 a tempo �ssato, allora la sua

proiezione y = P (y) sarà la soluzione del problema di minimo:

ming(y)=0

g(y) (9)

con g(y) =1

2(y − y)T (y − y).

Tale problema di minimo verrà trattato attraverso moltiplicatore di lagrange λ per

il vincolo e risolto con un metodo iterativo di Newton sempli�cato per evitare di dover

risolvere un vero sistema lineare ad ogni iterazione [2]. Questo è legittimo in quanto ci

aspettiamo di essere in ogni caso vicini alla soluzione della minimizzazione.

5

2 Struttura del codice

Il codice è strutturato su due librerie. La prima libreria risolve sistemi di ODE classici

senza il trattamento di discontinuità. Essa risulta del tutto indipendente dalla seconda

e può essere usata come libreria autonoma. La seconda libreria risolve sistemi di ODE

DRH. Tale seconda libreria importa i �le relativi alla prima e utilizza i solutori ivi

implementati per avanzare in tempo con il campo di forzanti de�nito dall'analisi per il

trattamento delle discontinuità presentata nella sezione precedente.

2.1 Inclusioni e test

La libreria per la risoluzione di sistemi di ODE classici senza il trattamento di discon-

tinuità importa strutture dati e funzioni delle librerie Boost (versione utilizzata 1.49.0)

ed Eigen (versione utilizzata 3.0.5). In particolare le Eigen sono state utilizzate per

la risoluzione di sistemi lineari con il metodo di Newton nel caso si utilizzino metodi

impliciti di risoluzione. La libreria per sistemi di ODE DRH ovviamente importando

la precedente richiede essa stessa le medesime librerie. Essa inoltre usa direttamente le

Boost. In�ne entrambe utilizzano per il passaggio dei dati una libreria di parser dis-

tribuita assieme ad esse ed utilizzata anche in LifeV (i dati potranno essere passati con

due modalità di�erenti spiegate più nel dettaglio in seguito).

Il codice è stato testato con ValGrind per correggere la presenza di memory leaks.

2.2 Libreria per sistemi di ODE

La prima libreria implementa i metodi di Eulero esplicito, Eulero implicito, Runge-

Kutta del quarto ordine, Heun, Crank-Nicolson, Adams-Bashforth e Adams-Moulton

per la risoluzione di sistemi di equazioni di�erenziali ordinarie. Si riporta in �gura 1 il

corrispondente diagramma UML.

Input/Output

La libreria richiede in input la �nestra temporale in cui si vuole approssimare numeri-

camente la soluzione del sistema (tempo iniziale e tempo �nale), il passo temporale che

si vuole utilizzare, il numero di variabili, le forzanti del sistema, le condizioni iniziali e

il nome del �le su cui si vogliono salvare i risultati. Si deve inoltre speci�care il metodo

desiderato per la risoluzione.

Tali variabili di ingresso saranno passabili sia attraverso un �le di testo esterno che

direttamente dall'interno del codice. Ovviamente la lettura da �le permette di non dover

ricompilare il codice se si cambiano i parametri di simulazione. Per quanto riguarda le

forzanti esse, nel caso di passaggio da �le, saranno passate come stringhe e trattate medi-

ante il parser. Nel caso di passaggio da codice interno, saranno trattate mediante funtori

(puntatori a funzioni). Come si spiegherà più nel dettaglio nella prossima sottosezione,

6

Figure 1: Diagramma UML per la libreria per sistemi di ODE.

il codice è fatto in modo tale che, la classe che implementa la risoluzione del sistema

risulti del tutto trasparente al fatto che le forzanti siano trattate in uno o nell'altro

modo. Tale di�erente trattazione verrà implementata esternamente a tale classe.

Se il metodo richiesto è implicito sarà necessario passare anche le derivate esatte delle

forzanti rispetto alle variabili coinvolte. Sempre per metodi impliciti è possibile de�nire

i parametri di tolleranza e numero massimo di iterazioni da utilizzare per il metodo

di Newton implementato anch'esso ex novo. Se tali parametri non vengono de�niti si

utilizzeranno valori di default. Ancora una volta, anche le derivate sono passabili sia

direttamente mediante �le che dal codice.

Per metodi multistep (Adams-Bashforth e Adams-Moulton) sarà possibile cambiare

il metodo ad un passo da usare per i primi passi del metodo (quando non sono a dispo-

sizione abbastanza soluzioni precedenti per usare il corrispondente multistep).

Nel momento in cui viene lanciata la simulazione numerica per la risoluzione del

sistema, viene controllato che i dati forniti siano corretti (il tempo iniziale minore del

�nale, se implicito che si siano caricate le derivate....) e se i dati non sono corretti viene

scritto a video un messaggio di errore di�erente a seconda del tipo di dato non corretto

e la simulazione si arresta.

In output la libreria restituisce un �le di testo con i valori delle soluzioni ed un

�le GnuplotView che può essere lanciato da main per creare direttamente un gra�co

dell'andamento delle variabili di stato in tempo. Vengono stampate a video le soluzioni.

Inoltre alcune funzioni ausiliarie permettono di esportare la struttura dati contenente le

soluzioni come variabile.

7

Speci�che codice

La classe centrale della libreria è Solver che implementa il risolutore astratto per il

sistema di ODE. In essa sono de�nite le strutture dati delle variabili coinvolte. In parti-

colare le soluzioni vengono memorizzate in un ublas vector di pile di double. La scelta di

tale struttura deriva dal fatto che i metodi risolutivi richiederanno per la maggior parte

la soluzione al passo precedente per l'avanzamento in tempo o �no a tre passi precedenti

per i metodi multistep; per questo motivo è stata scelta la struttura a pila che rende

disponibili gli ultimi valori immessi.

Le forzanti e le eventuali derivate sono salvate come ublas vector di puntatori ad una

classe functor. Tale classe è astratta e de�nisce l'operatore () per la valutazione del valore

della funzione al variare delle variabili spaziali e della variabile temporale. L'operatore

() utilizza una funzione GetValue che è la funzione de�nita pure virtual nella classe

functor. Le classi RealFunctor e ParsedFunctor ereditano dalla classe functor astratta,

de�niscono la corretta struttura dati a seconda rispettivamente del caso di passaggio

da funtori e di passaggio da stringa con Parser e fanno l'override corrispondente per la

funzione GetValue. Questo permette nella classe Solver di usare genericamente punta-

tori a functor che a seconda del tipo di passaggio corrisponderanno a puntatori ad uno o

all'altro tipo concreto. Visto che l'interfaccia all'utente delle due classi concrete derivate

da functor è la stessa, l'implementazione della classe Solver è del tutto trasparente al

tipo di funzioni passate. La scelta di mantenere due tipi di passaggi ha queste moti-

vazioni: in primo luogo, il passaggio da �le mediante Parser è certamente più comodo

per l'utilizzatore ma se la libreria venisse usata in un codice più esteso in cui i funtori

sono già disponibili per altri scopi (come avverrà anche nel nostro caso per la libreria

per ODE-DRH) è utile non costringere l'utente a riscriverli sotto forma di stringhe o di

scrivere un �le apposito; secondo, se il passaggio da �le è più comodo è altresì vero che

una funzione passata come funtore richiederà un tempo per la sua valutazione in coor-

dinate speci�che minore rispetto ad un'altra passata con Parser. Visto che nel codice

vanno fatte molte valutazioni delle forzanti, ciò risulta in un tempo di esecuzione minore

per il passaggio con funtori che in casi complessi può risultare molto utile.

Il metodo centrale della classe Solver è DoIt che permette di eseguire l'integrazione

numerica. Esso contiene la struttura comune a tutti i metodi per il calcolo della soluzione

ed utilizza il metodo AdvancingStep che è, nella classe Solver, pure virtual. Tale metodo

viene implementano nelle classi �glie che ereditano da Solver e de�niscono i singoli

metodi speci�ci di avanzamento in tempo elencati in precedenza. Oltre alla funzione

AdvancingStep, a seconda del tipo di risolutore utilizzato altri metodi saranno soggetti

ad overriding nelle classi �glie e altre variabili de�nite. In particolare per i metodi im-

pliciti si de�niscono le variabili di M_nmax_newton, M_toll_newton e M_DecompMet

che servono per de�nire rispettivamente il massimo numero di iterazioni, la tolleranza

e il metodo di decomposizione delle matrici per il metodo di Newton. Esse avranno un

valore di default modi�cabile attraverso le funzioni ChangeDecompMet e ChangeNew-

8

tonParam. Sempre negli impliciti si fa l'override anche della funzione NonLinearFunJac

che de�nisce a seconda del metodo lo jacobiano e la forzante del sistema lineare da ri-

solvere ad ogni passo di Newton. Viene inoltre fatto l'override della funzione DataisOK

che controlla i dati passati al problema prima di iniziare a risolverlo, per controllare

anche che le derivate siano state caricate. Per i metodi multistep (Adams-Bashforth e

Adams-Moulton) si aggiunge la variabile M_OneStepMet per il metodo ad un passo da

usare ai primi passi che di default è Eulero esplicito, ma può essere cambiato con la

funzione ChangeOneStep. Tale funzione prende valori da un enum OneStepSol che può

assumere i seguenti valori: OneStepExplicitEuler, OneStepImplicitEuler, OneStepHeun,

OneStepRK4, OneStepCrankNicolson.

L'istanziamento di oggetti delle classi concrete relative ai risolutori può certamente

avvenire direttamente dichiarando esplicitamente il nome della classe del risolutore cor-

rispondente. Viene però implementata un'altra classe speci�ca denominata Creator che

permette di istanziare un puntatore ad un oggetto del risolutore desiderato, ovviamente

sia attraverso passaggio da �le che da codice, utilizzando i metodi interni methods.

Tale classe restituisce un puntatore a Solver che ovviamente in base al tipo di risolutore

richiesto sarà un puntatore alla classe del risolutore concreto corrispondente. Inoltre tale

Creatore è un singleton, di tale determinata classe verrà creata una e una sola istanza. I

methods, oltre agli input de�niti in precedenza a seconda del tipo di passaggio dei dati,

richiede che sia speci�cato il metodo da creare de�nito come un enum SolverType che

può assumere i seguenti valori: MetExplicitEuler, MetImplicitEuler, MetCrankNicolson,

MetHeun, MetRK4, MetAdamsBashforth, MetAdamsMoulton.

Tutorial

Se si utilizza il passaggio da �le per prima cosa il �le di dati deve essere creato. Qui se

ne riporta un esempio.

# −∗− GetPot −∗−#−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−# Data f i l e f o r i n t e g r a t i o n

#−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−#Regole per l ' i n se r imento d i f o r z a n t i e s o g l i e e r i s p e t t i v e de r i v a t e

#−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−#fo r z a n t i : f0 , f 1 . . .

#de r i va t e f o r z a n t i : df0 /dx0 , df0 /dx1 . . . .

#Tutte l e f un z i on i s c r i t t e diversamente o mancanti verranno l e t t e

#dal Parser come 0 .

[ODE]

[ . / OutputFile ]

f i l ename = output

[ . . / ]

9

[ . / time ]

t_start = 0 .

t_end = 4 .

t_step = 0 .1

[ . . / ]

[ . / equat ion ]

num_eq = 2 .

CI = ' 1 . 2 . '

[ . . / ]

[ . / Imp l i c i tSo lve rParamete r s ]

nmax_newton = 40

tol l_newton = 5 . e−3

[ . . / ]

[ . / f unc t i on s ]

f 0 = t ∗ s i n ( x0 )f 1 = x0∗ cos ( t )

[ . . / ]

[ . / func t i on sDer i v ]

df0 /dx0 = t ∗ cos ( x0 )df0 /dx1 = 0 .

df1 /dx0 = cos ( t )

df1 /dx1 = 0 .

[ . . / ]

A tale punto si de�nisce la variabile GetPot per il �le, si istanzia i puntatore al risolutore

e si richiama il metodo DoIt per eseguire la simulazione

GetPot c l (" data " ) ;

s t r i n g s e c t i o n = "ODE/" ;

So lve r ∗ pluto=Creator : : GetInstance()−>method ( MetExpl ic i tEuler , c l , s e c t i o n ) ;

pluto−>DoIt ( ) ;

Ovviamente, non è questo il caso perchè si usa un Eulero esplicito, ma se si necessita di

funzioni ausiliarie prima del DoIt esse vengono richiamate. Alcuni esempi sono riportati:

pluto−>ReadForcerDeriv ( c l ) ;

pluto−>ChangeOneStep ( OneStepImpl i c i tEuler ) ;

pluto−>ChangeNewtonParam (50 , 1 . e−10);

A tale punto come post processing, se si vuole visualizzare a video i risultati e creare

l'immagine con gnuplot, basta fare:

pluto−>Pr intSo l ( ) ;system (" gnuplot GnuplotViewOde " ) ;

10

Se si utilizza il passaggio da funtori si devono per prima cosa creare le funzioni corrispon-

denti. Poi si creano i vettori e le strutture dati desiderate e si passano. Un esempio è il

seguente:

double t s t a r t =0;

double tend=4;

double t s t ep =0.1 ;

unsigned i n t numeq=2;

ublas : : vector<double> CI ( f ) ;

C[ 0 ]=1 ;

C[ 1 ]=2 ;

output f i l ename = "Prova " ;

ub las : : vector<RealFunctor∗> fo r c e r 2 ( f ) ;

f o r c e r 2 [0 ]=new RealFunctor ( fun1 , 2 ) ;

f o r c e r 2 [1 ]=new RealFunctor ( fun2 , 2 ) ;

So lve r ∗ paper ino=Creator : : GetInstance()−>method (MetAdamsBashforth ,

t s t a r t , tend , t s tep , numeq , CI , f o r c e r2 , output f i l ename ) ;

paperino−>DoIt ( ) ;

paperino−>Pr intSo l ( ) ;

2.3 Libreria per sistemi di ODE-DRH

La seconda libreria implementa il metodo descritto nella prima sezione per la soluzione

di sistemi di equazioni di�erenziali ordinarie con super�ci di discontinuità. Essa utilizza

i metodi di avanzamento in tempo della prima libreria e testa le condizioni per cambiare,

dove serve, il campo di forzanti a cui sono soggette le variabili. Non si considerano casi

in cui più soglie sono attive contemporaneamente. In �gura 2 i riporta il corrispondente

diagramma UML.

Input/Output

La libreria per sistemi di ODE DRH richiede in ingresso la �nestra temporale in cui

si vuole approssimare numericamente la soluzione del sistema (tempo iniziale e tempo

�nale), il passo temporale che si vuole utilizzare, il numero di variabili, le forzanti del

sistema sotto soglia, le forzanti del sistema sopra soglia, le derivate rispetto al tempo di

tali forzanti, il loro gradiente, le soglie, il gradiente delle soglie, le derivate rispetto al

tempo di tali soglie, le condizioni iniziali, e il nome del �le su cui si vogliono salvare i

risultati. Si deve inoltre speci�care il metodo desiderato per la risoluzione del problema

al di fuori dalle soglie e l'eventuale metodo di�erente su soglia. Per la risoluzione su

soglia si utilizzano metodi espliciti quindi se il metodo generale è esplicito viene attribuito

quello, se no viene attribuito di default Runge-Kutta del quarto ordine. Tale metodo sarà

modi�cabile mediante opportuna funzione. Nel caso di metodo generale multistep sarà

inoltre modi�cabile attraverso opportuna funzione il metodo ad un passo corrispondente.

I due metodi di passaggio per i dati del problema descritti per la libreria precedente sono

validi anche in questo caso e le funzioni sono trattate nei due modi descritti. Si noti

11

Figure 2: Diagramma UML per la libreria per sistemi di ODE-DRH.

inoltre che per il passaggio da �le, si attribuisce tramite parser di default la funzione

identicamente nulla quindi nel caso comune in cui questo sia corretto per molte delle

derivate del sistema considerato, non è necessario scrivere tali derivate sul �le (per

esempio nel caso autonomo le derivate in tempo semplicemente non si scrivono).

In outuput la libreria scrive a video gli istanti temporali, i corrispondenti valori delle

variabili e i �ag corrispondenti (se si trova la variabile sopra 1, sotto −1 o sulla soglia

0). Inoltre produce in automatico due �le, uno con la soluzione e il secondo che è il �le

da lanciare per ottenere l'immagine tramite gnuplot.

Speci�che codice

La classe centrale della libreria è OdeDrhSolver che implementa il risolutore di sistemi di

ODE DRH secondo la teoria di Dieci-Lopez estesa anche a casi non autonomi. Come nel

caso precedente si memorizza la soluzione in un ublas vector di pile. Si utilizzano inoltre

due membri M_ODESolver e M_ODESolver_on di tipo puntatore a Solver generico

che verranno utilizzati per avanzare in tempo di un passo al di fuori di soglia e con

soglia attiva rispettivamente. Ad ogni transizione si andranno a cambiare i functor

speci�ci di tali metodi attraverso il metodo pubblico ChangeForcer della classe Solver

e si selezionerà il risolutore da utilizzare. Ovviamente nel costruttore tali puntatori a

Solver generico andranno a puntare ad un oggetto della classe �glia del metodo Solver

a seconda del metodo di integrazione richiesto. Le funzioni richieste e date come input

saranno de�nite come ublas vector di puntatori a functor e trattate come nel caso della

libreria precedente a seconda di come vengono passate. In più alla libreria precedente

c'è inoltre una variabile M_time necessaria per memorizzare i tempi perchè, quando si

vanno a cercare gli istanti esatti di transizione dentro o fuori da soglie, le soluzioni non

saranno solo nei tempi da M_t_start a M_t_end con il passo de�nito ma vi saranno in

12

più anche i t∗ e le corrispondenti soluzioni a tali istanti di tempo. Inoltre vi è anche una

variabile M_�ag de�nita come un vettore di pile di interi che serve per memorizzare le

�ag delle soluzioni ai vari istanti (−1 se sotto soglia, 0 se su soglia, 1 se sopra soglia)

per ogni variabile. Per ogni variabile si sceglie di usare una pila in quanto il valore di

interesse per calcolare la soluzione successiva sarà sempre e solo l'ultimo, dove si trova

la soluzione al passo precedente.

Il metodo fondamentale per l'integrazione del sistema di ODE DRH è DoIt che

de�nisce la struttura generale dell'algoritmo e si occupa di creare e riempire i �le (sia

per gnuplot che per le variabili di soluzione). In esso si utilizza il metodo AdvancingStep

che e�ettivamente implementa la teoria per la risoluzione di tale tipo di sistema. Ad ogni

passo si attua l'avanzamento in tempo di uno step con il risolutore attivo a seconda delle

ultime �ag caricate. Si testa poi con una funzione WhereAmI la soluzione di arrivo. Se

eravamo inzialmente fuori da soglie e dopo l'avanzamento in tempo la situazione non è

cambiata non è necessario testare nulla e AdvancingStep restituisce il risultato ottenuto.

Se si attraversa una soglia si ricerca con ZeroGFind il tempo di arrivo sulla soglia e la

soluzione a tale tempo che vengono caricate nelle strutture dati opportune. A tale punto

si testano le condizioni e, in base alla teoria, si deduce come la soluzione deve proseguire

e si porta a compimento il passo temporale muovendosi con tale nuovo campo di forzanti

del tempo mancante. Se si è su soglia ad ogni passo si testano le condizioni della teoria

per capire se si deve uscire. Se si rimane su soglia si utilizza project_on per riproiettare

la soluzione numerica sulla soglia attiva. Se si deve uscire si calcola con ZeroGammaFind

il tempo in cui si deve uscire e si conclude il passo temporale con il nuovo campo. Se ad

un avanzamento in tempo si veri�ca una condizione non de�nita dalla teoria o si veri�ca

una condizione instabile e con soluzione non unica si restituisce un errore corrispondente

e si arresta la simulazione. In ogni caso in cui vi sono cambiamenti nei �ag M_usedforcer

permette di calcolare le nuove forzanti da usare.

Le funzioni A, B e C servono per calcolare i corrispondenti coe�cienti da testare

per la teoria di Dieci-Lopez (come de�nite in (3), (6) e (7)). Si noti come, se per

le A le funzioni passate esatte sono su�cienti per il calcolo, per il calcolo delle B e

C (± e 1, 2) servirebbero molte più derivate (derivate rispetto alle variabili di stato

�no al terzo ordine per le funzioni di soglia, derivate �no al secondo ordine per le

forzanti, derivate temporali �no al secondo ordine per le forzanti e derivate temporali

�no al secondo ordine per i gradienti di tutte le funzioni di soglia). Si è scelto di non

passare tali dati come funzioni esatte al costruttore ma di approssimarle con i dati che

si hanno attraverso un metodo del quarto ordine a di�erenze �nite implementato in

DFGradGamma, DFDgammaDt, DFdFidt (rispettivamente per il gradiente di gamma,

per le derivate temporali di soglie e le derivate temporali delle forzanti su soglia) e con

di�erenze �nite per derivate seconde per l'hessiano di gamma in DFHessGamma. La

scelta di approssimare numericamente tali derivate è fatta nella logica di rendere il codice

più utilizzabile da un utente esterno che non deve passare un numero di dati che in caso

13

contrario risulterebbe molto elevato. Inoltre si sottolinea che tale approssimazione con

tecniche di ordine elevato sarà accurata e che viene fatta unicamente se le condizioni del

secondo ordine sono necessarie (se capita che qualche A risulti nulla). Si aggiunge che

tali funzioni di approssimazione vengono utilizzate unicamente per il calcolo di B e C

di cui per l'analisi ci interessa unicamente il segno o il valore nullo.

Sono fornite inoltre funzioni ausiliarie come ChangeODEOneStep per cambiare nel

caso di metodo multistep il metodo ad un passo usato ai primi passi e ChangeNewtonPar

per cambiare i parametri di Newton nel caso di metodi implicito che si rifanno a quelli

della classe Solver della libreria precedente. Per un trattamento opportuno dei metodi

multistep si sottolinea che ad ogni cambiamento di forzanti da usare derivante da un

cambiamento della condizione della soluzione rispetto alle soglie, le variabili nel risolutore

Solver vengono ad essere resettate all'unico ultimo valore raggiunto.

Nel caso di soluzione su soglia attiva le forzanti corrispondenti devono essere de�nite

come combinazione lineare dei campi coinvolti come de�nito in (4). Per poterle trattare

anch'esse come functor come descritto precedentemente si è implementata una classe che

deriva da functor (come RealFunctor e ParsedFunctor) che si chiama ThresholdFunctor.

Essa carica i dati necessari e fa l'override del GetValue di functor in modo tale da de�nire

la alpha in modo corretto ed eseguire la combinazione lineare richiesta. Si noti come

il motivo per cui si usano unicamente solutori espliciti per casi in cui ci si muove su

soglia attiva è che in caso contrario servirebbero i gradienti di tali ThresholdFunctor da

passare al risolutore. Questo vorrebbe dire o passare le derivate corrispondenti da �le

che sarebbero però in numero molto elevato e non idoneo all'utilizzabilità del codice o

approssimarle numericamente. L'ipotesi di approssimarle numericamente è stata esclusa

in quanto introdurre un errore numerico sulla valutazione di tali derivate da usare nel

metodo di Newton per l'avanzamento in tempo avrebbe peggiorato notevolmente le

prestazioni del metodo numerico ed avrebbe reso vano l'utilizzo di un metodo implicito

per diminuire l'errore numerico derivante dall'integrazione.

Tutorial

Se si utilizza il passaggio da �le per prima cosa il �le di dati deve essere creato. Qui

se ne riporta un esempio.

# −∗− GetPot −∗−#−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−# Data f i l e f o r i n t e g r a t i o n

#−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−#Regole per l ' i n se r imento d i f o r z a n t i e s o g l i e e r i s p e t t i v e de r i v a t e

#−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−#Sog l i e : g0 , g1 , g2 . . . .

#Gradiente d e l l e s o g l i e : dg0/dx0 , dg0/dx1 , dg0/dx2 . . . .

# dg1/dx0 , dg1/dx1 , dg1/dx2 . . .

#Derivata in tempo d e l l e s o g l i e : dg0/dt , dg1/dt , . . .

#Forzant i s o t t o s o g l i a : f0down , f1down , f2down . . .

14

#Gradiente d e l l e f o r z a n t i s o t t o s o g l i a : df0down/dx0 , df0down/dx1 , . . . .

# df1down/dx0 , df1down/dx1 , . . . .

#Der ivate in tempo d e l l e f o r z a n t i s o t t o s o g l i a : df0down/dt , df1down/dt . . . .

#Le f o r z a n t i sopra s o g l i a s i s c r i vono n e l l a medesima maniera so s t i tu endo

#down con up .

#Tutte l e f un z i on i s c r i t t e diversamente o mancanti verranno l e t t e

#dal Parser come 0 .

[ODE]

[ . / OutputFile ]

f i l ename = PROVASCOTTICOMPLETO

[ . . / ]

[ . / time ]

t_start = 0 .

t_end = 10 .

t_step = 0.01

[ . . / ]

[ . / equat ion ]

num_eq = 3 .

CI = ' 1 . 0 . 0 . '

[ . . / ]

[ . / Imp l i c i tSo lve rParamete r s ]

nmax_newton = 40

tol l_newton = 5 . e−3

[ . . / ]


f0down = −2.∗(x0^2)f0up = −2.2∗( x0^2)f1down = −x1+x1^2f1up = −0.2∗( x0^2)f2down = ( x0^2)−2.∗x1+(x1^2)f2up = −1.∗(x0^2)

[ . . / ]

[ . / df_dt ]

[ . . / ]

[ . / GradF ]

df0down/dx0 = −4.∗x0df0up/dx0 = −4.4∗x0df1down/dx1 = −1.+2.∗x1df1up/dx0 = −0.4∗x0

15

df2down/dx0 = 2 .∗ x0df2down/dx1 = −2.+2.∗x1df2up/dx0 = −2.∗x0

[ . . / ]

[ . / th r e sho ld ]

g0 = −0.9∗x0g1 = x1−0.1∗x0g2 = x2−0.5∗x0[ . . / ]

[ . / GradG ]

dg0/dx0 = −0.9

dg1/dx0 = −0.1

dg1/dx1 = 1 .

dg2/dx0 = −0.5

dg2/dx2 = 1 .

[ . . / ]

[ . / dG_dt ]

[ . . / ]

[ . . / ]

A tale punto si de�nisce il GetPot per il �le, si crea l'oggetto della classe e si lancia la

simulazione. I risultati verranno stampati a video e salvati in un �le. Si usa inoltre il

�le per gnuplot salvato in automatico dal codice per generare il gra�co delle soluzioni

in funzione del tempo.

GetPot c l (" data " ) ;

s t r i n g s e c t i o n = "ODE/" ;

OdeDrhSolver pippo (MetRK4, c l , s e c t i o n ) ;

pippo . DoIt ( ) ;

system (" gnuplot GnuplotView " ) ;

Come nel caso della prima libreria, se si vogliono utilizzare funzioni ausiliarie per cam-

biare i parametri di Newton per esempio si usa la funzione corrispondente prima di

lanciare la simulazione. Nel caso si usi un passaggio delle forzanti da RealFunctor, essi

devono essere creati e passati da codice:

ublas : : vector<double> CI ( 3 ) ;

CI [ 0 ]= 1 . ; CI [ 1 ]= 0 . ; CI [ 2 ]= 0 . ;

ub las : : vector<RealFunctor∗> ForcerDown ( 3 ) ;

ForcerDown [0 ]=new RealFunctor ( f0down , 1 ) ;


16


ub las : : vector<RealFunctor∗> ForcerUp ( 3 ) ;

ForcerUp [0 ]=new RealFunctor ( f0up , 1 ) ;



ub las : : vector<RealFunctor∗> GradForcerDown ( 9 ) ;

GradForcerDown [0 ]=new RealFunctor ( df0downdx0 , 1 ) ;


GradForcerDown [2 ]=new RealFunctor ( zero , 1 ) ;







ub las : : vector<RealFunctor∗> GradForcerUp ( 9 ) ;

GradForcerUp [0 ]=new RealFunctor ( df0updx0 , 1 ) ;


GradForcerUp [2 ]=new RealFunctor ( zero , 1 ) ;







ub las : : vector<RealFunctor∗> dForcerDowndt ( 3 ) ;

dForcerDowndt [0 ]=new RealFunctor ( zero , 1 ) ;



ub las : : vector<RealFunctor∗> dForcerUpdt ( 3 ) ;

dForcerUpdt [0 ]=new RealFunctor ( zero , 1 ) ;



ub las : : vector<RealFunctor∗> g ( 3 ) ;

g [0 ]=new RealFunctor ( g0 , 1 ) ;



ub las : : vector<RealFunctor∗> GradG ( 9 ) ;

GradG[0]=new RealFunctor ( dg0dx0 , 1 ) ;

GradG[1]=new RealFunctor ( zero , 1 ) ;








ub las : : vector<RealFunctor∗> dGdt ( 3 ) ;

dGdt [0 ]=new RealFunctor ( zero , 1 ) ;



17

s t r i n g output f i l ename = "Prova " ;

unsigned i n t num_eq=3;

OdeDrhSolver paper ino (MetAdamsMoulton , 0 . , 1 0 . , 0 . 0 1 , num_eq , CI ,

ForcerDown , ForcerUp , GradForcerDown ,

GradForcerUp , dForcerDowndt , dForcerUpdt , g ,GradG , dGdt , output f i l ename ) ;

paper ino . DoIt ( ) ;

system (" gnuplot GnuplotView " ) ;

dove ovviamente le funzioni utilizzate per de�nire i RealFunctor devono essere create

precedentemente e incluse.

18

3 Casi test

3.1 Libreria ODE

Si mostrano due risultati relativi alla libreria di ODE per sistemi di ODE senza discon-

tinuità per tre variabili di stato. Il primo esempio è stato preso dalla tesi della Dr.ssa

Scotti [3], il secondo è un caso utilizzato in quanto si aveva a disposizione la soluzione

esatta per veri�care la veridicità della simulazione. Nel primo caso le forzanti erano le

seguenti:

# SENZA SOGLIE

f0 = −2.∗(x0^2)+(x1^2)

f1 = ( x0^2)−x1−2.∗( x1^2)f2 = ( x0^2)+x1+2.∗( x1^2)

Nel secondo caso si usano le seguenti forzanti:

# SISTEMA ODE 2

f0down = −x0−1.+t+t^2

f1down = x0−(x1 )^2.+( exp (−2.∗ t ))∗((− exp (2∗ t ))∗(−1+ t )∗ t++ (1−2∗ t+t ^2)∗ t^2−(exp ( t ))∗(1−3∗ t+t ^2))

f2down = x1^2.−( exp(−2∗ t ))∗(1−2∗ t+t ^2)∗ t^2+(t−1)∗ cos ( t)+ s i n ( t )

0 0.2 0.4 0.6 0.8 1−0.25

−0.2

−0.15

−0.1

−0.05

0

t

x0,x

1,x

2

SISTEMA ODE 2

x0

x1

x2

0 1 2 3 4 5 60

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

t

x0,x

1,x

2

SENZA SOGLIE

x0

x1

x2

Figure 3: Risultati relativi all'esempio per ODE non vincolate della Scotti (sx) e della soluzione esatta

(dx).

Per il caso del sistema ODE 2 avendo le soluzioni esatte si è potuto veri�care che

gli errori fossero in accordo con la teoria relativa a ciascun metodo di integrazione. In

particolare per la prima variabile di stato, essendo quadratica e indipendente dalle altre,

i metodi del secondo ordine risultano esatti.

19

3.2 Libreria ODEDRH-Autonomi

Si mostrano alcuni risultati ottenuti con la libreria di ODE DRH nel caso con funzioni

autonome sia per la soglia che per le forzanti. I primi due esempi sono a due variabili

di stato con soglie costanti. Le corrispondenti forzanti e soglie sono riportate in seguito.

Nelle �gure sono inoltre mostrati dei dettagli in cui si vede al passaggio della soglia

dove si andrebbe proseguendo con il campo precedente per tutto il passo e dove si va

controllando le condizioni e prendendo la direzione giusta. Si può vedere così l'e�etto e

l'importanza del calcolo del t∗.

1 2 3 4 5 61

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6

x0

x1

AUTONOMO 1

4 4.005 4.01 4.015 4.02

4

4.002

4.004

4.006

4.008

4.01

4.012

4.014

4.016

x0

x1

AUTONOMO 1−DET

5 5.002 5.004 5.006 5.008 5.01 5.012 5.014

4.8

4.802

4.804

4.806

4.808

4.81

4.812

4.814

4.816

4.818

x0

x1

AUTONOMO 1−DET

Figure 4: Risultati relative all'esempio autonomo 1 con dettagli sull'attraversamento di ciascuna

soglia.

# AUTONOMO 1


f0down = 1

f0up = 0 .5

f1down = 1

f1up = 0 .8

[ . . / ]


g0 = x0−5.005.

g1 = x1−4.005.

[ . . / ]

20

1 2 3 4 5 6 71

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6

x0

x1

AUTONOMO 2

4 4.005 4.01 4.015 4.02

4

4.002

4.004

4.006

4.008

4.01

x0

x1

AUTONOMO 2−DET

Figure 5: Risultati relative all'esempio autonomo 2 con dettagli sulla caduta sulla soglia.

#AUTONOMO 2


f0down = 1

f0up = 0 .5

f1down = 1

f1up = −0.8

[ . . / ]


g0 = x0−5.005.

g1 = x1−4.005.

[ . . / ]

Nell'ultimo esempio autonomo si mostra un caso preso dall'articolo di Dieci-Lopez

con una unica soglia (o due soglie con una non raggiungibile) e due variabili di stato.

La soluzione dopo un transitorio iniziale entra in un ciclo come mostrato in �gura. Si

mostrano in tale caso anche esempi di soluzioni ottenibili con la libreria di ODE utilizzata

per lo stesso problema con campo discontinuo.

# ESEMPIO DIECI


f0down = x1

f0up = x1

f1down = −x0+1/(1.2−x1 )

f1up = −x0−1/(0.8+x1 )

[ . . / ]


g0 = 2 .

g1 = x1−0.2

[ . . / ]

21

−1.5 −1 −0.5 0 0.5 1 1.5−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

x0

x1

PROVA DIECI

−1.13 −1.12 −1.11 −1.1 −1.09 −1.08 −1.07 −1.06 −1.05

0.05

0.1

0.15

0.2

0.25

0.3

x0

x1

PROVA DIECI

−0.1 −0.09 −0.08 −0.07 −0.06 −0.05 −0.04

0.18

0.2

0.22

0.24

0.26

0.28

0.3

x0

x1

PROVA DIECI

−0.95 −0.945 −0.94 −0.935 −0.93 −0.925 −0.92

0.194

0.196

0.198

0.2

0.202

0.204

0.206

x0

x1

PROVA DIECI CONF

0.96 0.98 1 1.02 1.04

0.188

0.19

0.192

0.194

0.196

0.198

0.2

0.202

0.204

x0

x1

PROVA DIECI CONF

Figure 6: Risultati relative all'esempio del paper di Dieci-Lopez. Dettagli all'incontro con le soglie e

confronto con la soluzione ottenibile la libreria di ODE.

22

3.3 Libreria ODEDRH-Non Autonomi

Si mostrano ora esempi non autonomi. I primi tre esempi hanno la non autonomia

unicamente nelle forzanti. Il quarto avrà la non autonomia sia nelle forzanti che in

una soglia allo scopo di testare anche tale aspetto del codice. Anche in questo caso si

riportano gli esempi di confronto con la soluzione ottenuta con la libreria di ODE. In tal

caso per la libreria di ODE si usa un passo un decimo di quello utilizzato con la libreria

di ODE discontinue. Si riportano forzanti e soglie corrispondenti.

#NON AUTONOMO 1


f0down = 1

f0up = 0 .5

f1down = 1

f1up = −0.8+0.1776∗ t[ . . / ]


g0 = x0−5.005.

g1 = x1−4.005.

[ . . / ]

0 2 4 6 8 10 121

2

3

4

5

6

7

x0

x1

NON AUTONOMO 1

5.5 5.52 5.54 5.56 5.58 5.6 5.62 5.64 5.66

4.004

4.005

4.006

4.007

4.008

4.009

x0

x1

NON AUTONOMO 1−DET

5.15 5.2 5.25 5.3 5.35 5.4

3.998

4

4.002

4.004

4.006

4.008

4.01

4.012

4.014

4.016

x0

x1

NON AUTONOMO 1−CONF

Figure 7: Risultati relative all'esempio non autonomo 1 con dettagli sull'uscita dalla soglia e confronto

con la soluzione ottenibile con libreria per ODE classiche.

23

#NON AUTONOMO 2


f0down = 1

f0up = −0.5

f1down = 1

f1up = −0.8+0.06512∗( t ^2)[ . . / ]


g0 = x0−5.005.

g1 = x1−4.005.

[ . . / ]

1 2 3 4 5 6 71

2

3

4

5

6

7

8

x0

x1

NON AUTONOMO 2

5 5.002 5.004 5.006 5.008 5.014.061

4.062

4.063

4.064

4.065

4.066

4.067

4.068

4.069

4.07

x0

x1


4.5 4.52 4.54 4.56 4.58 4.6 4.62

3.995

4

4.005

4.01

4.015

4.02

4.025

4.03

4.035

x0

x1

NON AUTONOMO 2−DEF

5.0035 5.004 5.0045 5.005 5.0055 5.006 5.0065 5.0074.055

4.06

4.065

4.07

4.075

4.08

x0

x1




#NON AUTONOMO 3


f0down = 1

f0up = −0.5+0.099∗ tf1down = 1

f1up = −0.8+0.06512∗( t ^2)[ . . / ]


24

g0 = x0−5.005.

g1 = x1−4.005.

[ . . / ]

1 2 3 4 5 6 70

2

4

6

8

10

12

14

16

18

20

x0

x1

NON AUTONOMO 3

5.0049 5.005 5.0051 5.0052 5.0053 5.0054

4.56

4.58

4.6

4.62

4.64

4.66

4.68

4.7

4.72

4.74

x0

x1


5.003 5.004 5.005 5.006 5.007 5.008 5.009

4.3

4.4

4.5

4.6

4.7

4.8

4.9

5

x0

x1




#NON AUTONOMO 4


f0down = 1

f0up = 0 .6

f1down = 1

f1up = −0.8

[ . . / ]


g0 = x0−5.005.

g1 = x1−2.003−0.5∗ t .[ . . / ]

25

2 3 4 5 6 7 8 9

1

2

3

4

5

6

7

NON AUTONOMO 4

x1

x0

3 3.005 3.01 3.015 3.02

3

3.005

3.01

3.015

x0

x1


4.995 5 5.005 5.01 5.015

4

4.002

4.004

4.006

4.008

4.01

4.012


x0

x1



3.4 Libreria ODEDRH-Ordine superiore

Si riportano in�ne dei risultati che coinvolgono l'utilizzo di condizioni di ordine superiore.

Tali esempi corrispondono al medesimo caso a tre variabili di stato riportato per il caso

test per la libreria di ODE ma con soglie. Inoltre il primo esempio è autonomo, il secondo

ha non autonomia nelle soglie del problema.

#SECONDO ORDINE AUTONOMO


f0down = −2.∗(x0^2)+(x1^2)

f0up = −2.2∗( x0^2)+1.1∗( x1^2)f1down = ( x0^2)−x1−2.∗( x1^2)f1up = −0.2∗( x0^2)+0.1∗( x1^2)f2down = ( x0^2)+x1+2.∗( x1^2)f2up = −1.∗(x0^2)+0.5∗( x1^2)[ . . / ]


26

g0 = −0.9∗x0g1 = x1−0.1∗x0g2 = x2−0.5∗x0[ . . / ]

#SECONDO ORDINE NON AUTONOMO


f0down = (−2.∗( x0^2)+(x1^2))∗ l og ( t+2)f0up = (1.1+0.11∗ t )∗ ( ( −2 .∗ ( x0^2)+(x1^2))∗ l og ( t+2))+x0 ∗0 .11f1down = ( ( x0^2)−x1−2.∗( x1^2))∗ l og ( t+2)f1up = (0.1+0.01∗ t )∗ ( ( −2 .∗ ( x0^2)+(x1^2))∗ l og ( t+2))+x0 ∗0 .01f2down = ( ( x0^2)+x1+2.∗( x1^2))∗ l og ( t+2)f2up = (0.5+0.05∗ t )∗ ( ( −2 .∗ ( x0^2)+(x1^2))∗ l og ( t+2))+x0 ∗0 .05[ . . / ]


g0 = x0−(1.+ t /10 . )∗1 . 1∗ x0g1 = x1−(1.+ t /10 . )∗0 . 1∗ x0g2 = x2−(1.+ t /10 . )∗0 . 5∗ x0[ . . / ]

0 2 4 6 8 10 120

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

t

x0,x

1,x

2

SECONDO ORDINE − AUTONOMO

x0

x1

x2

0.09 0.1 0.11 0.12 0.13

0.07

0.075

0.08

0.085

0.09

0.095

t

x0,x

1,x

2

SECONDO ORDINE −AUTONOMO−DET

x0

x1

x2

0 2 4 6 8 10 120

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

t

x0,x

1,x

2

SECONDO ORDINE− NON AUTONOMO

x0

x1

x2

Figure 11: Risultati relative ai casi con condizione di ordine superiore con dettaglo sull'ingresso alla

soglia e caso non autonomo.

27

Bibliography

[1] Filippov A.F. Di�erential equations with discontinuous right hand sides. Mathemat-

ics and Its applications, 1988.

[2] Lopez L.. Dieci L. Sliding motion in discontinuous di�erential systems: Theory and

a computational apporach. Georgia Institute of Technology, 2008.

[3] Scotti A. A numerical model for generation, retention and expulsion of hydrocarbons

from source rocks. PhD thesis, Politecnico di Milano, 2010.

28

LIBRERIA C++ PER L'INTEGRAZIONE DI SISTEMI DI ODE...

Documents

Transcript of LIBRERIA C++ PER L'INTEGRAZIONE DI SISTEMI DI ODE...