Integrazione geometrica numericaper problemi stocastici di evoluzione

Raffaele D’AmbrosioDISIM - Universita dell’Aquila

Congresso GNCSMontecatini, 11–13 febbraio 2020

Progetto GNCS 2018:“Approssimazione numerica diproblemi di evoluzione: aspettideterministici e stocastici”;

Progetto GNCS 2019:“Problemi di evoluzione e lorodiscretizzazione: questioni distabilita lineare e non lineare”.

Aspetti numerici nei problemi evolutivi (e.g., problemistocastici, sistemi discontinui, ODEs matriciali, equazionicon ritardo, PDEs che generano fronti d’onda periodici, ...);

vari inserti di algebra lineare numerica, ottimizzazionenumerica, teoria dei sistemi dinamici;

numerose applicazioni (e.g., oscillatori chimici accoppiati,dinamica delle popolazioni, matrix completion nei sistemi diraccomandazione, biomatematica, teoria dei grafi, ...).

Partecipanti non strutturati:Alessia Ando (Udine, 2018 e 2019), Asma Farooq (Trieste 2019), Da-

vide Liessi (Udine, 2018), Martina Moccaldi (Salerno, 2018), Francesca

Scarabel (Helsinki, 2018 e 2019), Simone Spada (Trieste, 2018 e 2019).

Pubblicazioni nel biennio 2018–2019

Andreotti, E., Edelmann, D., Guglielmi, N., Lubich, C. Constrained graph partitioning viamatrix differential equations, SIAM Journal on Matrix Analysis and Applications 40(1),1–22 (2019).

Banjai, L., Lopez-Fernandez, M. Efficient high order algorithms for fractional integrals andfractional differential equations, Numerische Mathematik 141(2), 289–317 (2019).

Boccarelli, A., Esposito, F., Coluccia, M., Frassanito M.A., Vacca, A., Del Buono, N.Improving knowledge on the activation of bone marrow fibroblasts in MGUS and MMdisease through the automatic extraction of genes via a nonnegative matrix factorizationapproach on gene expression profiles, Journal of Translational Medicine 16(1),217 (2018).

Bohinc, K., Rescic, J., Spada, S., May, S., Maset, S. Influence of added salt on the surfaceinduced ordering of nanoparticles with discretely distributed charges, Journal of MolecularLiquids 294,111134 (2019).

Breda, D., Liessi, D. Approximation of eigenvalues of evolution operators for linear renewalequations, SIAM Journal on Numerical Analysis 56(3), 1456–1481 (2018).

Breda, D., Menegon, G., Nonino, M. Delay equations and characteristic roots: Stability andmore from a single curve, Electronic Journal of Qualitative Theory of Differential Equations,89 (2018).

Cesarone, F., Pepe, P., Guglielmi, N. Solution by sampled-data control of a consensusproblem: An approach by stabilization in the sample-and-hold sense, Proceedings of theAmerican Control Conference 8430969, 6493–6498 (2018).

Pubblicazioni nel biennio 2018–2019

Cardone, A., Conte, D., D’Ambrosio, R., Paternoster, B. Stability issues for selectedstochastic evolutionary problems: A review, Axioms 7(4),91 (2018).

Cardone, A., Conte, D., D’Ambrosio, R., Paternoster, B. Collocation methods for volterraintegral and integro-differential equations: A review, Axioms 7(3),45 (2018).

Cardone, A., Conte, D., D’Ambrosio, R., Paternoster, B. On quadrature formulas foroscillatory evolutionary problems, International Journal of Circuits, Systems and SignalProcessing 12, 58–64 (2018).

Cardone, A., Conte, D., Paternoster, B. Two-step collocation methods for fractionaldifferential equations, Discrete and Continuous Dynamical Systems - Series B 23(7),2709–2725 (2018).

Cardone, A., D’Ambrosio, R., Paternoster, B. A spectral method for stochastic fractionaldifferential equations, Applied Numerical Mathematics 139, 115–119 (2019).

Casalino, G., Castiello, C., Del Buono, N., Mencar, C. A framework for intelligent Twitterdata analysis with non-negative matrix factorization, International Journal of WebInformation Systems 14(3), 334–356 (2018).

