Libreria C++ per la riduzione di modello di problemi ...forma/Didattica/ProgettiPacs/... ·...

Libreria C++ per la riduzione di modello di problemi

simplettici

Nicolo Ripamonti

Milano, Italia

Abstract

Questa libreria ha lo scopo di fornire un toolbox completo per la riduzionedi modello riguardante sistemi simplettici, lineari e non. Una delle questionicritiche riguardo la riduzione di modello e la conservazione della strutturasimplettica anche per il modello ridotto: la conservazione di questa proprietapermette di formulare stime di stabilita, basate sulla definizione di Lyapunov,garantendo significativi vantaggi in termini di stabilita per simulazioni dilunga durata. Il metodo implementato prende il nome di Cotangent Lift, unmetodo di riduzione a basi ridotte, e una variante del DEIM che preserva lastruttura simplettica del problema viene presa in considerazione.

Keywords: MOR, Proper Orthogonal Decomposition (POD), HyperbolicProblems, Hamiltonian Systems, Symplectic integration

November 28, 2017

1

Contents

1 Introduzione 3

2 Model Order Reduction 4

3 Proper Orthogonal Decomposition (POD) 63.1 Trattamento termine nonlineare . . . . . . . . . . . . . . . . . 9

4 Proper Symplectic Decomposition (PSD) 114.1 Riduzione di modello simplettica . . . . . . . . . . . . . . . . 144.2 Trattamento termine nonlineare . . . . . . . . . . . . . . . . . 18

5 Implementazione 195.1 Compilare, installare ed usare la libreria . . . . . . . . . . . . 195.2 Estensione delle funzionalita di Eigen . . . . . . . . . . . . . . 205.3 Struttura della libreria . . . . . . . . . . . . . . . . . . . . . . 245.4 Sistema dinamico . . . . . . . . . . . . . . . . . . . . . . . . . 25

5.4.1 Problemi trattati e differenze finite . . . . . . . . . . . 285.5 ODE solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

5.5.1 Funzione stepper . . . . . . . . . . . . . . . . . . . . . 315.5.2 Funzione di integrazione e observer pattern . . . . . . . 35

5.6 Riduzione di modello . . . . . . . . . . . . . . . . . . . . . . . 39

6 Test numerici 416.1 Equazione lineare delle onde . . . . . . . . . . . . . . . . . . . 426.2 Equazione di Sine-Gordon . . . . . . . . . . . . . . . . . . . . 44

7 Conclusioni e futuri miglioramenti 48

2

1. Introduzione

La maggior parte dei modelli fisici sono espressi tramite equazioni aderivate parziali (PDE) multiparametriche: questo tipo di descrizione deiproblemi e comune a diverse branche della scienza e dell’ingegneria, daltrasferimento stazionario o turbolento di massa e calore, all’applicazione inambito di meccanica dei fluidi o acustica e a problemi di finanza. In questicasi i parametri in ingresso possono caratterizzare il problema, distinguen-dolo dagli altri, in termini di coefficienti fisici, condizioni al contorno, datiiniziali o domini sotto studio. Di frequente, piu che alla soluzione numerica inquanto tale si e interessati a funzioni di essa, come la massima temperaturaregistrabile in un dato ambiente, il coefficiente di drag di un dato veicolo o lamassima caduta di pressione all’interno di un condotto: tipici casi a questoproposito sono dati dai processi di ottimizzazione e design industriali o dallavalutazione di rischio, accumunati dal dover risolvere un gran numero di voltelo stesso problema, variando i parametri che lo caratterizzano. In questi casinon solo l’accuratezza e di vitale importanza ma rilevante e soprattutto lacapacita di fornire una risposta veloce al cambiamento dei parametri, ideal-mente real-time. Per quanto efficienti implementazioni di metodi numerici eavanzamenti in ambito hardware abbiano garantito grandi risultati da questopunto di vista, il dover risolvere ripetutamente problemi computazionalmenteonerosi non permette di sfruttare un approccio diretto. Andremo quindi adescrivere un approccio a questo tipo di problemi basato sulla costruzione diun sistema di basi ad-hoc per ottenere sistemi ridotti e piu veloci da risol-vere senza sacrificare l’accuratezza, piu comunemente noto come riduzione dimodello tramite basi ridotte. L’idea non e quella di sostituire metodi direttiper la risoluzione di PDE ma piuttosto di utilizzare metodi che sfruttanoun certo numero di soluzioni gia calcolate tramite metodi onerosi. A com-plicare ulteriormente lo scenario considerato, sottolineamo come per certeclassi di problemi, principalmente iperbolici e caratterizzata dall’assenza dicontributi dissipativi, la stabilita dei sistemi ottenuti non e sempre garantita.Di questo gruppo di PDE fanno parte i sistemi simplettici, caratterizzati dauna paricolare struttura dell’equazione e dalla conservazione della quantitacomunemente nota come Hamiltoniano.In questo lavoro, una libreria che implementa metodi basati sulla decom-posizione a valori singolari viene considerata. In particolar metodo, e stataimplementata la variante Cotangent Lift del metodo POD che, sotto maggiorivincoli sulla forma della matrice di base, permette di conservare la delicata

3

struttura del problema. Per quest’ultimo algoritmo, anche una estensione alcaso di problemi nonaffini e non lineari viene considerata.Alcuni esempi, lineari e non, sono usati per validare le tecniche descritte.Gli esperimenti numerici considerati riguardano sia riduzioni in cui l’unicoparametro preso in considerazione e il tempo, sia problemi piu complessiaventi un parametro addizionale da tenere in conto nel processo di riduzione.Questo report e fondamentalmente diviso in due parti. La prima dedicataprincipalmente alle motivazioni per cui la riduzione di modello e vantaggiosain certi contesti e alla descrizione delle tecniche numeriche utilizzate. La sec-onda si propone come una descrizione delle scelte implementative fatte perla libreria in questione.

2. Model Order Reduction

Consideriamo un modello dinamico finito dimensionale dipendente dalparametro ω. La forma piu comune attraverso cui viene espresso e quella disistema di equazione differenziali ordinarie (ODEs) del primo ordine

ddt

x (t, ω) = f (x, t, ω) ,

x (0, ω) = x0.(1)

La natura del parametro puo essere molto varia, potendo esso influire sul datoiniziale x0, sulla geometria del sistema o sull’equazione f (x, t, ω). Alcuniesempi per quest’ultimo caso, tratti dal mondo della fluidodinamica, sonoil numero di Reynolds ed il numero di Prandtl. In (1), x ∈ Rn rapresentalo stato del sistema, ω ∈ Γ e un vettore contenente tutti i parametri diinteresse appartenente all’insieme compatto Γ ∈ Rd e f : Rn × R × Γ 7→Rn e una generica funzione dello stato del sistema, anch’essa dipendentedal vettore dei parametri. Un sistema del tipo (1) puo essere il risultatodella semidiscretizzazione di una equazione a derivate parziali. Centrale nelladefinizione della teoria della riduzione di modello e lo spazio contenente tuttele soluzioni di (1) al variare del vettore dei parametri ω. Esso viene definitocome

M = x (t, ω) |ω ∈ Γ, t ≥ 0 ⊂ Rn. (2)

Tranne in rari e fortunati casi, un’espressione esplicita della soluzione di (1)non e disponibile e quindi e necessario un integratore numerico per ottenereuna soluzione approssimante. Lo spazio approssimanteMδ, in cui sono con-tenute le soluzioni approssimanti al variare del parametro, in principio non

4

coincide con M. Un’ipotesi ragionevole e che l’errore di integrazione di (1)possa essere reso arbitrariamente piccolo e che quindi Mδ sia una ragionev-ole approssimazione di M. Questo ci permette di confondere i due spazi elavorare sullo spazio che possiamo avere a disposizione, ovveroMδ, indican-dolo conM. L’assunzione fondamentale riguardante la riduzione di modellotramite basi ridotte, argomento che trattiamo in questo report, e che la di-mensione di M sia di molto inferiore a n, ovvero che lo spazio generato daun ristretto numero di funzioni di base approssimi lo spazioM con un errorepiccolo. Queste basi verrano chiamate da questo momento in poi basi ridottee indicate con viki=1. Lo spazio da esse generate prende il nome di spazioridotto Mrb. E possibile defire l’approssimabilita di un problema attraversola k-dimensione di Kolmogorov [1]. Definiamo prima

E (M,Mrb) = supu∈M

infurb∈Mrb

‖u− urb‖2 . (3)

La k-dimensione di Kolmogorov di Mrb rispetto a M e quindi data da

dk (M) = infMrb

E (M,Mrb) , (4)

dove il minimo e considerato tra tutti i sottospazi k-dimensionali di M. Laquantita (4) viene quindi usata per indicare quanto bene M possa essereapprossimato da un sottospazio di dimensione k. Se dk decresce molto rap-idamente con k significa che lo spazio puo essere approssimato in manieraefficiente da un numero ristretto di funzioni di base, dove per i casi migliori ladecrescita puo essere anche esponenziale. I vari algoritmi che descriviamo neirispettivi capitoli si pongono principalmente l’obiettivo di trovare lo spazioMrb che minimizza E per un dato numero di basi, in particolare trovando unsistema di basi ridotte dello spazio stesso.

Supponendo di avere a disposizione un set di basi ridotte di dimensionek, con k n, la soluzione nello spazio ridotto viene approssimata come

x ≈ Vy, (5)

dove V ∈ Rn×k e la matrice avente come colonne le basi ridotte e y ∈ Rk

contiene le coordinate di x rispetto alla base ridotta. Sostituendo (5) in (1)otteniamo il seguente sistema sovradeterminato

Vd

dty = f (t,Vy, ω) . (6)

5

Consideriamo quindi una proiezione di Galerkin [2] e, richiedendo che V siainvertibile secondo Moore-Penrose, abbiamo

d

dtx = (V)† f (t,Vy, ω) , (7)

il quale rappresenta un sistema di k equazioni in k incognite, ovvero i coef-ficienti della soluzione completa rispetto alla base ridotta. Nel caso in cui lamatrice V risulti ortogonale abbiamo in particolare che

V† = VT . (8)

La soluzione di questo sistema, se k risulta essere molto minore di n, presentaun costo significativamente minore di quella del sistema originario. Dal mo-mento che la presenza di termini non lineari nell’aggregato f potrebbe portaread un costo computazionale ancora dipendente da n per il modello ridottoridotto, riscriviamo (7) come

d

dtx = (V)† fnl (t,Vy, ω) + (V)†A (t, ω) Vy, (9)

in maniera tale da evidenziare entrambi i contributi, lineare e non, i qualiverrano trattamenti in maniera differente.

3. Proper Orthogonal Decomposition (POD)

La tecnica POD [3][4][5], conosciuta anche come decomposizione di Karhunen-Loeve o analisi delle componenti principali (PCA) ha come idea centralequella di fornire una base ortonormale che permetta di rappresentare i datiraccolti in modo ottimale in ottica di minimi quadrati. Un altro grande van-taggio di questa tecnica e che, basandosi su di una decomposizione modalecompletamente dipendente dai dati raccolti, non richiede alcuna conoscenza apriori del processo che ha generato i dati. Questo e utile quando le conoscenzefisiche del fenomeno sono minime e non e possibile identificare un criterio chepermetta di trovare delle basi per lo spazio ridotto. Inoltre, sempre in uncontesto di scarsa conoscenza del fenomeno che si sta studiando, permette difornire alcune informazioni fondamentali sul processo stesso: ad esempio, nelcaso dello studio delle equazioni di Navier Stokes, ha permesso di avere unaconoscenza piu approfondita di alcune strutture molto localizzate in spazio etempo che caratterizzazione la transizione da un regime all’altro di turbolenza

6

[6]. Questa tecnica, combinata con la tecnica di proiezione di Galerkin,precedentemente menzionata, e stato uno dei primi metodi introdotti perla riduzione di sistemi dinamidici di dimensione potenzialmente infinita asistemi molto piu piccoli e quindi facilmente integrabili. Il fatto che questometodo ricerchi fondamentalemnte spazi ridotti lineari, o affini, anziche sot-tospazi curvi dello spazio di partenza lo rende estremamente efficiente da unpunto di vista meramente implementativo. Le sue applicazioni sono moltopiu generali della riduzione di modello per sistemi dinamici: molteplici sonogli esempi di suoi utilizzi per analisi di dati e compressione di dati in variebranche dell’ingegneria e della scienza piu in generale. Applicazioni del PODincludono image processing [7], modellizzazione e controllo di sistemi chimici[8][9][10], modelli di turbolenza [10], controllo di fluidi ed ingegneria del ventoper nominarne alcuni. La formulazione standard, che tratteremo in questoreport, e stata piu volte modificata per adattamenti all’applicazione presain considerazione. Uno dei motivi per l’adattamento del metodo e quellodi conservare alcune strutture note del sistema fisico o, piu in generale, disfruttare alcune conoscenza a priori del problema.

