attempt at quacrature rules

Author: smidl

docs/src/lecture_12/lecture.md

@@ -166,15 +166,33 @@ In multivariate setting, the same problem is typically solved with the aim to re
 
 One of the most popular approaches today is based on cubature rules approximating the Gaussian in radial-spherical coordinates.
 
-decomposition of the covariance matrix 
+## Cubature rules
+
+Consider Gaussian distribution with mean ``\mu`` and covariance matrix ``\Sigma`` that is positive definite with square root ``\sqrt\Sigma``, such that ``\sqrt\Sigma^T \sqrt\Sigma=\Sigma``. The quadrature pints are:
+```math
+X_q = \mu .+ \sqrt\Sigma Q
+```
+where ``Q=[q_1,\ldots q_{2d}]`` are constant vectors
+```math
+Q = \sqrt{d} [ I_d -I_d]
+```
+with associated weights
 ```math
+w_i = \frac{1}{2d}, i=1,\ldots,sd
+```
+where ``d`` is dimension of the vectors.
 
+The quadrature points are propogated through the non-linearity and the resulting Gaussian distribution is:
+```math
+\begin{align}
+x' & \sim N(\mu',\Sigma')\\
+\mu' & = \frac{1}{2d}\sum_{j=1}^{2d} x_i\\
+\Sigma &= \frac{1}{2d}\sum_{j=1}^{2d} (x_i-\mu')^T (X_i-\mu')
+\end{align}
 ```
 
-# Manipulating ODEs
 
-So far, we have considered first-order ODEs. Many ODEs are defined in higher order form, e.g. 
+## Containing uncertainty in vectors
 
+With vectorized uncertainty, it is now harder to assign which variable is now the first, second, etc.
 
-- check Chris Raucausacc latest blog
-  https://julialang.org/blog/2021/10/DEQ/

docs/src/lecture_12/ode.jl

@@ -98,5 +98,34 @@ plot!(MX[2,1:30:end],label="y",color=:red)
 savefig("LV_Measurements2.svg")
 
 
+
 # Plot receipe
-# plot(Vector{GaussNum})
+# plot(Vector{GaussNum})
+
+using LinearAlgebra
+function solve(f,x0::AbstractVector,√Σ0, θ,dt,N)
+  const n = length(x0)
+  const n2 = 2*length(x0)
+  const Qp = sqrt(n2)*[I(n) -I(n)]
+
+  X = hcat([zero(x0) for i=1:N]...)
+  S = hcat([zero(x0) for i=1:N]...)
+  X[:,1]=x0
+  √Σ = √Σ0
+  Σ = √Σ* √Σ'
+  S[:,1]= diag(Σ)
+  for t=1:N-1
+    Xp = x0 .+ √Σ * Qp
+    for i=1:n2 # all quadrature points
+      Xp.=X[:,t] + dt*f(Xp[:,i],θ)
+    end
+    mXp=mean(Xp,dims=2)
+    X[:,t+1]=mXp
+    Σ=(Xp.-mXp)*(Xp.-mXp)'/sqrt(n2)
+    S[:,t+1]=diag(Σ)
+    √Σ = chol(Σ)
+  end
+  X,S
+end
+
+## Extension to arbitrary