Cicone, A., Guglielmi, N., Protasov, V.Y. Linear switched dynamical systems on graphs,Nonlinear Analysis: Hybrid Systems 29, 165–186 (2018).

Citro, V., D’Ambrosio, R. Long-term analysis of stochastic θ-methods for damped stochasticoscillators, Applied Numerical Mathematics (2019).

Pubblicazioni nel biennio 2018–2019

Colombo, A., Del Buono, N., Lopez, L., Pugliese, A. Computational techniques to locatecrossing/sliding regions and their sets of attraction in non-smooth dynamical systems,Discrete and Continuous Dynamical Systems - Series B 23(7), 2911–2934 (2018).

Conte, D., D’Ambrosio, R., Moccaldi, M., Paternoster, B. Adapted explicit two-step peermethods, Journal of Numerical Mathematics 27(2), 69–83 (2019).

Conte, D., D’Ambrosio, R., Paternoster, B. On the stability of ϑ-methods for stochasticvolterra integral equations, Discrete and Continuous Dynamical Systems - Series B 23(7),2695–2708 (2018).

Conte, D., Paternoster, B., Moradi, L., Mohammadi, F. Construction of exponentially fittedexplicit peer methods, International Journal of Circuits, Systems and Signal Processing 13,501–506 (2019).

D’Ambrosio, R., Moccaldi, M., Paternoster, B., Rossi, F. Adapted numerical modelling ofthe Belousov–Zhabotinsky reaction, Journal of Mathematical Chemistry 56(10), 2876–2897(2018).

D’Ambrosio, R., Moccaldi, M., Paternoster, B., Rossi, F. Stochastic numerical models ofoscillatory phenomena, Communications in Computer and Information Science 830, 59–69(2018).

D’Ambrosio, R., Moccaldi, M., Paternoster, B. Numerical preservation of long-termdynamics by stochastic two-step methods, Discrete and Continuous Dynamical Systems -Series B 23(7), 2763–2773 (2018).

Pubblicazioni nel biennio 2018–2019

D’Ambrosio, R., Moccaldi, M., Paternoster, B. Parameter estimation inIMEX-trigonometrically fitted methods for the numerical solution of reaction–diffusionproblems, Computer Physics Communications 226, 55–66 (2018).

D’Ambrosio, R., Paternoster, B. Multivalue collocation methods free from order reduction,Journal of Computational and Applied Mathematics, in press (2019).

Del Buono, N., Elia, C., Garrappa, R., Pugliese, A. Preface: “structural dynamical systems:Computational aspects”, Discrete and Continuous Dynamical Systems - Series B 23(7), i(2018).

Dieci, L., Eirola, T., Elia, C. Periodic orbits of planar discontinuous system underdiscretization, Discrete and Continuous Dynamical Systems - Series B 23(7), 2743–2762(2018).

Dieci, L., Elia, C. Smooth to discontinuous systems: A geometric and numerical method forslow-fast dynamics, Discrete and Continuous Dynamical Systems - Series B 23(7), 2935–2950(2018).

Dieci, L., Papini, A., Pugliese, A. Coalescing points for eigenvalues of banded matricesdepending on parameters with application to banded random matrix functions, NumericalAlgorithms 80(4), 1241–1266 (2019).

Elia, C., Maroto, I., Nunez, C., Obaya, R. Existence of global attractor for anonautonomous state-dependent delay differential equation of neuronal type,Communications in Nonlinear Science and Numerical Simulation 78,104874 (2019).

Pubblicazioni nel biennio 2018–2019

Esposito, F., Gillis, N., Del Buono, N. Orthogonal joint sparse NMF for microarray dataanalysis, Journal of Mathematical Biology 79(1), 223–247 (2019).

Fazzi, A., Guglielmi, N., Markovsky, I. An ODE-based method for computing theapproximate greatest common divisor of polynomials, Numerical Algorithms 81(2), 719–740(2019).

Guglielmi, N., Manetta, M. Stability of gyroscopic systems with respect to perturbations,Springer INdAM Series 30, 253–266 (2019).