Siano x (ti, ωj) con i = 1, . . . ,m e j = 1, . . . , n un campionamento finito,chiamato piu comunemente insieme di snapshots, dello spazio M.Se assum-iamo di disporre di una base ridotta, possiamo generare la matrice V. Eccoquindi che l’operatore di proiezione da Rn sullo spazio ridotto viene constru-ito come VVT . La base ricercata tramite il metodo POD, di dimensione k,rappresenta la soluzione del problema di ottimizzazione

minV∈Rn×k

∥∥S−VVTS∥∥F

t.c.VTV = Ik, (10)

dove S e la matrice degli snapshots contenente gli snapshots x (ti, ωj) comesue colonne. ‖·‖F e la norma di Frobenius e Ik e la matrice identita di di-mensione k. La matrice degli snapshots S ha un numero di colonne s parial numero di snapshots considerati ed un numero di righe n uguale alla di-mensione dello spazio completo. Mentre la seconda quantita risulta esserefissata, la prima dipende dal problema che si sta analizzando. Geometri-camente sappiamo che s ≤ n, perche non possiamo sperare di ottenere piudi n vettori linearmente indipendenti. Non esiste un valore ottimale per p,ma la raccolta di snapshots andrebbe in ogni caso concentrata negli istantitemporali caratterizzati da un cambio nelle dinamiche del sistema [11].

Il perche la soluzione di (10) rappresenti l’ottimalita in senso quadraticoe dato dal modo in cui posso esprimere un elemento di uno spazio di Hilbert.

7

Sia X uno spazio di Hilbert separabile dotato di prodotto interno (·, ·) edi una base ortonormale φi. Ciasun elemento di f (x, t) ∈ X puo essererappresentato come

f (x, t) =∑i

ai (t)φi (x) =∑i

(f (x, t) , φi (x))φi (x) . (11)

I coefficienti ai sono comunemente chiamati coefficienti di Fourier . L’ideae quella di cercare una base ortonormale φi tale che tale che le medie deicoefficienti di Fourier siano ordinati come

〈a21 (t)〉 ≥ 〈a2

2 (t)〉 ≥ . . . , (12)

dove 〈·〉 un operatore di media. Ritornando allo spazio M, possiamo quindiinterpretare, grazie all’ortonormalita di V, l’errore

∥∥S−VVTS∥∥F

come lasomma dei quadrati dei coefficienti dei termini non considerati e con la tecnicaPOD ricerchiamo la matrice V tale per cui questo errore e il minore tratutte le matrici appartenenti a Rn×k. Per la maggior parte delle applicazionipratiche, l’utilizzo di pochi elementi permette di rappresentare anche piu del99% del contenuto in termini di quadrati dei coefficienti.

Nonostante il problema (10) sembri un generico problema di minimiz-zazione lineare, il teorema di Schimdt-Mirsky-Eckart-Young ne caratterizzala soluzione come la decomposizione a valori singolari (SVD) troncata dellamatrice Sw. L’approssimazione a valori singolari, troncata a k termini, dellamatrice S e data da

σ1u1vT1 + · · ·+ σkukv

Tk , (13)

dove σi, ui e vi sono rispettivamente i valori singolari, i vettori singolarisinistri e i vettori singolari destri di S.

In [ Algoritmo 1 ] viene sintetizzato l’algoritmo POD cosı presentato,

Algorithm 1 POD

Require: Snapshots y (t0) , . . . ,y (tn)1: Definisci Y = [y (t0) , . . . ,y (tn)]2: Calcola la SVD di Y, Y = UΣV∗

3: Considera le prime k colonne della matrice U = U (:, 1 : k).1

4: Definisci la matrice di riduzione POD come ΨPOD = U

1E stata usata la notazione Matlab per la selezione di righe e colonne di una matrice.

8

basato sulla decomposizione a valori singolari della matrice degli snapshots.Nella sua prima versione, esso era basato sulla definizione di una matrice dicorrelazione tra i data raccolti. Tuttavia, cio che accomuna tutti i metodi diproiezione POD e la ricerca di un operatore di proiezione che minimizzi ladistanza tra la soluzione completa e la sua proiezione sullo spazio ridotto.

3.1. Trattamento termine nonlineare

Il sistema ridotto (9) ci permette di illustrare l’inefficienza che si presentaquando dobbiamo risolvere un sistema nonlineare ottenuto tramite il metodoPOD-Galerkin. Il termine lineare viene efficacemente semplificato nella forma

(V)†A (t, ω) V = A (t, ω) , (14)

dove A (t, ω) ∈ Rk×k. Le matrici A (t, ω) possono essere precomputate du-rante la fase offline e, durante l’integrazione numerica del modello ridotto, ilcosto computazionale dipende solo da k, numero delle basi scelte, e non dan, dimensione del sistema originale. Il contributo nonlineare invece, dato da

(V)† fnl (t,Vy, ω) , (15)

dipende anche da n. In particolar modo, se valutare la nonlinearita richiedeun numero di operazioni pari a O (α (n)), con α (·) funzione di n, allora va-lutare la nonlinearita nel modello ridotto costa O (α (n) + 4nk). In caso dinonlinearita, quindi, il modello ridotto ha almeno lo stessa complessita delmodello di partenza e viene fallito l’obiettivo principale della riduzione dimodello. A livello implementativo quello si presenta un problema di lift-ing bottleneck, ovvero, nonostante la soluzione ridotta viva in uno spazio k-dimensionale, e necessario tornare allo spazio completo (lifting) per valutarela nonlinearita. Un problema simile si presenta anche quando la riduzione dimodello non riguarda sistemi dinamici, ma semplicemente equazioni dipen-denti da un parametro. In questo caso e richiesta la valutazione della matriceJacobiana che dipende da n.Diversi metodi sono stati introdotti per la risoluzione di questo problema[12][13][14][15][16]. In questo report ci concentriamo sul DEIM [15], forse ilpiu famoso dei metodi diquesta famiglia, e che rappresenta una variante disc-reta dell’EIM [13]. Quest ultimo e stato introdotto per la approssimazionedi funzioni parametriche non affini [17][18], presenti in operatori completi,prima di effettuare la riduzione di modello. Allo stesso modo anche il DEIMha una formulazione piu generale del setting di riduzione di modello di cui

9

abbiamo parlato fino ad ora. Data f (ω), funzione vettoriale del parametro,il metodo DEIM viene usato per la sua approssimazione e valutazione.Nel caso di nonlinearita polinomiali di alto ordine, o non-polinomiali, lariduzione di sistema prende il nome di iper-riduzione. L’idea del metodoDEIM e molto simile a quanto visto con il metodo POD per la riduzione dellasoluzione del sistema (1), ma applicata alla funzione non lineare f (Vy, t, ω).Essa viene proiettata su di un sottospazio di dimensione m dello spazio di di-mensione n a cui appartiene il termine nonlineare, con m n. Introduciamoquindi

Mfnl = g (t,x, ω)| t ∈ R,x ∈ Rn, ω ∈ Γ (16)

come lo spazio in cui vive il termine nonlineare. Esso svolgera esattamente lostesso ruolo cheM ha svolto per la tecnica POD, ovvero sara l’obiettivo dellaapprossimazione con uno spazio di dimensione fortemente ridotta. Suppo-nendo di disporre delle basi dello spazio ridotto ed indicandole con u1, . . . , um,un generico elemento f (t,x, ω) diMfnl viene rappresentato su questo spaziotramite

f (t,x, ω) ≈ Uc (t,x, ω) , (17)

dove U presenta come colonne le basi dello spazio ridotto e c sono i coeffi-cienti di f rispetto a queste basi. Supponiamo quindi di disporre di un setp1, . . . , pmdi m indici da 1 a n e di rappresentare con

P = [ep1 , . . . , epm ] (18)

la matrice avente come i-esima colonna la pi-esima colonna della matriceidentita. L’effetto della moltiplicazione tra P ed il termine f e quello di se-lezionare le componenti corrispondenti agli indici precedentemente introdotti.Moltiplicando (17) per P T otteniamo un sistema non piu sovradeterminatodella forma

P T f (t,x, ω) ≈ P TUc (t,x, ω) . (19)

Si puo dimostrare [19] che la matrice P TU risulta essere invertibile e quindi(19) puo essere usata per trovare il valore di c. L’approssimazione di f equindi data da

f(t,x, ω) ≈ Uc = U(P TU)−1P T f(t,x, ω), (20)

tecnica comunemente nota come Discrete Empirical Interpolation. Il sistemacosı ridotto ha quindi la forma

d

dtx = (V)†A(ω)Vy + V†U(P TU)−1P T f(t,x, ω). (21)

10

In questo caso, durante la fase offline del problema, la matrice V†U(P TU)−1

puo essere pre-computata e, dato che P T f e un vettore di m componenti,(21) non dipende dalla dimensione del sistema originale.La matrice U viene originata seguendo lo stesso procedimento utilizzato perV, applicato alla matrice degli snapshots delle nonlinearita, che puo esseresalvata durante la computazione o ottenuta a posteriori a partire da S.La computazione degli indici p1, . . . , pm richiede invece un procedimento iter-ativo di selezione delle componenti che limiti l’errore di interpolazione. Datala matrice di proiezione U, il primo indice di interpolazione p1 e dato dallacomponente di u1 piu grande in valore assoluto. I restanti indici consideratisono dati dalle componenti di maggiore ampiezza nel residuo r = ul −Uc.La prova dell’invertibilita di P TU deriva dal fatto che il DEIM garantisce unresiduo non nullo ad ogni passo. Lo pseudocodice per questo algoritmo vienefornito in [ Algoritmo 2 ].

Algorithm 2 DEIM

Require: Sistema di basi del termine nonlineare u1, . . . , um.1: Prendi p1 come l’indice della componente piu grande di u1.2: U← [u1].3: P ← [p1].4: for i← 2 to m do5: Risolvi P TUc = P Tui per trovare c.6: ri ← ui −Uc.7: Prendi pi come l’indice della componente piu grande di ri.8: U← [u1 . . . ui]9: P ← [p1 . . . pi]

10: end for

4. Proper Symplectic Decomposition (PSD)

Nei primi decenni dopo l’applicazione di metodi a basi ridotte per PDE, ilfocus principale della ricerca nel settore era quello di migliorare l’accuratezzanell’approssimazione del sistema di partenza. Il metodo POD, con il relativotrattamenti per i termini nonlineari, rappresenta un chiaro esempio di questiapproccio: la matrice di base ridotta viene definita come la soluzione di unproblema di ottimizzazione riguardante una funzione errore della proiezionedella matrice degli snapshots. Diverse stime di convergenza si possono trovare,

11

per il caso lineare e non, in [18]. Molto piu di recente l’attenzione si e spostatasulla conservazione di strutture geometriche della equazione di partenza,garantendo solidi risultati in termini di stabilita del sistema. Il prezzo piuoneroso da pagare, rispetto ai metodi sopra citati, l’aggiunta di ulteriorivincoli alle funzioni errore precedentemente descritte. I primi risultati inquesto settore sono stati raggiunti da Lall, Krysl e Marsden in Lall et al.[20], dove la struttura geometrica considerata quella Lagrangiana. Comespesso accade per questo tipo di sistemi, la descrizione della dinamica delproblema e espressa in termini di una coppia data da un elemento di unospazio di configurazioni e da un elemento del tangente spazio delle fasi. Ilcaso piu semplice di struttura Lagrangiana, usata in meccanica, e costituitadalla posizione e velocita in R3. Nel caso di sistemi lineari, si dimostra checonsiderando una proiezione dello spazio delle configurazioni ed un liftingdi questa proiezione allo spazio delle fasi, si ottiene un sistema ridotto chepreserva l’originale struttura Lagrangiana. In questo lavoro ci concentriamosulle strutture Hamiltoniane, un tipo diverso di strutture, comunque legatea quella Lagrangiana. In termini di ODE, un sistema Hamiltoniano canonicopresenta la forma descritta nella sequente definizione.

Definizione 1. (Sistema Simplettico)Un sistema Hamiltoniano finito-dimensionale viene descritto dalla relazione

z = J2n∇zH (z) ,

z (0) = z0,(22)

dove H : R2n 7→ R e la funzione Hamiltoniana, z0 e la condizione iniziale e

J2n =

[0n In−In 0n

]. (23)

Introducendo le coordinate canoniche z =(qT ,pT

)T, dove q e p rappre-

sentano le posizioni ed i momenti generalizzati, il sistema (22) puo essereriscritto come

q = ∇pH (q,p) ,

p = −∇qH (q,p) .(24)

