Some Magnus Methods

[ChrisRackauckas]

Hairer III

• [x] Example 7.4 page 123 MagnusGauss4

IMPROVED HIGH ORDER INTEGRATORS BASED ON THE MAGNUS EXPANSION https://link.springer.com/content/pdf/10.1023/A:1022311628317.pdf

• [ ] MagnusBlanes2
• [ ] MagnusBlanes4
• [ ] MagnusBlanes6
• [ ] MagnusBlanes8
• [x] MagnusGL6
• [x] MagnusNC6
• [x] MagnusGL8
• [x] MagnusNC8
• [ ] MagnusHybrid4
• [ ] MagnusHybrid6

(and adaptive!)

https://arxiv.org/pdf/0810.5488.pdf

• [ ] MagnusGL4

[Biswajitghosh98]

MagnusGauss4 implemented in https://github.com/SciML/OrdinaryDiffEq.jl/pull/1175

[Biswajitghosh98]

MagnusGL6 implemented in https://github.com/SciML/OrdinaryDiffEq.jl/pull/1183 And MagnusNC6 implemented in https://github.com/SciML/OrdinaryDiffEq.jl/pull/1184

[Biswajitghosh98]

MagnusNC8 implemented in https://github.com/SciML/OrdinaryDiffEq.jl/pull/1188 and Magnus GL8 implemented in https://github.com/SciML/OrdinaryDiffEq.jl/pull/1186

[Biswajitghosh98]

Formally, given an IVP defined by the linear differential equation dX/dt = A(t)X , and X(t₀) = I , where I is a matrix containing the inital values of the dependent variables at t₀, we can obtain a solution of the IVP that can be written in the form X(t,t₀) = e^Ω(t,t₀).

Now, Ω is obtained as an infinite series Ω(t,t₀) = Σ Ω_k(t,t₀) , wherein the k^th term, or Ω_k in the series is a multiple integral of nested commutators containing k operators A(t).

What each of the implemeted algoritm does is, an algorithm of order n, say MagnusGL6, where n=6 numerically approximates the integral terms with quadratures so as to obtain the required accuracy of order O(dtⁿ⁺¹) for a sufficiently smooth matrix A(t)

To use any of these, one can simply :

• Define an update function for the coefficient matrix A(t)

  function update_func(A,u,p,t)
  A[1,1] = cos(t)
  A[2,1] = sin(t)
  A[1,2] = -sin(t)
  A[2,2] = cos(t)
  end

• Formulate the ODE Problem

  A = DiffEqArrayOperator(ones(2,2),update_func=update_func)
  prob = ODEProblem(A, ones(2), (1.0, 6.0))

• Use any of the listed algorithms to solve the problem, depending on the convergence order needed, which in our case was MagnusGL6()

  sol = solve(prob,MagnusGL6(),dt=1/10)

The generated solution object sol can now be handled separately.

NOTE : The MagnusBlanes and MagnusHybrid methods involved Taylor series expansion of of the Matrix A(t) about t = t₀+dt/2, and hence requiring the need to calculate k^th order derivatives. Hence these have been skipped for now