Guglielmi, N., Mason, O., Wirth, F. Barabanov norms, Lipschitz continuity andmonotonicity for the max algebraic joint spectral radius, Linear Algebra and ItsApplications 550, 37–58 (2018).

Guglielmi, N., Protasov, V.Y. On the closest stable/unstable nonnegative matrix and relatedstability radii, SIAM Journal on Matrix Analysis and Applications 39(4), 1642–1669 (2018).

Guglielmi, N., Scalone, C. Computing the closest real normal matrix and normalcompletion, Advances in Computational Mathematics, in press (2019).

Gyllenberg, M., Scarabel, F., Vermiglio, R. Equations with infinite delay: Numericalbifurcation analysis via pseudospectral discretization, Applied Mathematics andComputation 333, 490–505 (2018).

Lopez, L., Maset, S. Time-transformations for the event location in discontinuous ODEs,Mathematics of Computation 87(313), 2321–2341 (2018).

Pubblicazioni nel biennio 2018–2019

Markovsky, I., Fazzi, A., Guglielmi, N. Applications of polynomial common factorcomputation in signal processing, Lecture Notes in Computer Science (including subseriesLecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 10891, 99–106(2018).

Maset, S. Conditioning and relative error propagation in linear autonomous ordinarydifferential equations, Discrete and Continuous Dynamical Systems - Series B 23(7),2879–2909 (2018).

Sadeghpour, M., Breda, D., Orosz, G. Stability of Linear Continuous-Time Systems withStochastically Switching Delays, IEEE Transactions on Automatic Control 64(11), 8665935,4741–4747 (2019).

Scalone, C., Guglielmi, N. A gradient system for low rank matrix completion, Axioms 7(3),51 (2018).

Spada, S., Maset, S., Bohinc, K. Interaction between like-charged surfaces mediated byuniformly charged counter-nanoparticles, International Journal of Modern Physics B33(10),1950092 (2019).

Citro, V., D’Ambrosio, R., Di Giovacchino, S. A-stability preserving perturbation ofRunge-Kutta methods for stochastic differential equations, Appl. Math. Lett. (2019).

Citro, V., D’Ambrosio, R. Nearly conservative multivalue methods with extended boundedparasitism, Appl. Numer. Math. (2019).

Discretizzazioni structure-preserving

Problemi deterministici

Integrazione geometrico-numerica di problemi Hamiltonianicon E. Hairer (Ginevra); J. Butcher (Auckland);

sistemi dinamici discontinuicon L. Dieci (Georgia Tech), Fulbright; A. Scotti (Milano);

PDEs che generano fronti d’onda periodicicon B. Paternoster (Salerno).

Problemi stocastici

Hamiltonianicon D. Cohen (Umea);

non lineari (contrattivita in media, oscillatori smorzati)con E. Buckwar (Linz); S. Di Giovacchino, C. Scalone (L’Aquila);

con memoriacon D. Conte e B. Paternoster (Salerno).

Integrazione geometrica (deterministica)

Problemi Hamiltoniani

y(t) = J ∇H(y(t)), t ≥ 0,

y(0) = y0,J =

[0 −II 0

].

y(t) =

[p(t)

q(t)

]∈ R2d, momenti e coordinate generalizzate;

la funzione Hamiltoniana

H : R2d → R

e un integrale primo del problema, dunque

H(y(t)) = H(y0), t ≥ 0.

Infatti,

d

dtH(y(t)) = ∇H(y(t))Ty(t) = ∇H(y(t))TJ−1∇H(y(t)) = 0.

Integratori simplettici

I metodi Runge-Kutta simplettici conservano esattamenteHamiltoniane quadratiche;

se la soluzione numerica calcolata con un metodo Runge-Kutta diordine p giace in un insieme compatto, allora

H(yn) = H(y0) +O(hp),

per tempi lunghi (Benettin-Giorgilli, 1994);

metodi lineari multistep quasi-conservativia lungo termine (Hairer-Lubich, 2004);

metodi multivalue quasi conservativi a lungo termine (Butcher,D’Ambrosio, 2017; D’Ambrosio, Hairer, 2013, 2014).

Esempio: pendolo semplice

H(y) =p2

2− cos q.

Integratori simplettici

I metodi Runge-Kutta simplettici conservano esattamenteHamiltoniane quadratiche;

se la soluzione numerica calcolata con un metodo Runge-Kutta diordine p giace in un insieme compatto, allora

H(yn) = H(y0) +O(hp),

per tempi lunghi (Benettin-Giorgilli, 1994);

metodi lineari multistep quasi-conservativia lungo termine (Hairer-Lubich, 2004);

metodi multivalue quasi conservativi a lungo termine (Butcher,D’Ambrosio, 2017; D’Ambrosio, Hairer, 2013, 2014).

Esempio: pendolo semplice

H(y) =p2

2− cos q.

Pendolo semplice: spazio delle fasi

Runge-Kutta di Gauss a 2 stadi (simplettico)

Pendolo semplice: spazio delle fasi (ctd.)

Eulero esplicito (non simplettico)

Uri M. Ascher, Surprising computations, Appl. Numer. Math. (2012).

Pendolo semplice: spazio delle fasi (ctd.)

Eulero esplicito (non simplettico)

Uri M. Ascher, Surprising computations, Appl. Numer. Math. (2012).

Equazioni differenziali stocastiche

Problemi di Ito:

dX(t) = f(X(t))dt+ g(X(t)) dW (t), t ≥ 0,

f : Rd → Rd (drift), g : Rd → Rd×m (diffusione),

W (t) processo di Wiener m-dimensionale.Formulazione integrale:

X(t) = X(0) +

∫ t

0f(X(s))ds+

∫ t

0g(X(s))dW (s)︸ ︷︷ ︸

integrale di Ito∗

.

*Su una rete uniforme {0 < t1 < t2 < · · · ≤ tn}, e definito come

limn→∞

n∑j=1

g(X(tj))(W (tj+1)−W (tj)).

Incrementi di Wiener distribuiti come√h · N (0, 1), con h = tj+1 − tj .

Integrazione geometrica (stocastica)

Integrazione geometrica (stocastica)

Invarianza di leggi asintotiche nella discretizzazione di sistemilineari (Schurz, 1999);

conservazione a lungo termine della densita stazionaria dioscillatori lineari smorzati, con rumore additivo (Burrage, Lythe,

2007, 2009); D’Ambrosio, Moccaldi, Paternoster, 2018; Citro, D’Ambrosio, 2020);