Prima di descrivere metodi di riduzione di equazioni (24), introduciamoalcune nozioni di base riguardanti i sistemi Hamiltoniani e la geometria sim-plettica.

12

Sia M un manifold regolare e Ω : M×M 7→ R, una 2-form chiusa, nondegenere e anti-simmetrica definita su M. Nella maggior parte delle appli-cazioni pratiche, il manifoldM e dato da uno spazio vettoriale di dimensione2n. La coppia (M,Ω) prende il nome di spazio vettoriale simplettico. De-finendo un sistema di coordinate canoniche suM come (q1, . . . , qn, p1, . . . , pn)la forma Ω puo essere scritta come

Ω =n∑i=1

dqi ∧ dpi, (25)

dove ∧ rappresenta il prodotto esterno. Per semplicita assumiamo che ilmanifold sia definito su R, e quindi Ω =M2n. Questo fa si che, dati v1,v2 ∈R2n , otteniamo una versione facilmente computabile della 2-forma

Ω (v1, v2) = vT1 J2nv2. (26)

Il risultato (26) e stato fortemente utilizzato in [21] per lo sviluppo di metodidi riduzione simplettica di tipo greedy.La seguente definizione introduce l’elemento fondamentale per la construzionedi operatori che preservano la struttura generica.

Definizione 2. (Mappe Simplettiche)Siano (M1,Ω1) e (M2,Ω2) due spazi vettoriali lienari e simplettici di dimen-sioni 2n1 e 2n2, rispettivamente. Una mappa lineare φ : M1 7→ M2 si dicesimplettica se

Ω1 = φ∗Ω2 (27)

dove φ∗Ω2 e il pullback di Ω2 dato da φ,i.e., per ogni z1, z2 ∈M1

Ω1 (z1, z2) = Ω2 (φ (z1) , φ (z2)) (28)

Nel caso di spazio finito dimensionale, che e quello in considerazione nelcaso di sistemi di ODE, abbiamo che (2) puo essere riespressa in forma ma-triciale come

ΦTJ2n1Φ = J2n2 (29)

dove Φ ∈ R2n1×2n2 rappresenta l’azione dell’operatore φ. Piu in generale,una matrice Φ che soddisfa (29) viene definita come matrice simplettica. Ilconcetto di matrice, o trasformazione simplettica, e piu generale del ramodella riduzione di modello: ad esempio, nella sua versione nonlineare, viene

13

sfruttato nell’ambito dell’integrazione numerica.Definiamo come Ψt il flusso dell’ODE (24). E possibile dimostrare che questoflusso e un morfismo simplettico, ovvero conserva la 2-form Ω, l’Hamiltonianodel sistema H e il volume del flusso. Fin dai pionieristici lavori di Feng Kangand Ruth, l’obiettivo e stato quello di costruire flussi numerici Ψh

t che ap-prossimassero Ψt e che, in particolar modo, conservassero la simpletticitadell’operatore. Restringendoci alla classe degli integratori di tipo Runge-Kutta, in particolar modo menzioniamo il metodo Stormer-Verlet, partic-ularmente utilizzato nell’ambito di applicazioni chimiche e astronomiche, egli schemi di collocazione di Gauss. Nei nostri esperimenti numerici con-sidereremo il metodo del punto medio implicito, il piu semplice esempio diintegratore di Gauss, per la sua versatilita: esso puo essere anche applicatoad un problema non simplettico, rimuovendo alcuni biasses che potrebberorendere non equo il confronto tra i vari metodi di riduzione. Sottolineamo chela maggior parte dei comuni metodi di tipo Runge-Kutta non sono simplet-tici e quindi non preservano la struttura geometrica del problema originale:in [22] viene proposta una condizione di simpletticita sui coefficienti dellatabella di Butcher che descrive il metodo.Tornando al semplice caso di mappe lineari, definiamo il concetto di inversadi una mappa simplettica.

Definizione 3. (Matrice Inversa Simplettica)L’inversa simplettica di una matrice Φ ∈ R2n1×2n2, indicata con Φ+, edefinita come

Φ+ = JT2n2ΦTJ2n1 . (30)

Inoltre (Φ+)T

e anche essa una matrice simplettica e vale Φ+Φ = I2n2.

E facile verificare che ΦΦ+ e una matrice idempotente, un requisito fon-damentale per il ruolo di proiettore su di uno spazio di dimensione ridotta.

4.1. Riduzione di modello simplettica

Ora che le basi di riguardanti gli spazi geometrici necessari sono statedefinite, ci concentriamo sulla riduzione di modello. Come nel caso di PODe DMD, l’assunzione fondamentale e che esista uno spazio ridotto Mrb, didimensione contenuta 2k, che approssimi lo spazio originale, di dimensione2n. In questo contesto, aggiungiamo la richiesta che sia possibile definireuna struttura simplettica su di esso, per avere (Mrb,Ωrb). Definita una base

14

simplettica Φrb, possiamo approssimare la soluzione del sistema (22) tramite

z ≈ Φrby. (31)

Da questo, riscriviamo il termine destro del problema come

∇zH (Φrby) = ΦTrb∇yH (Φrby) . (32)

Considerando quindi la derivata temporale di (31) e la caratterizzazione dellematrici simplettiche abbiamo che, partendo dalla relazione inversa di (31),

y = Φ+rbz

= Φ+rbJ2n∇zH (z)

= J2kΦT∇xH (x)

= J2k∇yH (Φy) . (33)

L’equazione (33) rappresenta il problema ridotto di (22) e mantiene la strut-tura simplettica se e solo se la matrice di basi ridotte Ψ e simplettica.Il problema adesso si sposta sulla computazione di una matrice simpletticaΨ che sia anche una base dello spazio da ridurre. Diverse tecniche sono stateproposte negli ultimi anni, sia basate sulla decomposizione a valori singolarisia su algoritmi greedy. In questo report consideriamo il metodo conosciutocome cotangent lift, il cui nome deriva dalla riformulazione che si puo faredi (quasi) tutti i sistemi simplettici come spazi cotangenti di uno spazio diconfigurazione.Partendo da un caso generico, consideriamo due spazi vettoriali R e Q didimensioni n e k, rispettivamente, e k ≤ n. Introduciamo W = T ∗R eV = T ∗Q, gli spazi cotangenti relativi a R e Q. Dato uno spazio vettoriale,il relativo spazio cotangente e dato da dalle coppie formate da punti dellospazio vettoriale e i correspondenti duali del piano tangente. Allo stessomodo, gli spazi tangenti relativi a R e Q, denotati con TR e TQ, sono for-mati dalle coppie di punti dello spazio e dai relativi piani tangenti. In questasezione l’idea e quella di fornire una descrizione generica dietro all’algoritmocotangent lift, senza, dove e possibile, scendere in dettagli riguardanti risul-tati avanzati di geometria differenziale. Consideriamo due funzioni linearif : R 7→ Q e π : Q 7→ R tali che π f = idR, ovvero la loro composizionecorrisponde all’operatore identita definito su R. Nel contesto di riduzione dimodello, f rappresenta una riduzione dell’originale spazio di configurazione.

15

Da questo, considerando che un generico spazio simplettico finito dimension-ale e rappresentabile come spazio cotangente di uno spazio vettoriale (cherappresenta nel nostro caso lo spazio delle configurazioni), il passo successivoe quello di sollevare la relazione tra R e Q ai relativi spazi cotangenti, dacui il nome cotangent lift. Il passo intermedio e quello di sollevare la re-lazione allo spazio tangente: si puo dimostrare che, data la linearita di f ,il suo tangent lift f∗ coincide con f stesso. Intuitivamente, questo si puospiegare considerando che lo spazio tangente e formato dai piani tangentiai punti dello spazio delle configurazioni: essendo esso generato da derivatedirezionali e sapendo che queste non influenzano operatori lineari, abbiamoil risultato sopra descritto.Introduciamo quindi il pairing naturale tra spazi cotangenti e tangenti rel-ativamente a punti degli spazi di configurazioni, indicandoli con 〈·, ·〉q :=T ∗qQ×TqQ 7→ R per q = f (r) ∈ Q e 〈·, ·〉r := T ∗r R×TrR 7→ R per r ∈ R. Inatural pairings derivano dalla definizione di spazio cotangente: considerandoun termine dello spazio cotangente ed uno dello spazio tangente, essi restitu-iscono come risultato l’applicazione del primo sul secondo. Il cotangent lift el’unico operatore f ∗ : V 7→W che soddisfa

〈f∗|q (p) ,w〉r = 〈p, f∗r (w)〉q. (34)

Introducendo un sistema di coordinate su R e Q, la relazione (34) puo essererappresentata in forma matriciale. Usando coordinate canoniche, consideri-amo r, w ∈ Rk, q, pq ∈ Rn come vettori di coefficienti di elementi di R e Q.Le mappe lineari f e π sono date dalle matrici Φ ∈ Rn×k e ΨT ∈ Rn×k. Lacondizione di identita sulla composizione degli operatori viene quindi espressacome ΨTΦ = Ik. Ricordando che il tangent lift di una funzione lineare rimanelineare abbiamo che f∗r (w) = Φw. Da questo, abbiamo che, identificando scon la componente appartenente allo spazio cotangente del cotangent lift dif valutato in p,

s = ΦTp. (35)

Consideriamo quindi il vettore x = [q; p] ∈ R2n, rappresentante i coefficientidi un elemento dello spazio simplettico da ridurre (chiamato W): se la mappaf rappresenta una riduzione dello spazio originario in uno spazio simpletticoridotto (indicato con V) tale che z = [r;w] ∈ R2k e la proiezione di x, allora

16

abbiamo che [rs

]= z = f (x) = f

([qp

])=

[ΨT 00 ΦT

]︸︷︷︸

B+

[qp

]. (36)

Indicando con B l’inversa simplettica di B, essa puo essere scritta come

B = JT2nB+J2k =

[Φ 00 Ψ

]. (37)

Data la proprieta sul prodotto tra Φ e ΨT e facilmente dimostrabile cheB+B = I2k. Questa e una condizione sufficiente per far si che B sia simplet-tica.Consideriamo lo scenario in cui Ψ = Φ e quindi abbiamo che il vincoloΨTΦ = ΦTΦ = Ik rappresenta una condizione di ortogonalita sulla matriceΦ. Da un punto di vista della riduzione questo fa si che, essendo Φ unamatrice di proiezione, il suo spazio colonna debba rappresentare una basesia per p che per q. Un possibile algoritmo per la costruzione di una buonamatrice in termini di accuratezza e basato sulla decomposizione a valorisingolari. Definiamo come SCL la matrice degli snapshots data da

SCL = [pij, qij], (38)

dove gli indici i e j servono per listare gli snapshots in funzione del numerodi frame per tempo e del numero di parametri considerati per il problemamulti-parametrico. Consideriamo il problema di ottimizzazione

minΦ∈Rn×k,ΦTΦ=Ik

‖SCL −ΦΦTSCL‖. (39)

Grazie ancora al teorema di Schmidt-Mirsky-Eckart-Young, possiamo usareuna SVD di SCL per ottenere la base Φ. Il metodo cosı ottenuto vienechiamato Cotangent Lift ed e riassunto nel seguente pseudocodice [ Algoritmo3 ].

La stabilita di un sistema simplettico e intesa nel senso della seguentedefinizione

Definizione 4. Consideriamo un sistema dinamico della forma (1) e sup-poniamo che xe sia un punto di equilibrio del sistema tale che f(xe)=0. xe eLyapunov stabile se, per ogni ε > 0, possiamo trovare un δ > 0 tale per cuiper ogni traiettoria φt, se ‖φt − xe‖ ≤ δ, allora abbiamo che ‖φt − xe‖ ≤ εper 0 ≤ t ≤ ∞.

17

Algorithm 3 Cotangent Lift

Require: Snapshots qi (t0) , . . . ,qi (tn) ,pi (t0) , . . . ,pi (tn)1: Definisci la matrice degli snapshots SCL2: Calcola la SVD di SCL, SCL = ΦΣV∗

3: Definisci la matrice di proiezione

B =

[Φ 00 Φ

]

In particolar modo, una forte indicazione circa la stabilita e data da

Definizione 5. Un punto di equilibrio e Lyapunov stabile se esiste una fun-zione scalare W : Rn 7→ R tale che ∇W (xe) = 0, ∇2W (xe) e definita positivae, in un intorno di xe, abbiamo che d

dtW ≤ 0.

Si dimostra che nel caso di sistemi simplettici questo ruolo viene giocatodalla funzione Hamiltoniana ( opportunamente aggiustata di segno ). Eccoquindi che, mantenendo la struttura simplettica del problema anche nel casoridotto, questa proprieta di stabilita vale per la soluzione ridotta.

