search

Neuroblox / Spectral-Dynamic-Causal-Modeling

2 stars 2 forks source link

integration_approximation.jl produces rather poor approximations #1

Closed david-hofmann closed 2 years ago

david-hofmann commented 2 years ago

In integration_approximation.jl I have implemented the way SPM12 handles linearization of the right hand side of an ODE and numerical approximation of the solution of a non-linear ODE. Testing on a linear system gives exact results as expected. However, already weakly non-linear systems, like the generalized Lotka-Volterra equations, produces substantial deviations. (Same for the Lorenz system, both models are in the code, ready to be tested)

The numerical approximation is simply: x_approx[i, :] = expv(dt, J_x(x_approx[i-1, :]), x_approx[i-1, :]) Where J_x is the Jacobian w.r.t. the right hand side of the ODE system.

In the MATLAB implementation they provide another routine that suggests to approximatively integrate based on: dx = (expm(dt*J) - I)*inv(J)*f where f is (supposedly) the value of the right hand side of the ODE. I'll test this next unless someone has deeper insights and a different suggestion.

ChrisRackauckas commented 2 years ago

Show where in the code and what exactly you're trying to do here. That update looks like it's matching a first order exponential Rosenbrock method https://diffeq.sciml.ai/stable/solvers/split_ode_solve/#OrdinaryDiffEq.jl-2, specifically IIRC it might be LawsonEuler in Rosenbrock form. But that's only a good idea if f(u,t) = f(t) and even then I would suggest using a higher order Magnus method, so I'd like to know a bit more about why this integration is using this form.

david-hofmann commented 2 years ago

I'll show it tomorrow at the meeting, that's easier

ChrisRackauckas commented 2 years ago

I think this is the same issue as https://github.com/Neuroblox/Spectral-DCM/issues/2#issuecomment-952082271 .