metodi partizionati per oscillatori lineari con rumore additivo(Melbo, Higham, 2004);

oscillatore di Duffing stocastico (Welfert, 2017);

Runge-Kutta per Hamiltoniane stocastiche separabili (Burrage

2012, 2014).

integrazione geometrica di problemi Hamiltoniani stocastici.

C. Chen, D. Cohen, R. D’Ambrosio, A. Lang, Drift-preserving numericalintegrators for stochastic Hamiltonian systems, Adv. Comp. Math.(2020).

Sistemi Hamiltoniani stocastici

Per Hamiltoniane della forma

H(p(t), q(t)) =1

2pTp+ V (q), t ≥ 0,

con V : Rd → R potenziale sufficientemente regolare, consideriamo

dq(t) = p(t) dt

dp(t) = −V ′(q(t)) dt+ Σ dW (t).

con Σ ∈ Rd×m.

Legge invariante di traccia (Burrage, 2014)

E [H(p(t), q(t))] = E [H(p(t0), q(t0))] +1

2Tr(

Σ>Σ)t

Crescita lineare della funzione Hamiltonian attesa.

Questa proprieta si conserva lungo discretizzazioni del sistema?

Sistemi Hamiltoniani stocastici

Per Hamiltoniane della forma

H(p(t), q(t)) =1

2pTp+ V (q), t ≥ 0,

con V : Rd → R potenziale sufficientemente regolare, consideriamo

dq(t) = p(t) dt

dp(t) = −V ′(q(t)) dt+ Σ dW (t).

con Σ ∈ Rd×m.

Legge invariante di traccia (Burrage, 2014)

E [H(p(t), q(t))] = E [H(p(t0), q(t0))] +1

2Tr(

Σ>Σ)t