4.2. Trattamento termine nonlineare

Consideriamo quindi un generico sistema Hamiltoniano e la sua formaridotta, considerando la matrice di proiezione simplettica B. Dividiamo lafunzione Hamiltoniana in due contributi H = H1 +H2 tali che ∇H1y = Lye ∇H2 = g(y), dove L ∈ R2n×2n e una matrice costante e g e una funzionenonlineare. Il sistema ridotto ha la forma

d

dty = A†J2nLA︸︷︷︸

L

y + A†J2ng(Ay) (40)

e, applicando il metodo DEIM nel caso di nonlinearita, otteniamo

d

dty = A†J2nLA︸︷︷︸

L

y + A†J2nU(P TU)−1P Tg(Ay)︸︷︷︸N(y)

. (41)

Una generica scelta di U comporta la perdita della simpletticita del sistema,con relative conseguenze in termini di stabilita sul lungo termine. E possibile

18

dimostrare che il sistema ottenuto e Hamiltoniano se andiamo a considerareU = (A†)T . Come primo passo prendiamo

g(Ay) = ∇zH2(z) = (A+)†∇yH2(Ay) (42)

dove abbiamo usando la chain rule nella derivazione. Sostituendo (42) nellaformulazione di N(y) abbiamo che

N(y) = A†J2nU(P TU)−1P T (A+)†∇yH2(Ay). (43)

Prendendo U = (A+)T abbiamo che

N(y) = A+J2n(A+)†∇yH2(Ay) = J2k∇yH2(Ay), (44)

dove nell’ultimo passaggio abbiamo usato il fatto che (A+)T sia una matricesimplettica. L’equazione (44) afferma che il termine nonlineare ha la formasimplettica richiesta e quindi U = (A†)T e una condizione sufficiente per ilmetodo DEIM di mantenere la simpletticita del sistema.

5. Implementazione

5.1. Compilare, installare ed usare la libreria

La libreria e costituita quasi interamente da header files e inline func-tions e viene sfruttutata la libreria Eigen 3.3.4 per le sue funzionalita relativeall’algebra lineare e l’interfaccia gnuplot-iostream per visualizzare i dati du-rante la simulazione. Sebbene queste librerie non richiedano un’installazionedi alcun tipo di per se, gnuplot-iostream richiede la libreria Boost, moltocomune e disponibile nativamente sui compilatori di diverse distribuzioniUnix. Sebbene la maggior parte delle librerie Boost risulti essere header-only, in questo caso abbiamo bisogno di Boost.IOStreams principalmente, edi altre librerie secondarie, che invece devono essere costruite separatamente.Per ottenerle nella versione piu aggiornata, come suggerito dal relativo sito,bisogna:

• Scaricare da SourceForge la versione piu recente della libreria (fileboost 1 61 0.tar.bz2 al momento).

• Nella directory in cui vogliamo installare Boost, eseguire

ta r −−bzip2 −xf /path/ to / boo s t 1 61 0 . ta r . bz2

19

http://www.sharelatex.com

Alternativamente, la versione disponibile nella universe repository, usata pertestare il codice, e ottenibile tramite

sudo apt−get i n s t a l l l i bboo s t−a l l−dev

Per utilizzare la libreria in un proprio codice, bisogna linkare le correttelibrerie Boost, includendo come flags

−l b o o s t i o s t r e a m s −l boo s t sy s t em − l b o o s t f i l e s y s t e m

e assicurarsi che esse si trovino nel path corretto. Per quanto riguarda Eigen,gnuplot-iostream e la libreria per la riduzione di modello, esse sono fornitenella sottocartella include del progetto. Per semplificare l’utilizzo della li-breria, alcuni file, costituiti da include di header files, sono stati realizzati everrano considerati nel seguito.Per compilare i problemi test realizzati e possibile utilizzare il Makefile for-nito nella cartella del progetto, per il quale, tramite il comando make verragenerato un object file per ciascuno dei test contenuti nelle sottocartelle disrc e posizionato nella cartella build. La fase successiva di linkage, genererail relativo eseguibile nella cartella bin. E possibile compilare e linkare i testsingolarmente, andando ad utilizzare i Makefiles contenuti nelle rispettivecartelle dei test: operativamente, quello che viene fatto dal Makefile gen-erale e chiamare ricorsivamente i Makefiles degli esempi. Il comando makeclean rimuove gli object files dalla cartella build. Il compilatore testato egcc5.4.0 con optimization flags -O3 e -march=native: sconsigliamo di uti-lizzare -ffast-math perche abbiamo notato una degradazione importantedell’accuratezza nella decomposizione a valori singolari che utilizzeremo nelseguito.Prima di descrivere in dettaglio le parti fondamentali della libreria, descrivi-amo brevemente alcune funzioni template, implementate per ampliare quelloche viene gia offerto da Eigen.

5.2. Estensione delle funzionalita di Eigen

C++ non permette di aggiungere metodi ad una classe gia esistente, quindisi e deciso di creare un header file contenente alcune funzioni template einline, non fornite nativamente da Eigen, e utilizzate di frequente nel codice.Di seguito forniamo una descrizione delle piu importanti. Per facilita dicomprensione del codice, abbiamo introdotto i seguenti typedefs

1 typedef Matrix<double,Dynamic,Dynamic> Matrix_Real;

2 typedef Matrix<complex<double>,Dynamic,Dynamic> Matrix_Complex;

20

3 typedef Matrix<double,Dynamic,1> ColVector_Real;

4 typedef Matrix<double,1,Dynamic> RowVector_Real;

5 typedef Matrix<long int,Dynamic,1> ColVector_Int;

per indicare diversi tipi di data storage.Una funzionalita che Eigen non fornisce nativamente, e quella di poter leggeree scrivere su file il contenuto di una matrice. Questa funzione e invece impor-tante nella riduzione di modello, dove il salvare i risultati delle simulazioni delmodello completo per poi assemblarle in seguito come matrice degli snapshotsrisulta essere importante. Per questo motivo abbiamo definito le funzioniwrite binary e read binary. Esse accettano come parametri una stringarappresentante l’indirizzo del file in cui andare a scrivere/leggere, realtivo allaposizione dell’eseguibile, e la referenza all’oggetto da salvare/sovrascrivere.

1 template <class T>

2 void write_binary( string Filename, const T& matrixToWrite )

3

4 ofstream out( Filename, ios::out | ios::binary | ios::trunc );

5 typename T::Index rows=matrixToWrite.rows(),

6 cols=matrixToWrite.cols();

7 out.write((char*) (&rows), sizeof(typename T::Index));

8 out.write((char*) (&cols), sizeof(typename T::Index));

9 out.write((char*) matrixToWrite.data(),

10 rows*cols*sizeof(typename T::Scalar));

11 out.close();

12 ;

13

14 template <class T>

15 void read_binary( string Filename, T& matrixToRead )

16

17 ifstream in( Filename, ios::in | ios::binary);

18 typename T::Index rows=0, cols=0;

19 in.read((char*) (&rows), sizeof(typename T::Index));

20 in.read((char*) (&cols), sizeof(typename T::Index));

21 matrixToRead.resize( rows, cols );

22 in.read((char *) matrixToRead.data(),

23 rows*cols*sizeof(typename T::Scalar));

24 in.close();

25 ;

21

Un’altra importante feature richiesta dalla riduzione di modello e quelladi poter accedere ad una matrice di Eigen tramite un sottoinsieme di indici.Nell’algoritmo relativo al DEIM, frequentemente viene richiesto, P TU o P T f ,ovvero prodotti tra la trasposta della matrice degli indici e la base ridotta oil termine nonaffine. Mentre nel primo caso si avrebbe un semplice prodottosparse-dense matriciale, nel secondo si dovrebbe conoscere l’intero vettoref , rendendo quindi inutile la procedura di interpolazione quando in realtabasterebbe un accesso tramite indici. Eigen permette di fare cio tramiteuna speciale classe chiamata CwiseNullarOp, utilizzata prevalentemente perdefinire matrici procedurali, come la matrice costante, la matrice identita ouna matrice randomica, ma che puo essere anche utilizzata per accedere aduna matrice tramite indici. Il metodo per accedere tramite indici ad unamatrice di Eigen viene descritto in un’apposita sezione della guida online elo riportiamo di seguito. L’idea e quella di definire un nullary functor persalvare come stato una referenza alla matrice a cui bisogna accedere, gli indiciriga e gli indici colonna. Viene poi creata la nullary expression che utilizzeral’operatore () del functor.

1 template<class ArgType, class RowIndexType, class ColIndexType>

2 class indexing_functor

3 const ArgType &m_arg;

4 const RowIndexType &m_rowIndices;

5 const ColIndexType &m_colIndices;

6 public:

7 typedef Matrix<typename ArgType::Scalar,

8 RowIndexType::SizeAtCompileTime,

9 ColIndexType::SizeAtCompileTime,

10 ArgType::Flags&RowMajorBit?RowMajor:ColMajor,

11 RowIndexType::MaxSizeAtCompileTime,

12 ColIndexType::MaxSizeAtCompileTime> MatrixType;

13 indexing_functor(const ArgType& arg,

14 const RowIndexType& row_indices,

15 const ColIndexType& col_indices)

16 : m_arg(arg), m_rowIndices(row_indices),

17 m_colIndices(col_indices)

18

19 const typename ArgType::Scalar& operator() (Index row,

20 Index col) const

22

21 return m_arg(m_rowIndices[row], m_colIndices[col]);

22

23 ;

24

25 template <class ArgType, class RowIndexType, class ColIndexType>

26 CwiseNullaryOp<indexing_functor<ArgType,RowIndexType,ColIndexType>,

27 typename indexing_functor<ArgType,RowIndexType,

28 ColIndexType>::MatrixType>

29 indexing(const Eigen::MatrixBase<ArgType>& arg,

30 const RowIndexType& row_indices,

31 const ColIndexType& col_indices)

32

33 typedef indexing_functor<ArgType,RowIndexType,ColIndexType> Func;

34 typedef typename Func::MatrixType MatrixType;

35 return MatrixType::NullaryExpr(row_indices.size(),

36 col_indices.size(), Func(arg.derived(),

37 row_indices, col_indices));

38

Altre funzioni di interesse in questa parte della libreria riguardano lagenerazione di una griglia di valori dati due vettori in ingresso, imitando lefunzionalita della funzione meshgrid di Matlab. Utilizzare il tipo di datoin input come template, permette inoltre di ricevere in ingresso oggetti piugenerali, come matrici, e svilupparli come vettori prima di generare la griglia,anche se questa funzionalita non viene sfruttata nel nostro codice. La parteimportante della nostra implementazione di meshgrid e l’utilizzo della fun-zione reshape, basata sull’oggetto mappa fornito da Eigen. Esso non e altroche una expression per definire la mappatura di una matrice in un’altra aventelo stesso numero di elementi ma in una disposizione differente.

1 template<typename T1,typename T2>

2 void reshape( T1& reshapeMatrix, T2& dataMatrix,

3 const long int numberOfRows,

4 const long int numberOfCols )

5

6 reshapeMatrix.resize(numberOfRows,numberOfCols);

7 reshapeMatrix = Map<T1>(dataMatrix.data(),

8 numberOfRows,numberOfCols);

9 ;

23

Nel caso di dimensioni non compatibili, un run time assert viene lanciatoda Eigen per segnalare il problema.

5.3. Struttura della libreria

Per descrivere la struttura scelta per la libreria e necessario definire primail workflow del processo di riduzione.

Sistemadinam-

ico

SolverODE

Analisi ecollezionedei dati

Riduzionedei dati

Ridurre una PDE tramite un metodo a basi ridotte significa andare a ridurreil relativo sistema di ODEs ottenuto dopo una discretizzazione spaziale.Bisogna quindi definire un oggetto che rappresenti il sottostante sistema di-namico, andando principalmente a definire la derivata temporale della vari-abile incognita.Questo oggetto si deve interfacciare con l’integratore numerico, il quale rap-presenta la parte piu costosa in termini computazionali dell’intero processoe quindi deve essere il piu efficiente possibile.Allo stesso tempo viene data all’utente la possibilita di analizzare i dati inreal-time, andando da una semplice rappresentazione di un sottoinsieme dellasoluzione al monitoraggio di alcune grandezze di interesse, come ad esempiola funzione Hamiltoniana nel caso dei sistemi simplettici.Una volta raccolti i dati, tramite il processo di riduzione si genera il sistemadinamico ridotto, il quale e strutturalmente uguale al sistema di partenzae quindi deve poter essere risolto ed analizzato seguendo lo stesso procedi-mento. La spiegazione della libreria segue l’andamento di questo workflow,partendo quindi dal modo in cui viene definito un sistema dinamico.

24

5.4. Sistema dinamico

