Galerkin-based Solution of Random Differential Equation

This tutorial demonstrates how random differential equations can be solved using polynomial chaos expansions (PCE).

Theory

A random differential equation is an ordinary differential equation that has random parameters, hence its solution is itself a (time-varying) random variable. Perhaps the simplest non-trivial example is the following scalar, linear ordinary differential equation

\[\dot{x}(t) = a x(t), \quad x(0) = x_{0},\]

where $a$ is the realization of a Gaussian random variable $\mathsf{a} \sim \mathcal{N}(\mu, \sigma^2)$ with mean $\mu$ and variance $\sigma^2$. Arguably, for every realization $a$ we can solve the differential equation and obtain

\[x(t) = x_0 \mathrm{e}^{a t},\]

from which we find that

\[\ln (x(t)) = \ln (x_0) + at \sim \mathcal{N}(\ln(x_0) + \mu t, (\sigma t)^2).\]

In other words, the logarithm of the solution is normally distributed (so-called log-normal distribution).

We'd like to obtain this result numerically with the help of PCE. The first step is to define the (truncated) PCE for the random variable $\mathsf{a}$

\[\mathsf{a} = \sum_{i=0}^{L} a_i \phi_i,\]

where $a_i$ are the so-called PCE coefficients, and $\phi_i$ are the orthogonal basis polynomials. As the solution to the random differential equation is itself a random variable, we treat $x(t)$ as the realization of the random variable $\mathsf{x}(t)$, and define its PCE

\[\mathsf{x}(t) = \sum_{i=0}^{L} x_i(t) \phi_i.\]

The question is how to obtain the unknown PCE coefficients $x_i(t)$ from the known PCE coefficients $a_i$ relative to the orthogonal basis polynomials $\phi_i$. This can be done using Galerkin projection, which is nothing else than projecting onto the orthogonal basis. Think of a three-dimensional space, in which you have placed some three-dimensional object. If you know project the silhouett of the object onto every axis of the three-dimensional space, then you are doing a Galerkin projection. With PCE the concept is equivalent, but the imagination has a harder time. The first step for Galerkin projection is to insert the PCEs

\[\sum_{i=0}^{L} \dot{x}_i(t) \phi_i = \sum_{j=0}^{L} a_j \phi_j \sum_{k=0}^{L} x_k(t) \phi_k;\]

the second step is to project onto every basis polynomial $\phi_m$ for $m = 0, 1, \dots, L$, and to exploit orthogonality of the basis. This gives

\[\dot{x}_m(t) \langle \phi_m, \phi_m \rangle = \sum_{j=0}^{L} \sum_{k=0}^{L} a_j x_k(t) \langle \phi_l \phi_k, \phi_m \rangle \quad m = 0, 1, \dots, L.\]

Of course, the initial condition must not be forgotten:

\[x_0(0) = x_0, \quad x_m(0) = 0 \quad m = 1, \dots, L.\]

If we can solve this enlarged system of ordinary random differential equations, we can reconstruct the analytic solution.

Practice

We begin by defining the random differential equation

x0 = 2.0
μ, σ = -0.5, 0.05
tend, Δt = 3.0, 0.01
(3.0, 0.01)

Next, we define an orthogonal basis (and its quadrature rule) relative to the Gaussian measure using PolyChaos. We choose a maximum degree of L.