Crescita lineare della funzione Hamiltonian attesa.

Questa proprieta si conserva lungo discretizzazioni del sistema?

Sistemi Hamiltoniani stocastici

Per Hamiltoniane della forma

H(p(t), q(t)) =1

2pTp+ V (q), t ≥ 0,

con V : Rd → R potenziale sufficientemente regolare, consideriamo

dq(t) = p(t) dt

dp(t) = −V ′(q(t)) dt+ Σ dW (t).

con Σ ∈ Rd×m.

Legge invariante di traccia (Burrage, 2014)

E [H(p(t), q(t))] = E [H(p(t0), q(t0))] +1

2Tr(

Σ>Σ)t

Crescita lineare della funzione Hamiltonian attesa.

Questa proprieta si conserva lungo discretizzazioni del sistema?

Conservazione numerica della legge di traccia

In Burrage (2014) emerge (dalla sperimentazione) che

la perturbazione stocastica di metodi simplettici non conserva lalegge di traccia;

lo stesso accade per metodi energy-preserving.

V. Citro, R. D’Ambrosio, Long-term analysis of stochastic theta-methodsfor damped stochastic oscillators, Appl. Numer Math. (2020).

Esempio: singolo processo di Wiener

dq(t) = p(t) dt

dp(t) = −V ′(q(t)) dt+ σ dW (t).

La parte lineare delle σ-espansioni di p e q contiene il σ√t (termine

secolare).

Conservazione numerica della legge di traccia

C. Chen, D. Cohen, R. D’Ambrosio, A. Lang, Drift-preserving numerical integrators forstochastic Hamiltonian systems, Adv. Comp. Math. (2020).

Ψn+1 = pn + Σ∆Wn −h

2

∫ 1

0V ′(qn + shΨn+1)ds,

qn+1 = qn + hΨn+1, (?)

pn+1 = pn + Σ∆Wn − h∫ 1

0V ′(qn + shΨn+1) ds.

Teorema

Se V ∈ C1(Rd), lo schema (?) soddisfa la legge di traccia

E [H(pn, qn)] = E [H(p(t0), q(t0))] +1

2Tr(

Σ>Σ)tn,

per ogni tn = nh punto di rete.

Conservazione numerica della legge di traccia

Teorema

Supponiamo V ′ Lipschitz continuo. Allora, esiste h∗ > 0 tale che, perogni 0 < h ≤ h∗, lo schema numerico (?) converge con ordine 1 inmedia quadratica,(

E[|q(tn)− qn|2

])1/2+(E[|p(tn)− pn|2

])1/2 ≤ Ch,ove la costante C e indipendente da h e da n.

Altre proprieta

Analisi dell’errore in senso debole e forte;

analisi dell’errore in relazione allo stimatore del valore medio(Monte Carlo, Monte Carlo multilevel, ...);

analisi dell’errore Hamiltoniano.

Test numerici (pendolo stocastico)

H(p, q) =1

2p2 − cos(q), σ = 0.25, (p0, q0) = (1,

√2).

0 1 2 3 4 5

Time

0.3

0.35

0.4

0.45

0.5

0.55

0.6

Energy

EM

DP

BEM

Exact

0 5 10

Time

0.3

0.4

0.5

0.6

0.7

Energy

DP

Exact

Test numerici (doppia buca di potenziale)

H(p, q) =1

2p2 +

1

4q4 − 1

2q2, ε = 0.5, (p0, q0) = (

√2,√

2).

0 10 20 30 40 50

Time

0

2

4

6

8

Energy

DP

Exact

Test numerici (Henon-Heiles con doppio rumore)

H(p, q) =1

2

(p2

1 + p22

)+

1

2

(q2

1 + q22

)+ α

(q1q

22 −

1

3q3

1

),

Σ =

(0.2 00 0.2

), p0 =

(11

), q0 =

(√3

1

), α = 1/16.

0 10 20 30 40 50

Time

3

3.5

4

4.5

5

Energy

DP

Exact

Conclusions

Primi passi nell’integrazione geometrico-numerica di SDEs

Conservazione di leggi invarianti

Stabilita non lineare