Il piu generico sistema dinamico che abbiamo a disposizione prende laforma

dx

dt= f(x(t), t), (45)

dove l’evoluzione di una generica quantita x viene espressa in funzione dellasua derivata temporale f , la quale verra chiamata diverse volte per la risoluzionedi (45) tramite un integratore numerico. Possibile scelte implementative perla rappresentazione di f sono lambda functions, puntatori a funzioni e func-tors, ma nel codice abbiamo preferito usare questi ultimi. I motivi principalidi questa scelta sono:

• Stato. Essendo fondamentalmente classi in cui viene definito l’operatore(), permettono di aggiungere membri della classe utili per la com-putazione di f stessa, come ad esempio matrici di rigidezza e altrioperatori per il calcolo di contributi nonlineari, passati all’oggetto so-lamente durante la sua costruzione.

• Polimorfismo. Dovendo il solver ODE essere in grado di risolvereun generico sistema dinamico, risulta utile introdurre un sistema di-namico come classe puramente astratta, avente il ruolo di interfaccia,e poi sfruttare il polimorfismo per passare tramite referenza/puntatoreoggetti di classi derivate in cui vengono definite funzionalita aggiuntive( ed in particolare in cui viene definito l’operatore () ).

• Performance. I functors vengono frequentemente messi inline dalcompilatore, mentre solitamente lo stesso non avviene per puntatoria funzioni ( vedi comparazione tra qsort e std::sort ), riducendoeventuali overhead nel primo caso. Non avendo a disposizione unaimplementazione per ognuna delle tecniche, questo criterio pero rimanepuramente teorico e non e stato verificato.

La semplice interfaccia e quindi data da

1 template<class Container>

2 class DynamicalSystem

3

4 typedef Eigen::Block<Container> matrixBlock;

5 public:

25

6 virtual ~DynamicalSystem() ;

7 virtual void operator()

8 (const double time,

9 Container& inputState,

10 matrixBlock&& outputState) = 0;

11 ;

dove l’operatore () deve accettare in ingresso lo stato del sistema (inputState),il tempo a cui si vuole calcolare f (time) e l’output in cui si vuole valutaref (outputState). Il passaggio per referenza dell’input viene fatto perche, incasi generali, lo stato potrebbe avere grandi dimensioni e quindi una copiasarebbe stata inefficiente. La spiegazione circa la natura di matrixBlock esul perche venga usata la move semantic viene rimandata alla sezione sulsolver ODE. Un modo per definire l’operatore () nella classe figlia e quellodi andare a scrivere una classe derivata di DynamicalSystem in cui vienefatto l’override di operator() come nel caso dei problemi di Lotka-Volterrae Kepler analizzati come test cases ( si vedano i file lotka volterra.hpp ekepler.hpp ).Dal momento che la riduzione di modello riguarda prevalentemente PDE dis-cretizzate in sistemi di ODE, abbiamo deciso di implementare classi appositeper rappresentare problemi lineari della forma

dx

dt= Ax(t) (46)

e nonlineari della forma

dx

dt= Ax(t) + g(x(t), t) (47)

che sono i tipici problemi ottenuti tramite metodo delle linee applicato ad unaPDE. Questa differenziazione e legata al modo in cui il processo di riduzioneviene applicato a sistemi lineari e non.Partendo dal caso lineare, riportiamo di seguito parte della sua definizione

1 template<class Container, class Affine>

2 class generalDynamicalSystem : public DynamicalSystem<Container>

3


5 public:

6 generalDynamicalSystem(Affine affineOperatorUser);

26

7 void operator()(const double time,


9 matrixBlock&& outputState) ;

10 private:

11 Affine affineOperator;

12 ;

L’unico membro aggiunto da questa classe e affineOperator, il qualerappresenta la matrice A del sistema dinamico. Il parametro template Affinee necessario perche la matrice A potrebbe essere sparsa, come succede aseguito di una discretizzazione usando elementi o differenze finite, oppuredensa, come quella che normalmente si ottiene a seguito di una riduzione dimodello ma l’unico motivo per cui viene utilizzata e quello di moltiplicarelo stato input per ottenere l’output ( operazione che non dipende dalla suanatura )

1 outputState.noalias() = affineOperator*inputState;

Il metodo noalias() viene utilizzato perche Eigen, in caso di moltipli-cazione tra matrici, suppone sempre che si possa essere nel caso in cui lastessa matrice possa apparire a destra e a sinistra dell’operatore di assegna-mento: viene quindi generato un temporary dell’oggetto sulla destra che poiviene assegnato a quello sulla sinistra, per evitare di sovrascrivere un valorementre si sta ancora calcolando ( problema di aliasing ). Nel nostro caso,non essendoci aliasing, il metodo noalias ci permette di evitare la creazionedell’oggetto temporaneo.Nel caso di sistema nonlineare ( o piu generalmente di sistema nonaffine )abbiamo a che fare con un’equazione del tipo (47) e quindi abbiamo bisognodi implementare anche il termine g. In questo caso, utilizziamo un ulterioreparametro template per ottenere un codice generico e parte della definizionedella relativa classe e data da

1 template<class Container, class Affine, class nonAffine>

2 class generalDynamicalSystemNonaffine :

3 public DynamicalSystem<Container>

4


6 public:

7 generalDynamicalSystemNonaffine(Affine affineOperatorUser,

27

8 nonAffine nonaffineOperatorUser);

9 void operator()(const double time,


11 matrixBlock&& outputState) ;

12 private:

13 Affine affineOperator;

14 nonAffine nonaffineOperator;

15 ;

dove l’operator () viene overridato dalla classe base per ottenere

1 outputState.noalias() = affineOperator * inputState

2 + nonaffineOperator( inputState, time );

Il nonaffineOperator e stato implementato come un semplice functor,la cui interfaccia e rappresentata dalla classe nonaffineOperator data ininclude/Dynamical System/nonaffine contribution.hpp e, sfruttando ilpolimorfismo, l’utente puo definire la propria classe derivata per implementareil contributo nonaffine ( si veda come esempio la definizione dell’operatore diSine-Gordon nella relativa cartella di esempio ).

5.4.1. Problemi trattati e differenze finite

Prima di procedere con la descrizione del solver ODE, ci soffermiamobrevemente sui problemi piu complessi trattati nei casi tests. Questi sonol’equazione lineare del calore in 1D, data da

∂u

∂t=∂2u

∂x2, (48)

l’equazione delle onde lineare in 2D∂p∂t

= q,∂q∂t

= c2∆p,(49)

e l’equazione nonlineare di Sine-Gordon in 2D∂p∂t

= q,∂q∂t

= c2∆p+ sin(p),(50)

dove negli ultimi due casi abbiamo utilizzato la formulazione simplettica. Intutti i casi e stato necessario discretizzare la derivata spaziale per ottenere

28

un sistema di ODE e questo e stato fatto, per semplicita, tramite un metodobasato su differenze finite, ambientato in un dominio rettangolare e period-ico. Nella nostra implementazione, ci siamo avvalsi della tecnica numerica digenerazione dello stencil descritta in [23] per ottenere uno schema di ordinearbitrario. Supponiamo che lo stencil per un determinato punto xi com-prenda tutti i punti xj tali che i − l ≤ j ≤ i + r, con l e r dati dal numerodi punti alla sinistra e destra, rispettivamente, e scelti dall’utente. Ad unnumero maggiore di punti di stencil ( in totale sono l + r + 1 ) corrispondeuna discretizzazione piu accurata. Data una funzione h(x), l’espansione diTaylor della sua derivata nel punto di indice i sulla griglia e data da

hi+m = hi +(m∆x)

1!hx +

(m∆x)2

2!hxx +

(m∆x)3

3!hxx + . . . , (51)

per m = −l, . . . , r. Moltiplicando ciascuna di queste espansioni per un coeffi-ciente am e sommando tutti questi contributi otteniamo la seguente equazione

r∑m=−l,m6=0

amhi+m −

(r∑

m=−l,m6=0

am

)hi =

(r∑

m=−l,m6=0

mam

)∆x

1!

∂h

∂x|i

+

(r∑

m=−l,m6=0

m2am

)∆x2

2!

∂2h

∂x2|i

+

(r∑

m=−l,m6=0

m3am

)∆x3

3!

∂3h

∂x3|i

+ · · · (52)

E chiaro che il coefficiente della k-esima derivata e dato da bk =∑r

m=−l,m6=0 mkam.

L’equazione (52) ci permette di determinare i l + r coefficienti am a secondadella derivata e dell’ordine desiderati. Se ad esempio vogliamo i coefficientidello stencil per la derivata del primo ordine al quarto ordine di accuratezzasettiamo b1 = 1 e b2,3,4 = 0. Genereremo quindi quattro equazioni e, di con-seguenza, almeno quattro punti saranno necessari per avere una soluzioneunivoca. Consideriamo il caso in cui r = 1 e l = 3 per ottenere uno schemadecentrato. Il sistema di equazioni da risolvere per i coefficienti prende quindila forma

−3 −2 −1 1(−3)2 (−2)2 (−1)2 (1)2

(−3)3 (−2)3 (−1)3 (1)3

(−3)4 (−2)4 (−1)4 (1)4

a−3

a−2

a−1

a1

=

1000

(53)

29

L’estensione al caso di ordine arbitrario non richiede spiegazioni. Questoalgoritmo viene implementato nella funzione computeStencil contenuta nelfile include/FD/finite difference.hpp.Una volta che abbiamo una funzione per calcolare lo stencil, possiamo fi-nalmente calcolare la matrice che rappresenta l’operatore di discretizzazionedella derivata in 1D. Questa viene fornita come output dalla funzione periodicFD1D,contenuta nello stesso file di computeStencil. Questa funzione semplice-mente cicla le righe della matrice di discretizzazione, tenendo conto dellecondizioni periodiche per i nodi vicini al bordo. Le uniche particolarita diquesta funzione sono che la matrice viene memorizzata secondo un ordina-mento per righe, contrariamente a quanto avviene di default per le matricidi Eigen

1 SparseMatrix<double,RowMajor> matrixFD1D(numberOfNodes,

2 numberOfNodes);

Il vantaggio di questo approccio e che viene reso piu efficiente qualsi-asi prodotto matriciale della forma matrixFD1D ∗ B, con B ordinata percolonne, e permette di allocare memoria preventivamente per la matricetramite reserve, dal momento che si conosce il numero esatto di elementiper riga ( pari alla lunghezza dello stencil ).Per ottenere la discretizzazione del secondo ordine, e sufficiente moltiplicarela matrice di discretizzazione della derivata prima per meno la trasposta dellamatrice stessa, come viene fatto nella funzione periodicFD1D Laplacian.La matrice di discretizzazione della derivata nel caso 2D viene ottenuta evi-tando numerosi loops che renderebbero il codice poco intuibile all’utente:sfruttiamo invece la possibilita di ottenere lo stesso risultato usando il prodottodi Kronecker tra la matrice di discretizzazione nel caso 1D e la matrice iden-tita. Il cambiare l’ordine dei termini in questo prodotto permette di costru-ire la matrice relativa alla direzione x e y: questo e quello che sostanzial-mente differenzia periodicFD2D Y e periodicFD2D X. La costruzione delLaplaciano in 2D sfrutta la stessa idea del caso 1D e viene eseguita inperiodicFD2D Laplacian.

5.5. ODE solver

La funzione f in (45) descrive la ODE, e sia x che f , per come e imple-mentata questa parte della libreria, possono essere vettori o matrici, grazieall’implementazione template. Una ODE risulta completa quando viene for-nita anche una condizione iniziale, andando quindi a definire un problema ai

30

valori iniziali (IVP). Il solver implementato e basato su due blocchi fonda-mentali:

• Una funzione stepper che calcola lo stato x(t+ dt) a partire da x(t)

• Una funzione integratrice che applica iterativamente la funzione step-per.

5.5.1. Funzione stepper

Prima di descrivere l’implementazione dello stepper, descriviamo breve-mente il modo in cui questo procedimento viene svolto. Il risolutore e basatosulla famiglia di metodi nota come Runge-Kutta. Dato il sistema dinamico(45), la soluzione numerica xn+1 ≈ x(t+ dt) viene ottenuta come

xn+1 = xn + dts∑i=1

biki, (54)

con ki i-esimo stadio, s numero di stadi e bi relativi coefficienti. Per otteneregli stadi e necessario risolvere il sistema, nel caso piu generale implicito, datodalle s equazioni

ki = f(xn + dts∑j=1

aijkj, tn + cidt). (55)

Questo tipo di impostazione e comune alla maggior parte ( non tutti ) deimetodi di Runge-Kutta. Frequentemente, per descrivere il metodo, vengonoforniti i coefficienti bi, ci e aij, raggruppati nelle tavole di Butcher

