Integrals

[brownadder]

The aim is to use Magnus expansion up to 4th order for non-stochastic. We need mu0 and mu1 line integrals and Lambda(f,g), Lambda(g,f) integrals.

mu0(f, t, h, method='gl7') $\approx \int_0^h f(t+s) ds =: \mu_0^f$ with 7 Gauss-Legendre points
mu1(f, t, h, method='scipy') $\approx \int_0^h (s-h/2) f(t+s) ds =: \mu_1^f$ with scipy (we need to prescribe a specific accuracy to scipy, not too high)

nu(f, g, t, h, method='gl3') $\approx \int_0^h \int_0^\zeta f(t+\zeta) g(t+\xi) d \xi d \zeta =: \nu^{f,g}$. In many cases, we need the other way around too, e.g. cross product formulae, in which case we can compute nu(g, f, t, h, method='gl3').

Guannan can you please correct if the above order is incorrect?

Doing so separately in different functions keeps us flexible, makes it easy to do unit testing etc. However, we can also make things more efficient by exposing the Gauss--Legendre knots explicitly, e.g. by creating a function quadweights which returns the knots and weights for GL quadrature (for example)

knots, weights = quadweights('gl7')

Then we sample f and g at these knots

fv = [f(t+k*h) for k in knots]
gv = [g(t+k*h) for k in knots]

Then passing fv instead of f,

mu0(fv, t, h, method='gl7')

f should no longer need to be sampled inside the method. The method should detect whether we are passing a function handle or a list of numbers.

We have discarded weights, which will be computed or used internally again (ideally not re-computed). Similarly for nu method. This way we sample each f at 7 GL knots and then pass them on to these functions.

[brownadder]

Proposed change in signature:

mu0(fv, timegrid, n, method='gl7')

Here we are specifying the full timegrid (as a list) and asking for integration over $[t_n, t_{n+1}]$. This might make syntax easier later.

For convergence test timegrid = [0,h], n=0 should do.