Altri operatori (SVIEs, SFDEs, oscillatori stocastici)

E. Buckwar, R. D’Ambrosio, Exponential mean-square stability properties of stochasticlinear multistep methods, submitted.

C. Chen, D. Cohen, R. D’Ambrosio, A. Lang, Drift-preserving numerical integratorsfor stochastic Hamitlonian systems, Adv. Comput. Math. (2020).

V. Citro, R. D’Ambrosio, Long-term analysis of stochastic theta-methods for dampedstochastic oscillators, Appl. Numer Math. (2020).

D. Conte, R. D’Ambrosio, B. Paternoster, On the stability of theta-methods forstochastic Volterra integral equations, Discr. Cont. Dyn. Sys.-B (2018).

R. D’Ambrosio, S. Di Giovacchino, Mean-square contractivity of stochastic θ-methods,submitted.

R. D’Ambrosio, M. Moccaldi, B. Paternoster, Numerical preservation of long-termdynamics by stochastic two-step methods, Discr. Cont. Dyn. Sys.-B (2018).

Assegni di ricerca a L’Aquila

PRIN2017-MIUR“Structure preserving approximation of evolutionary problems”

Primo bando: giugno 2020

Tematiche: integrazione structure-preserving di problemi di evoluzionedeterministici e stocastici

Grazie per l’attenzione!

ε-espansioni

Idea classica in ambito deterministico (e.g., problemisingolarmente perturbati; vedi Hairer, Wanner).Esempio (didattico):

x2 + 2εx− 1 = 0, ε� 1.

Ansatz: x =

∞∑n=0

xnεn. Sostituendo nell’equazione e raccogliendo le

potenze di ε, si ottiene

ε0 : x20 − 1 = 0

ε1 : 2x0(x1 + 1) = 0

ε2 : x21 + 2x0x2 + 2x1 = 0

x ≈ xε = ±1− ε± ε2

2

ε = 10−1 ε = 10−3

|x+ − x+ε | 1.24e-5 1.25e-13

ε-espansioni

Consideriamo il problema test scalare[dqdp

]=

[0 −11 0

]y dt+

[σ0

]dW,

corrispondente al problema del secondo ordine

q = −q + σξ(t),

ove ξ(t) e un rumore Gaussiano. Generalmente, la letteratura consideracome problema test

q = −q + ε√t

con ε ∼ N (0, σ2). Ansatz:

q(t) = q0(t) + εq1(t) + ε2q2(t) + . . .

ε-expansions (ctd.)

Sostituendo nell’equazione e raccogliendo le potenze di ε, si ottiene

ε0 : q0 + q0 = 0

ε1 : q1 + q1 =√t

...

Dunque,

q(t) = cos(t) + sin(t)− ε(√

π

2C

(√2

π

√t

)cos(t) +

√π

2S

(√2

π

√t

)sin(t) +

√t

)+ . . .

p(t) = cos(t)− sin(t) + ε

√π

2

(C

(√2

π

√t

)sin(t)− S

(√2

π

√t

)cos(t)

)+ . . .

ove S e C sono le funzioni integrali di Fresnel

S(x) =

∫ x

0sin(π

2t2)dt,

C(x) =

∫ x

0cos(π

2t2)dt.

ε-espansioni (ctd.)

q(t) = cos(t) + sin(t)− ε(√

π

2C

(√2

π

√t

)cos(t) +

√π

2S

(√2

π

√t

)sin(t) +

√t

)+ . . .

p(t) = cos(t)− sin(t) + ε

√π

2

(C

(√2

π

√t

)sin(t)− S

(√2

π

√t

)cos(t)

)+ . . .

√t: termine secolare;

nell’integrazione a lungo termine, la parte secolare e dominante;

nell’integrazione a lungo termine, la funziona Hamiltonianamodificata (i.e., quella associata alla ε-espansione) differiscesignificativamente da quella originale;

la conservazione di leggi invarianti, dunque, non e unacaratteristica generale dei metodi stocastici.

Figure: Hamiltonian modificata attesa (rosso) VS Hamiltoniana attesa (blu)per ε = 0.01; 0.1; 0.2; 0.5.