c1 a11 a12 . . . a1s

c2 a21 a22 . . . a2s

......

.... . .

...

cs as1 as2 . . . ass

b1 b2 . . . bs

=c A

b

Esempi di metodi simplettici che rispettano questa classificazione sono imetodi di Gauss-Legendre come il punto medio implicito

12

12

1

31

ed il metodo di Gauss del quarto ordine

12− 1

6

√3 1

414− 1

6

√3

12

+ 16

√3 1

4− 1

6

√3 1

412

12

L’idea dietro all’implementazione e quella, quindi, che l’utente possa forniregli elementi della tavola di Butcher, che vengono poi usati per risolvere ilsistema di equazioni degli stadi, e infine venga fatto uno step in avanti dellasoluzione. Sottolineiamo come il problema degli stadi non implichi necessari-amente la soluzione di un sistema nonlineare: in caso di metodo esplicito,ciascuno stadio richiede solamente di conoscere gli stadi precedenti e quindidobbiamo solamente valutare s volte il termine f . Nel caso di sistema implic-ito abbiamo implementato un semplice metodo di punto fisso per risolvere ilrelativo problema.Utilizzando una strategia gia vista per altre parti della libreria, definiamoun’interfaccia per lo stepper, data da una classe puramente astratta in cui an-diamo ad implementare alcune funzionalita comuni a metodi espliciti ed im-pliciti, i quali poi nelle rispettive classi computeranno gli stadi in maniera di-versa. La classe base e stepperOde, contenuta in include/ODE Solver/stepperOde.hpp,la quale presenta:

• Un costruttore che accetta in input le tavole di Butcher ed il dato in-iziale, usato per inizializzare la soluzione numerica del problema. Nellanostra implementazione abbiamo due membri, Solution e oldSolution,che vengono usati per contenere la soluzione all’istante t e t+dt: per essiverra usata la move semantic, alla base del metodo swap() di Eigen apartire dalla major release 3.3 ed evitare la creazione di oggetti tempo-ranei. Inoltre, anziche salvare un numero arbitrario di stadi, creiamodue matrici aventi lo stesso numero di colonne della soluzione ed unnumero di righe pari al numero di righe della soluzione per il numerodi stadi. Utilizzando la feature delle block expressions data da Eigen,possiamo poi lavorare su sottomatrici senza dover richiedere copie deidati.

• Due metodi initializeStages e resetSolver, i quali servono nel casosi voglia usare lo stesso solver per un problema avente un dato inizialediverso, anche in termini di numero di nodi utilizzati.

32

• Un metodo stageEvaluation, usato per valutare le equazioni (55).

• Un metodo doStep (virtuale) che viene usato per calcolare la soluzionenumerica dopo un time step.

Ciascuna delle due classi derivate stepperOdeExplicit e stepperOdeImplicit,implementa un metodo updateStages per ottenere gli stadi aggiornati. Men-tre nel caso esplicito semplicemente viene chiamato il metodo stageEvaluationdella classe base


2 void stepperOdeExplicit<Container>::updateStages( System& inputSystem,

3 const double Time,

4 const double deltaTime )

5

6 this->stageEvaluation( inputSystem, this->Stages,

7 this->Stages, Time, deltaTime );

8

nel caso implicito il procedimento e solo leggermente piu complesso. Ven-gono definiti come membri aggiuntivi della classe altri due contenitori distadi, tmpStages e oldtmpStages, usati per fare un update e swap dellasoluzione del sistema implicito (55) dopo ogni iterata del metodo del puntofisso.