using PolyChaos
L, Nrec = 6, 40
opq = GaussOrthoPoly(L; Nrec=Nrec, addQuadrature=true)
GaussOrthoPoly{Vector{Float64}, GaussMeasure, Quad{Float64, Vector{Float64}}}(6, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [1.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0  …  30.0, 31.0, 32.0, 33.0, 34.0, 35.0, 36.0, 37.0, 38.0, 39.0], GaussMeasure(PolyChaos.w_gaussian, (-Inf, Inf), true), Quad{Float64, Vector{Float64}}("golubwelsch", 39, [-11.289716044476338, -10.31334144275699, -9.501306853782872, -8.773438698067654, -8.099152213152408, -7.46268237786604, -6.854552519790957, -6.268491549685282, -5.700062454271861, -5.1459645440941815  …  5.145964544094184, 5.700062454271859, 6.268491549685282, 6.854552519790963, 7.46268237786604, 8.099152213152397, 8.77343869806765, 9.501306853782866, 10.313341442756998, 11.289716044476336], [9.443344575063092e-29, 2.7847875052254124e-24, 7.592450542206571e-21, 5.370701458462944e-18, 1.4861295877333035e-15, 1.998726680623688e-13, 1.491720010448911e-11, 6.748646364787844e-10, 1.9698729205993654e-8, 3.8841182837133097e-7  …  3.884118283713288e-7, 1.9698729205993754e-8, 6.748646364787816e-10, 1.4917200104488884e-11, 1.9987266806237089e-13, 1.4861295877332983e-15, 5.370701458462905e-18, 7.592450542206324e-21, 2.784787505225327e-24, 9.443344575062904e-29]))

Now we can define the PCE for $\mathsf{a}$ and solve the Galerkin-projected ordinary differential equation using DifferentialEquations.jl.

using DifferentialEquations

a = [ convert2affinePCE(μ, σ, opq); zeros(Float64,L-1) ] # PCE coefficients of a
xinit = [ x0; zeros(Float64,L) ] # PCE coefficients of initial condition

t2 = Tensor(2, opq); # \langle \phi_i, \phi_j \rangle
t3 = Tensor(3, opq); # \langle \phi_i \phi_j, \phi_k \rangle

# Galerkin-projected random differential equation
function ODEgalerkin(du,u,p,t)
   du[:] = [ sum( p[j+1]*u[k+1]*t3.get([j,k,m])/t2.get([m,m]) for j=0:L for k=0:L) for m=0:L ]
end