2 void stepperOdeImplicit<Container>::updateStages( System& inputSystem,



5

6 oldtmpStages.swap(this->Stages);

7 double Error(Precision+1);

8 int Iter(0);

9 do

10 this->stageEvaluation( inputSystem, tmpStages,

11 oldtmpStages, Time, deltaTime );

12 Error = (tmpStages-oldtmpStages).norm();

13 Iter++;

14 oldtmpStages.swap(tmpStages);

15 while( (Error > Precision) && (Iter < MaxIter) );

33

16 this->Stages.swap(tmpStages);

17

Ad ogni time step, come guess iniziale della soluzione del sistema implic-ito degli stadi consideriamo Stages, contenente gli stadi aggiornati al timestep precedente. Ad ogni iterata, l’update degli stadi viene calcolato tramiteil metodo stageEvaluation fino a che l’errore, in termini di quadrato dellanorma Euclidea risulta essere inferiore a Precision o viene raggiunto il nu-mero massimo di iterazioni concesse, evitando in ogni caso copie tramiteswap. Una volta finito questo processo, gli stadi aggiornati vengono mossinel membro Stages.Passiamo quindi ad analizzare quello che succede nella funzione stageEvaluation,per spiegare il motivo per cui abbiamo precedentemente introdotto matrixBlocknella dichiarazione dell’operator () della classe per gestire il sistema dinam-ico.


2 void stepperOde<Container>::stageEvaluation( System& inputSystem,

3 Container& Stage,

4 Container& oldStage,



7

8 Container StagesCombination(numberOfRows,numberOfCols);

9 for( int i(0); i<numberOfStages; i++ )

10

11 StagesCombination = Solution;

12 for( int j(0); j<numberOfStages; j++ )

13 StagesCombination += deltaTime*A(i,j)

14 *oldStage.block(j*numberOfRows,0,

15 numberOfRows,numberOfCols);

16

17 inputSystem( Time+c(i)*deltaTime, StagesCombination,

18 Stage.block(i*numberOfRows,0,numberOfRows,numberOfCols));

19

20

34

Il loop interno serve per andare a calcolare il contributo

xn + dt

s∑j=1

aijkj (56)

che viene usato da inputSystem per valutare f e resituire ki che, per comesono salvati gli stadi, corrisponde ad un blocco della matrice Stage. Sfortu-natamente, il metodo block<> restituisce una expression che non puo essereconvertita in una referenza non costante ad un oggetto di tipo matrice: questovuol dire che se si prova a passare un blocco come referenza, il compilatorefallisce non potendo scrivere il risultato su di un oggetto temporaneo. Invecequello che vogliamo fare e proprio scrivere direttamente sul blocco, senzadover passare per una struttura intermedia che richiederebbe una copia. Unmodo alternativo, compatibile con C++03 e che non richiede la move se-mantic, sarebbe quello di passare il blocco come const reference e poi usareconst cast per modificarlo.

5.5.2. Funzione di integrazione e observer pattern

Una volta definito il modo di fare uno step, l’unico scopo della funzioneintegratrice e quello di iterare la sua chiamata fino alla fine dell’intervallotemporale considerato. L’utente pero, potrebbe essere interessato a moni-torare in qualche modo la soluzione all’avanzare del tempo, plottandola us-ando un’animazione, controllando una funzione della soluzione stessa, comead esempio l’Hamiltoniano o campionando e salvando alcuni snapshots chepoi verranno utilizzati per la generazione di un sistema di basi ridotto. Aquesto proposito, introduciamo il concetto di observer. Anche se quello im-plementato non e propriamente un observer pattern, ne prende ispirazionee utilizza la stessa idea di fondo. L’idea e quella che esista un oggetto carat-terizzato da uno stato e un insieme di observers che reagiscono in qualchemodo ad un cambiamento di quello stato: il reagire, nel nostro caso, significasalvare il dato, calcolare l’Hamiltoniano e/o plottare una superficie nel mo-mento in cui la soluzione e avanzata di un certo numero di iterazioni. Laversione classica di questo pattern consiste in una lista di observers e in unsubject, ovvero l’oggetto osservato: compito del subject e anche quello diavvisare gli observers di un cambiamento di stato, evitando che questi locontrollino ripetutamente.Nella nostra implementazione degli observers sfruttiamo ancora una volta ilpolimorfismo, fornendo una classe astratta che funge da interfaccia, in modo

35

da poter controllare cosa gli observers possano monitorare della soluzione.Essa e definita in include/Observer.hpp con il nome di Observer e rapp-resenta fondamentalmente un functor


2 class Observer

3

4 public:

5 virtual void operator()( Container& Solution,


7 const double deltaTime,

8 const long int Iteration ) = 0;

9 ;

Gli observers implementati fino ad ora sono elencati in Figura 1 Nel

Observer< Container >

ObserverCollector< Container > ObserverPlot< Container >

ObserverCollectorNonlinearity< Container > ObserverPlot1D< Container > ObserverPlot2D< Container > ObserverPlotSymplectic< Container >

ObserverPlot2D_reduced< Container > ObserverPlot2DHamiltonian< Container >

ObserverPlot2DHamiltonian_reduced< Container >

Figure 1: Observers disponibili.

seguito andremo a descriverne brevemente alcuni, partendo da quelli chesi occupano della raccolta dei dati. L’ObserverCollector si occupa dicampionare la soluzione e di salvarla su file alla fine della simulazione perpoterne riutilizzare i dati nel processo di riduzione. Essa viene costruitaricevendo come parametri in input il dato iniziale come referenza ed un in-sieme di parametri legati alla discretizzazione temporale, in modo da poterallocare preventivamente la memoria richiesta per l’intero campionamento.Purtroppo questo non e compatibile con un’estensione tempo adattiva delsolver e quindi andrebbe creato un altro observer in quel caso.Per essere salvati, i dati devono essere riscritti sotto forma di vettore e suc-cessivamente copiati nella relativa colonna della matrice di collezione: dalmomento che i container di Eigen garantiscono continuita di memoria, e pos-sibile saltare il primo passaggio, reinterpretando l’oggetto matrice tramiteuna mappa.


2 void ObserverCollector<Container>::operator()( Container& Solution,

36


4 const double deltaTime,

5 const long int Iteration )

6

7 Map<ColVector_Real> tmpStorage(Solution.data(),Solution.size());

8 if ( Iteration % samplingRate == 0 )

9

10 Storage.col(round(Iteration/samplingRate)) = tmpStorage;

11

12

Una mappa di Eigen richiede espressamente un oggetto derivato dallaclasse MatrixBase di Eigen e la quantita di memoria occupata per esserecostruita: se, come nel nostro caso, come parametro template consideriamoun vettore colonna, la mappa avra quasi tutte le funzionalita di un normalevettore senza dover copiare i dati dell’elemento passato al suo costruttore.Il metodo saveFile, specificando la directory in cui salvare i dati e il nomeche si vuole loro assegnare, salva il file chiamando write binary nella formarichiesta.La classe ObserverCollectorNonlinearity aggiunge alla classe appena de-scritta la possibilita di calcolare gli snapshots delle nonlinearita sfruttando ilnonaffineOperator.Gli observers che si occupano di plottare i dati sfruttano la libreria gnuplot-iostream per interfacciare C++ con gnuplot. Sebbene i contenitori di arraye matrici di librerie come Blitz++ e Armadillo sono supportati, lo stessonon vale per quelli di Eigen: l’unica soluzione e stata quella di convertire idati in vettori di tuple e poi passare quelli a gnuplot. Questo procedimentonon rappresenta un eccessivo costo aggiuntivo, specialmente se comparato albottleneck dato dal modo in cui gnuplot gestisce il plot dei dati. Non es-sendo nativamente legato a C++, gnuplot copia temporanemente files di datisu disco in formato binario e solo successivamente questi vengono plottati:mantenere la struttura del plot (labels, titolo, dominio, etc.) garantisce uncerto speedup ma consigliamo comunque di non richiedere una frequenza diplotting troppo alta per evitare eccessivi rallentamenti.L’ultima classe implementata nel sistema di observers e quella relativa all’analisidell’Hamiltoniano: essa fondamentalmente riceve in ingresso una referenza alfunctor functionUser, il quale viene valutato sulla soluzione ad ogni aggior-namento. Il motivo della referenza e ancora una volta legato al polimorfismo:

37

https://github.com/dstahlke/gnuplot-iostream

https://github.com/dstahlke/gnuplot-iostream

l’utente puo definire una classe derivata in cui viene concretamente imple-mentata la valutazione della quantita da investigare.Nel classico design observer, il subject ha il compito di tenere una lista degliobserver e chiamarli e cancellarli quando necessario. Nel nostro caso sempli-fichiamo questo aspetto sostituendo al subject la funzione integrateODE,facendo coincidere la chiamata e la fine del compito degli observers conlo scope della funzione e chiamandoli attraverso il redifinito operator ().Per passare gli observers alla funzione sfruttiamo i variadic templates, inmodo che il loro numero possa essere non noto a priori. Esistono diversestrategie per fare l’unpack del numero incognito dei parametri, con la piu co-mune e semplice da implementare basata sulla ricorsione della chiamata allafunzione, andando a definire un caso baso in cui nessun parametro viene pas-sato. Nel nostro caso questa strada non e percorribile: avendo un ciclo while

all’interno della funzione, non possiamo riscrivere ricorsivamente questa partedel codice, altrimenti andremmo a risolvere piu volte (ed inutilmente) lostesso problema. Una soluzione viene fornita da std::initializer list edal perfect forwarding scrivendo

1 template<class Container, typename... Fs>

2 void integrateOde(

3 stepperOde<Container>& stepper,

4 DynamicalSystem<Container>& inputSystem,

5 Container& state,

6 double startTime,

7 double endTime,

8 double deltaTime,

9 Fs&&... observer )

10

11 stepper.initializeStages(inputSystem,startTime);

12 long int iteration(0);

13 auto list = (std::forward<Fs>(observer)(stepper.getSolution(),

14 startTime,deltaTime,

15 iteration),0)... ;

16 while( startTime < endTime )

17

18 stepper.doStep(inputSystem,startTime,deltaTime);

19 startTime += deltaTime;

20 ++iteration;

38

21 list = (std::forward<Fs>(observer)(stepper.getSolution(),

22 startTime,deltaTime,

23 iteration),0)... ;

24

25 ;

In C++17, e possibile evitare l’introduzione di initializer list edusare le piu recenti fold expressions.

5.6. Riduzione di modello

Una volta che i dati della simulazione completa sono stati raccolti e salvaticon un certo identificativo, e possibile costruire il modello ridotto. Restrin-gendosi al caso di riduzione di modello tramite basi ridotte con trattamentodella nonlinearita simil-DEIM (o comunque basato su interpolazione) e pos-sibile definire alcune funzionalita che l’oggeto adibito a questo scopo deveavere. Esso deve essere in grado, innanzitutto, di assemblare la matrice deglisnapshots a partire da datafile relativi allo stesso problema e generare unsistema di basi ridotto. Siccome questa procedura puo, come nel nostrocaso, essere basata su SVD, risulta conveniente poter salvare e leggere dafile una base gia creata, dato il costo computazionale che pio avere una de-composizione. Inoltre, deve poter essere capace di, dato un problema con unelevato numero di incognite, generare il relativo sistema ridotto, in manierache esso abbia la stessa forma di sistema dinamico che viene numericamenterisolto tramite la procedura precedentemente descritta.Queste operazioni sono date in parte nella classe virtuale modelReduction

contenuta in model reduction.hpp. Alcune di esse, trasversali al metodo diriduzione scelto, sono state definite per questa classe, mentre la definizionedi altre viene lasciata a classi derivate. Il costruttore di questa classe ha ilcompito principale di identificare e fare una lista di tutti i file di dati salvati,in modo da sapere quali raccogliere successivamente per generare la matricedegli snapshots.


2 modelReduction<Container>::modelReduction(

3 string dataDirectory,

4 string equationName,

5 int dimension,

6 long int numberOfNodes,

7 int numberOfSnapshots )

39

8

9

10 DIR *dp;

11 struct dirent *dirp;

12 fileSolution = equationName + "_" + to_string(dimension) + "D_";

13 string fileProblemIdentifierComplete = fileSolution

14 + to_string(numberOfNodes)

15 + "_nodes_"

16 + to_string(numberOfSnapshots)

17 + "_snapshots_";

18 if ((dp = opendir(dataDirectory.c_str())) != NULL)

19 while ((dirp = readdir(dp)) != NULL)

20 string tmpFile = dirp->d_name;

21 auto res = mismatch(fileProblemIdentifierComplete.begin(),

22 fileProblemIdentifierComplete.end(),

23 tmpFile.begin());

24 if ( res.first == fileProblemIdentifierComplete.end() )

25

26 fileList.push_back(dataDirectory+dirp->d_name);

27

28

29 closedir (dp);

30 else

31 perror ("The directory is not in the right place ");

32

33

Per questa operazione si e resa necessaria la libreria boost::filesystem

per alcune sue funzionalita. La riduzione del generico sistema di partenza

dx

dt= Ax(t) + g(x(t), t) (57)

segue due strade diverse per quanto riguarda parte affine e non. Mentre per laparte affine essa viene fatta internamente alla classe modelReduction tramiteil metodo reduceAffineOpt, per la parte nonaffine definiamo una classe spe-ciale per la classe dei termini nonlineari, ovvero nonaffineOperatorReduced.La classe in questione riceve come parametri in input un oggetto della classenonaffineOperator per referenza e gli elementi necessari per la sua riduzione,

40

ovvero la matrice interpolante e la base interpolata: la referenza e necessariaancora per sfruttare il polimorfismo, dal momento che modelReduction rap-presenta solo un’interfaccia per la definizione del termine nonlineare. Questaimplementazione fa si che la classe adibita alla gestione di sistemi linearipossare servire per definire sia grandi sistemi sparsi ottenuti tramite metodicome differenze o elementi finiti e sistemi di dimensione molto contenuta ri-cavati dalla riduzione di modello.La parte fondamentale della riduzione viene svolta nelle classi CL e CL SDEIM,dove vengono computate le basi ridotte per problemi simplettici. Nel seguitoelenchiamo solamente le features fondamentali della seconda classe, adibitaal caso nonlineare, dal momento che la prima ne rappresenta un caso sem-plificato. La prima operazione che deve essere compiuta, una volta identi-ficati i files contenenti i dati da utilizzare per la riduzione, e quella di an-dare ad assemblare la matrice degli snapshots, operazione svolta dal metodoassembleSnapshotMatrix. Ricordando che la soluzione di un sistema di-namico per un valore del parametro αi viene salvata nel formato[

qαit1 qαit2 . . . qαitspαit1 pαit2 . . . pαits

]=

[qαi

pαi

](58)

con s numero di snapshots raccolti. Compito di assembleSnapshotMatrixe quindi quello di costruire la matrice degli snapshots

S = S ∪[qαi pαi g(qαi) g(pαi)

]con i = 1, . . . ,m cardinalita dell’insieme dei parametri. Una propria basesimplettica e poi generata tramite generateBasis, a cui segue l’implementazionedella variante simplettica del DEIM descritta nell’apposito paragrafo.Una nota a margine riguarda l’estensione degli observers al caso ridotto. Lasoluzione del generico problema ridotto (41) e data dai coefficienti rispettoalla base ridotta, non rispetto ai nodi scelti per discretizzare il domino.Questo rende necessario l’introduzione di ulteriori classi derivate, specificatedall’estensione reduced, che, accentando come parametro di costruzionela base U, ricostruiscano la soluzione originaria per poi chiamare il rela-tivo metodo definito nella classe parente e quindi plottare i dati, calcolarel’Hamiltoniano o salvare la soluzione.

6. Test numerici

Stabilita e accuratezza del metodo di riduzione considerato sono illus-trate tramite simulazioni numeriche di varianti delle equazioni delle onde.

41

Dopo aver definito l’Hamiltoniano per questa classe di equazioni, analizzer-emo prima il caso lineare con il tempo come unico parametro per testare ilmetodo in un caso semplice per poi spostarci ad un caso nonlineare con piuparametri coinvolti. Il motivo per cui le equazioni delle onde rappresentanoun valido test per un metodo di riduzione e dato dalla loro natura iperbol-ica: questo impedisce di usare alcune delle stime, basate sulla coercivita, chenormalmente i classici metodi di riduzione usano per garantire la stabilitadel problema in casi parabolici. In particolar modo, i particolari problemiconsiderati consistono in estensioni al caso bidimensionale dei test condottiin [19], dove si dimostra che la combinazione standard POD - DEIM non im-plica stabilita del problema ridotto, anche nel semplice caso lineare. Quandoun elevato numero di parametri influenza la soluzione, una valida alternativaal metodo considerato in questo lavoro viene fornito in [21]: sfortunatamentequesto porta ad una riduzione dell’accuratezza. In caso di uno spazio diparametri di dimensione limitata, una riduzione basata su SVD come quellacondotta dal Cotangent Lift permette di mantenere l’accuratezza dei metodiclassici e allo stesso tempo forzare la stabilita della soluzione.

6.1. Equazione lineare delle onde

Sia u = u(x, t). L’equazione delle onde lineare con costante di movimentoc e data da

∂2u

∂t2= c2∆u. (59)

nel dominio Ω = [0, 2π] × [0, 2π]. Considerando le posizioni generalizzateq = u e i momenti generalizzati p = ∂u

∂t, insieme a condizioni periodiche al

dominio, l’equazione Hamiltoniana associata a (59) e data da∂q∂t

= δHδp,

∂p∂t

= − δHδq,

(60)

dove l’Hamiltoniano e e

H(q, p) =

∫Ω

[1

2p2 +

1

2c2q2

x

]dx. (61)

E possibile dimostrare [19] che una discretizzazione a differenze finite constencil di dimensione arbitraria mantiene la struttura simplettica del prob-lema anche nella sua versione semi-discreta se viene impiegato uno schema

42

centrato. Dato n = nx × ny il numero totale di nodi considerati, definiamocome y ∈ R2n la soluzione del sistema di ODE

dy

dt= Ay =

[0n Inc2D 0n

] [qp

], (62)

con D operatore rappresentante il Laplaciano discreto. Nel nostro caso an-dremo a considerare uno schema a differenze finite caratterizzato da 5 puntiper stencil per ogni dimensione (9 punti per nodo). Il numero totale di nodin e pari a 40000, 200 per lato, e quindi la dimensione del problema e ugualea 2n = 80000. Da un punto di vista implementativo, putroppo, la libre-ria Eigen non permette di definire la matrice A tramite la versatile classeblock<>: il motivo e che, essendo A sparsa, identificare un blocco tramiteindici risulterebbe poco pratico ed efficiente. Lo stratagemma utilizzato estato quindi quello di giocare con gli orientamenti di una matrice sparsa:matrici orientate per righe (colonne) sono naturalemnte concatenabili perrighe (colonne) e questo viene usato per assemblare A pezzo per pezzo.L’integrazione di (62) viene condotta tramite il punto medio implicito, us-ando un time step uguale a dt = 0.001s per una simulazione di durata 5s:il relativo sistema implicito viene risolto con metodo iterativo che raggiungeconvergenza quando la norma l2 della differenza tra due iterate successive eminore di 10−9. Consideriamo come condizione iniziale

y0 = y(0) =

[q0

p0

]=

[4atan

(e− |x−π|+|y−π|√

0.75

)0

](63)

che origina un sistema periodico che genera un’onda che si propaga in manieracircolare ed interagisce con il bordo del dominio.Fissiamo la costante di moto c a 1 e la relativa riduzione coinvolgera sola-mente il parametro temporale. Raccogliamo uno snapshot ogni dt = 0.02s,tenendo quindi un rateo di campionamento pari a 120 rispetto al passo tem-porale di simulazione. La matrice degli snapshots avra quindi una dimensionepari a 80000×250. Una riduzione di questo tipo risulta utile quando il mod-ello originale e particolarmente oneroso da risolvere e il modello ridotto vieneusato per investigare le dinamiche per time steps molto piccoli.Il primo indicatore che andiamo ad osservare per valutare la qualita dellariduzione e il tempo richiesto per completare una simulazione ridotta com-parata al caso completo. Sulla macchina considerata per i test (Intel(R)Core(TM) i7-6700HQ CPU @ 2.60GHz) il problema originale viene risolto

43

in 26s considerando come unico observer quello che salva periodicamente lasoluzione: in Figura (2) sono riportati invece i tempi richiesti dal modelloridotto in funzione del numero di basi. Il guadagno risulta essere evidente:anche nel caso in cui 100 basi sono considerate, il modello ridotto e all’incirca200 volte piu veloce. Per capire se la riduzione risulta essere effettivamenteaccurata, consideriamo il variare delle media temporale della norma l2 delladifferenza tra soluzione completa e ridotta: questo risultato viene riportatoin Figura (3), sempre in funzione del numero di basi. Da questo notiamoche l’errore incomincia a calare sensibilmente a partire dalle 80 basi, valoreper cui l’approssimmazione introdotta dalla riduzione risulta piu che accetta-bile. Per trasmettere questa conclusione in termini qualitativi invitiamo adanalizzare alcuni snapshots delle soluzioni numeriche riportate in Figura (4):con 30 basi, nonostante il modello ridotto riesca a seguire l’andamento delmodello completo, si notano comunque delle discrepanze, mentre a 100 basie difficile riuscire a distinguere la riduzione dall’originale.

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

10 20 30 40 50 60 70 80 90 100

Tim

e [

s]

Number of basis []

Time

Figure 2: Tempo richiesto per risolvere il sistema ridotto in funzione del numero di basinel caso di equazione delle onde lineare. Viene considerato anche il tempo richiesto persalvare la soluzione campionata.

6.2. Equazione di Sine-Gordon

E possibile generalizzare la definizione di equazione delle onde andandoa comprendere un contributo nonlineare. La relativa PDE che otteniamo hala forma

∂2u

∂t2= c2∆u+ g(u). (64)

44

0.001

0.01

0.1

1

10 20 30 40 50 60 70 80 90 100

Avera

ge d

iscr

ete

l2

err

or

Number of basis []

Error

Figure 3: Andamento dell’errore in funzione del numero di basi per ridurre l’equazionedelle onde lineare. L’errore preso in considerazione e la media sull’intervallo temporaledella norma l2 discreta della differenza tra la soluzione del modello completo e la soluzionedel modello ridotto.

dove u ha lo stesso significato del caso lineare, ovvero displacement rispettoad un valore medio. In questo all’Hamiltoniano si aggiunge un contributodipendente dal nuovo termine nonlineare

H(q, p) =

∫Ω

[1

2p2 +

1

2c2q2

x +G(q)

]dx, (65)

con G′(q) = g(q). Consideriamo in particolare il caso in cui G(u) = 1−cos(u)e g(u) = sin(u) mantenendo le stesse impostazioni in termini di discretiz-zazione spaziale ed integrazione numerica: la stessa condizione iniziale vieneconsiderata, dando origine ad una soluzione nota onda solitaria periodica. Ilproblema in questione prende il nome di equazione di Sine-Gordon e vieneutilizzato per studi di campi relativistici e giunzioni di Josephson. Per ren-dere piu complessa l’operazione di riduzione, consideriamo il parametro cnon piu come fisso ma variabile nell’intervallo [0.5, 1.5]: questo influenza lavelocita di propagazione dell’onda, andando fisicamente a modellizzare la per-meabilita del materiale di propagazione. Questo intervallo viene discretizzatousando i punti 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, e le relative 6 simulazioni sono sal-vate nei corrispondenti datafile, ognuno del peso approssimativo di 300MB.La risultante matrice degli snapshots, creata durante la fase di riduzione, hadimensione 40000 × 6024, mettendo a dura prova la SVD fornita da Eigenche per decomporla impiega all’incirca 10 minuti sulla macchina test: il mo-

45

(a) Completo, t = 0s (b) 30 basi, t = 0s (c) 100 basi, t = 0s

(d) Completo, t = 2.5s (e) 30 basi, t = 2.5s (f) 100 basi, t = 2.5s

(g) Completo, t = 5s (h) 30 basi, t = 5s (i) 100 basi, t = 5s

Figure 4: Comparazione qualitativa tra diversi snapshots della soluzione del modello com-pleto e della soluzione del modello ridotto tramite Cotangent Lift per equazione delle ondelineare.

tivo e che il modulo SVD implementato in Eigen e stato supportato solo apartire da Eigen 3.3, e l’attuale algoritmo piu efficiente disponibile risultainferiore a quanto fornito da software concorrenti quali Armadillo, Python e

46

Matlab. Una volta ottenuta la base ridotta abbiamo considerato un valoredel parametro test, nel nostro caso c = 1, che non fosse compreso all’internodel gruppo di simulazioni archiviate, e abbiamo condotto la stessa analisiconsiderata per il caso lineare ( Figure 5, 6 e 7 ) ottenendo simili conclusionicirca l’errore di riduzione, che descresce in maniera considerevole a partireda circa 80 basi. Riguardo al tempo computazionale, abbiamo ancora che ilmodello risulta essere decisamente piu veloce del modello completo (39s) macon rapporti piu contenuti rispetto al caso lineare. Un’ulteriore analisi con-dotta in questo caso riguarda l’andamento dell’Hamiltoniano. In Figura 8,l’andamento della funzione hamiltoniana viene riportato nel caso del sistemaoriginale e quando vengono impiegate 40, 80 e 120 basi. La presenza di oscil-lazioni rispetto ad un valore medio e dovuta principalmente all’utilizzo di unmetodo interpolatorio nel trattamento della nonlinearita e diminuisce, comeci aspettiamo,all’aumentare delle basi. Un altro dato interessante, particolar-mente evidente nel caso 40 basi, riguarda il valore medio dell’Hamiltoniano,che risulta essere diverso rispetto al caso completo e converge ad esso piubasi vengono utilizzate. Questo e completamente in accordo con la teoriasottostante al metodo: il problema ridotto presenta, espresso in termini delmodello completo e delle basi ridotte, la forma

d

dty = J2kH(Ay) (66)

con A matrice simplettica. L’Hamiltoniano del sistema (66) risulta quindi

essere H(·) = H(A·), e nonH, ovvero l’Hamiltoniano del sistema di partenza.

47

0

0.2

0.4

0.6

0.8

1

1.2

20 40 60 80 100 120

Tim

e [

s]

Number of basis []

Time

Figure 5: Tempo richiesto per risolvere il sistema ridotto in funzione del numero di basi nelcaso di equazione nonlineare di Sine Gordon. Viene considerato anche il tempo richiestoper salvare la soluzione campionata.

0.001

0.01

0.1

1

20 40 60 80 100 120

Avera

ge d

iscr

ete

l2

err

or

Number of basis []

Error

Figure 6: Andamento dell’errore in funzione del numero di basi per ridurre l’equazionenonlineare di Sine Gordon. L’errore preso in considerazione e la media sull’intervallotemporale della norma l2 discreta della differenza tra la soluzione del modello completo ela soluzione del modello ridotto.

7. Conclusioni e futuri miglioramenti

La libreria implementata in C++ ha dimostrato di saper gestire in manieraefficiente la riduzione di modelli dinamici, simplettici e non, generando mod-elli ridotti integrabili in piccole frazioni del tempo richiesto alla soluzionedel modello completo, con speedup dell’ordine delle centinaia nel caso di

48

(a) Completo, t = 0s (b) 30 basi, t = 0s (c) 100 basi, t = 0s

(d) Completo, t = 2.5s (e) 30 basi, t = 2.5s (f) 100 basi, t = 2.5s

(g) Completo, t = 5s (h) 30 basi, t = 5s (i) 100 basi, t = 5s

Figure 7: Comparazione qualitativa tra diversi snapshots della soluzione del modello com-pleto e della soluzione del modello ridotto tramite Cotangent Lift per equazione nonlinearedi Sine Gordon.

riduzione lineare e decine nel caso nonaffine, con soluzioni stabili. Nonostantevari siano i metodi implementati per quanto riguarda l’integrazione numericae la riduzione di modello, la struttura secondo cui il codice e stato scritto si

49

-24.194966

-24.194966

-24.194965

-24.194965

-24.194965

-24.194965

-24.194965

0 1 2 3 4 5

Ham

ilto

nia

n

Complete-24.9

-24.8

-24.7

-24.6

-24.5

-24.4

-24.3

-24.2

-24.1

0 1 2 3 4 5

40 basis

-24.36-24.34-24.32

-24.3-24.28-24.26-24.24-24.22

-24.2-24.18

0 1 2 3 4 5

Ham

ilto

nia

n

Time[s]

80 basis-24.245

-24.24-24.235

-24.23-24.225

-24.22-24.215

-24.21-24.205

-24.2-24.195

-24.19

0 1 2 3 4 5

Time[s]

120 basis

Figure 8: Comparazione tra l’Hamiltoniano nel caso completo e le approssimazioni usandobasi ridotte.

presta a numerose estensioni. In primis, nonostante la descrizione tramitetavole di Butcher dei metodi Runge-Kutta sia comune ad un largo gruppodi metodi, non li riguarda propriamente tutti. Il famoso metodo Stromer-Verlet, usato frequentemente per calcoli astronomici e riguardanti reazionichimiche, di natura simplettica non rientra in questa categoria ma potrebbeessere facilmente definito tramite classe derivata da odeSolver. La stessatavola di Butcher potrebbe essere ampliata in maniera tale da implementareun metodo embedded, di un ordine inferiore, richiedendo un semplice vettoreb∗ aggiuntivo per introdurre un time step adattivo.Molto piu interessanti sono le possibili estensioni in termini di riduzione dimodello: data la struttura fornita in modelReduction, simili metodi possonoessere implementati in maniera da fornire un benchmark di confronto. Altripossibili metodi sono il complexSVD e il metodo basato su programmazionenonlineare descritti in [19] ed entrambi simplettici, oppure il dynamical mode

50

decomposition (DMD), recentemente introdotto per ridurre sistemi parabol-ici e le cui prestazioni nei confronti di sistemi iperbolici ed Hamiltoniani sonoancora sotto studio.

51

References

[1] M. Bachmayr, A. Cohen, Kolmogorov widths and low-rank approxima-tions of parametric elliptic PDEs * (2015).

[2] C. Rowley, T. Colonius, R. Murray, Model reduction for compressibleflows using POD and Galerkin projection, Physica D: Nonlinear Phe-nomena (2004).

[3] A. Quarteroni, G. Rozza (Eds.), Reduced Order Methods for Model-ing and Computational Reduction, Springer International Publishing,POD1, 2014.

[4] J. Hesthaven, G. Rozza, B. Stamm, Certified reduced basis methods forparametrized partial differential equations (2016).

[5] Y. Liang, H. Lee, S. Lim, W. Lin, K. Lee, Proper orthogonal decompo-sition and its applicationsPart I: Theory, Journal of Sound and (2002).

[6] G. Berkooz, P. Holmes, J. L. Lumley, THE PROPER ORTHOGONALDECOMPOSITION IN THE ANALYSIS OF TURBULENT FLOWS,Annu. Rev. Fluid Mech 25 (1993) 539–75.

[7] A. Rosenfeld, Digital picture processing (1976).

[8] S. Shvartsman, C. Theodoropoulos, Order reduction for nonlinear dy-namic models of distributed reacting systems, Journal of Process (2000).

[9] S. Shvartsman, I. Kevrekidis, Low-dimensional approximation and con-trol of periodic solutions in spatially extended systems, Physical ReviewE (1998).

[10] M. Graham, I. Kevrekidis, Alternative approaches to the Karhunen-Loeve decomposition for model reduction and data analysis, Computers& chemical engineering (1996).

[11] A. J. Siade, M. Putti, W. W.-G. Yeh, Snapshot selection for ground-water model reduction using proper orthogonal decomposition, WaterResources Research 46 (2010) 1–13.

[12] R. Everson, L. Sirovich, Karhunen–Loeve procedure for gappy data, J.Opt. Soc. Am. A 12 (1995) 1657–1664.

52

[13] M. Barrault, Y. Maday, N. C. Nguyen, A. T. Patera, An empiricalinterpolation’ method: application to efficient reduced-basis discretiza-tion of partial differential equations, Comptes Rendus Mathematique339 (2004) 667–672.

[14] P. Astrid, S. Weiland, K. Willcox, T. Backx, Missing Point Estima-tion in Models Described by Proper Orthogonal Decomposition, IEEETransactions on Automatic Control 53 (2008) 2237–2251.

[15] S. Chaturantabut, D. C. Sorensen, Nonlinear Model Reduction via Dis-crete Empirical Interpolation, SIAM Journal on Scientific Computing32 (2010) 2737–2764.

[16] K. Carlberg, C. Farhat, J. Cortial, D. Amsallem, The GNAT methodfor nonlinear model reduction: Effective implementation and applica-tion to computational fluid dynamics and turbulent flows, Journal ofComputational Physics 242 (2013) 623–647.

[17] M. a. Grepl, Y. Maday, N. C. Nguyen, A. T. Patera, Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations,ESAIM: Mathematical Modelling and Numerical Analysis 41 (2007)575–605.

[18] A. Manzoni, A. Quarteroni, G. Rozza, Shape optimization for viscousflows by reduced basis methods and free-form deformation, InternationalJournal for Numerical Methods in Fluids 70 (2012) 646–670.

[19] L. Peng, K. Mohseni, Symplectic model reduction of Hamiltonian sys-tems, SIAM Journal on Scientific Computing (2016).

[20] S. Lall, P. Krysl, J. E. Marsden, Structure-preserving model reductionfor mechanical systems, Physica D 184 (2003) 304–318.

[21] B. M. Afkham, J. S. Hesthaven, Structure Preserving Model Reductionof Parametric Hamiltonian Systems (2017).

[22] E. E. Hairer, C. Lubich, G. Wanner, Geometric numerical integration: structure-preserving algorithms for ordinary differential equations,Springer, 2006.

53

[23] H. Z. Hassan, A. A. Mohamad, G. E. Atteia, An algorithm for the finitedifference approximation of derivatives with arbitrary degree and orderof accuracy, Journal of Computational and Applied Mathematics 236(2012) 2622–2631.

54

Libreria C++ per la riduzione di modello di problemi ...forma/Didattica/ProgettiPacs/... ·...

Documents

Transcript of Libreria C++ per la riduzione di modello di problemi ...forma/Didattica/ProgettiPacs/... ·...