probgalerkin = ODEProblem(ODEgalerkin,xinit,(0,tend),a)
solgalerkin = solve(probgalerkin;saveat=0:Δt:tend)
t, x = solgalerkin.t, solgalerkin.u;
([0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09  …  2.91, 2.92, 2.93, 2.94, 2.95, 2.96, 2.97, 2.98, 2.99, 3.0], [[2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [1.9900252071385, 0.0009950126035692532, 2.4875315089174326e-7, 4.145885853203447e-11, 5.182355424144811e-15, 5.182372859052201e-19, 4.446367982010874e-23], [1.9801006575270281, 0.0019801006679350595, 9.90048283875986e-7, 3.3022172308019874e-10, 7.191587836603058e-14, 2.552297839257698e-16, -9.33641712389512e-19], [1.9702260956550872, 0.00295533916993449, 2.2164991440224486e-6, 1.1087782320980687e-9, 3.880756852155503e-13, 7.647078374450021e-16, -3.0601683025036038e-18], [1.9604012673415023, 0.003920802569373895, 3.920795659105214e-6, 2.6145693296523704e-9, 1.2695719817981672e-12, 1.4309116759425873e-15, -5.331499977669577e-18], [1.9506259196896107, 0.004876564830716209, 6.095699691025283e-6, 5.080409655788518e-9, 3.1388591538673724e-12, 2.5226473781782474e-15, -6.64782295583295e-18], [1.94089980108726, 0.005822699423081009, 8.73404503581089e-6, 8.734487546531226e-9, 6.524935767106835e-12, 4.675567823573074e-15, -5.8575450018776286e-18], [1.9312226612068097, 0.006759279320244441, 1.1828737423617146e-5, 1.3800365033883849e-8, 1.2063344567128056e-11, 8.892172692208805e-15, -1.7572954276259113e-18], [1.92159425100513, 0.007686377000639288, 1.537275451872634e-5, 2.0496977845828917e-8, 2.0496172479472406e-11, 1.6541808464828514e-14, 6.908074908161572e-18], [1.9120143227236033, 0.008604064447354927, 1.9359145919547708e-5, 2.9038635406327716e-8, 3.267205060960937e-11, 2.9360668422835965e-14, 2.1445493599785707e-17]  …  [0.4717682324476054, 0.06864220151579947, 0.004993736418123432, 0.00024219450240019178, 8.80990541705698e-6, 2.563666881018032e-7, 6.2003759829759125e-9], [0.4694494907228333, 0.0685395483131638, 0.0050034034535760006, 0.00024349724514500056, 8.887729361711572e-6, 2.5952015845416374e-7, 6.298100919935763e-9], [0.4671422619436801, 0.06843626325909682, 0.005012972845694051, 0.00024479844355019794, 8.965822463898332e-6, 2.626970759581227e-7, 6.396916976417835e-9], [0.46484648838474385, 0.06833235513846005, 0.0050224447642581525, 0.00024609805746775206, 9.044182650674616e-6, 2.658974600039844e-7, 6.496830428797749e-9], [0.46256211263138935, 0.06822783265777074, 0.005031819383033409, 0.0002473960470441747, 9.122807836618814e-6, 2.691213288279263e-7, 6.597847538468639e-9], [0.4602890775797432, 0.06812270444520165, 0.005041096879769488, 0.00024869237272052427, 9.20169592383053e-6, 2.7236869951200615e-7, 6.699974551841412e-9], [0.45802732643669736, 0.06801697905058114, 0.0050502774362006005, 0.0002499869952324037, 9.280844801930453e-6, 2.7563958798415703e-7, 6.803217700344562e-9], [0.45577680271990884, 0.06791066494539329, 0.005059361238045507, 0.00025127987560996076, 9.360252348060376e-6, 2.789340090181873e-7, 6.907583200424187e-9], [0.4535374502577972, 0.06780377052277767, 0.0050683484750075215, 0.000252570975177889, 9.439916426883263e-6, 2.822519762337845e-7, 7.013077253544102e-9], [0.45130921318954764, 0.06769630409752952, 0.005077239340774506, 0.00025386025555542666, 9.519834890583166e-6, 2.8559350209651094e-7, 7.119706046185711e-9]])

For later purposes we compute the expected value and the standard deviation at all time instants using PCE.

# an advantage of PCE is that moments can be computed from the PCE coefficients alone; no sampling required
mean_pce = [ mean(x_, opq) for x_ in x]
std_pce = [ std(x_, opq) for x_ in x]
301-element Vector{Float64}:
 0.0
 0.0009950126657575441
 0.0019801011629582794
 0.0029553408323057034
 0.003920806490165915
 0.004876572450342814
 0.005822712524166983
 0.006759300020585337
 0.007686407746251845
 0.00860410800561906
 ⋮
 0.06890642513273801
 0.06880511185670363
 0.06870317328013556
 0.06860061804929832
 0.068497454732652
 0.06839369182086613
 0.06828933772683328
 0.06818440078568275
 0.06807888925479469

We compare the solution from PCE to a Monte-Carlo-based solution. That means to solve the ordinary differential equation for many samples of $\mathsf{a}$. We first sample from the measure using sampleMeasure, and then generate samples of $\mathsf{a}$ using evaluatePCE. After that we solve the ODE and store the results in xmc.

using Statistics
Nsmpl = 5000
ξ = sampleMeasure(Nsmpl, opq)     # sample from Gaussian measure; effectively randn() here
asmpl = evaluatePCE(a, ξ, opq)     # sample random variable with PCE coefficients a; effectively μ + σ*randn() here
# or: asmpl = samplePCE(Nsmpl,a,opq)
xmc = [ solve(ODEProblem((u,p,t)->aa*u,x0,(0,tend));saveat=0:Δt:tend).u for aa in asmpl]
xmc = hcat(xmc...);
301×5000 Matrix{Float64}:
 2.0       2.0       2.0       2.0       …  2.0       2.0       2.0
 1.9899    1.98973   1.98964   1.98811      1.99041   1.98936   1.98875
 1.97984   1.97952   1.97933   1.97629      1.98086   1.97878   1.97756
 1.96984   1.96936   1.96908   1.96455      1.97136   1.96825   1.96644
 1.95989   1.95925   1.95888   1.95287      1.9619    1.95778   1.95538
 1.94998   1.94919   1.94873   1.94126   …  1.95249   1.94737   1.94438
 1.94013   1.93919   1.93864   1.92972      1.94312   1.93701   1.93344
 1.93033   1.92923   1.92859   1.91825      1.9338    1.92671   1.92256
 1.92058   1.91933   1.9186    1.90685      1.92452   1.91646   1.91175
 1.91087   1.90948   1.90866   1.89552      1.91529   1.90627   1.90099
 ⋮                                       ⋱                      
 0.455701  0.445038  0.438924  0.350743     0.491177  0.421386  0.385173
 0.453399  0.442753  0.43665   0.348658     0.488821  0.419145  0.383007
 0.451108  0.44048   0.434388  0.346586     0.486476  0.416915  0.380852
 0.448829  0.438219  0.432138  0.344526  …  0.484143  0.414697  0.37871
 0.446561  0.43597   0.429899  0.342478     0.48182   0.412492  0.37658
 0.444305  0.433732  0.427672  0.340442     0.479509  0.410297  0.374461
 0.44206   0.431505  0.425457  0.338418     0.477209  0.408115  0.372355
 0.439827  0.429291  0.423253  0.336407     0.47492   0.405944  0.37026
 0.437604  0.427087  0.42106   0.334407  …  0.472641  0.403785  0.368177

Now we can compare the Monte Carlo mean and standard deviation to the expression from PCE for every time instant.

[ mean(xmc,dims=2)-mean_pce std(xmc,dims=2)-std_pce]
301×2 Matrix{Float64}:
 0.0           0.0
 1.83207e-5   -3.11003e-5
 3.64283e-5   -6.18625e-5
 5.43246e-5   -9.22895e-5
 7.20115e-5   -0.000122384
 8.94907e-5   -0.000152149
 0.000106764  -0.000181588
 0.000123833  -0.000210702
 0.0001407    -0.000239496
 0.000157367  -0.000267971
 ⋮            
 0.00095933   -0.00199436
 0.000956853  -0.00199131
 0.00095437   -0.00198824
 0.000951881  -0.00198516
 0.000949385  -0.00198206
 0.000946885  -0.00197895
 0.000944378  -0.00197583
 0.000941866  -0.00197269
 0.00093935   -0.00196954

Clearly, the accuracy of PCE deteriorates over time. Possible remedies are to increase the dimension of PCE, and to tweak the tolerances of the integrator.

Finally, we compare whether the samples follow a log-normal distribution, and compare the result to the analytic mean and standard deviation.

logx_pce = [ log.(evaluatePCE(x_,ξ,opq)) for x_ in x]
[ mean.(logx_pce)-(log(x0) .+ μ*t) std.(logx_pce)-σ*t ]
301×2 Matrix{Float64}:
 -6.66134e-16   6.662e-16
  9.21397e-6   -1.56352e-5
  1.84279e-5   -3.12704e-5
  2.76419e-5   -4.69056e-5
  3.68559e-5   -6.25408e-5
  4.60698e-5   -7.8176e-5
  5.52838e-5   -9.38112e-5
  6.44978e-5   -0.000109446
  7.37118e-5   -0.000125082
  8.29258e-5   -0.000140717
  ⋮            
  0.00269075   -0.00456565
  0.00269997   -0.00458128
  0.00270919   -0.00459692
  0.0027184    -0.00461256
  0.00272762   -0.00462819
  0.00273683   -0.00464383
  0.00274605   -0.00465947
  0.00275526   -0.0046751
  0.00276447   -0.00469074