Awkward ODE problem

Specific Domains Modelling & Simulations
,
1

I have a differential equation of the kind

\dot q_1 = f(q_1, q_2)\\ A(q_1,q_2)\dot q_2=B(q_1,q_2)\dot q_1

where q_1, q_2 are vectors, A is a symmetric, positive-definite real matrix, and B is a real (non-square) matrix.

When I use Tsit5 or RK4, the time step dt fluctuates a lot and sometimes the solver goes into a loop (evaluating the r.h.s. over and over without moving forward in time). The problem goes away when I set \dot q_2=0.

So, I suspect that the solution to the second equation (the linear equation) is not good enough and thus the error estimate for adaptive methods is unreliable.

Any ideas what could be a better approach here? Am I using the wrong solvers?

2

Is this setup using mass matrices? That wouldn’t work with explicit methods so I’m a bit confused there.

How do the implicit or Rosenbrock methods do?

3

In practice I am inverting A (it should be invertible). I am not using the mass matrix solvers in OrdinaryDiffEq (yet)

EDIT: I should be more precise. I run the problem as

\dot q_2 = A^{-1}B\dot q_1

using ldiv from LinearAlgebra

4

The vectors q_1 and q_2 have combined 4000 components. I’m not sure implicit / Rosenbrock methods are feasible here. I have no idea how I would calculate the Jacobian in a practical amount of time.

5

try FBDF. it’s often really good.

6

Thanks, I tried ImplicitEuler with Krylov GMRES and autodiff=false which works (albeit very slowly).

I’ll try FBDF next

7

ImplicitEuler is essentially always slow. That’s kind of by definition though. FBDF, QNDF, etc. try the efficient methods if you need efficiency.

2 Likes
8

So, just to close this thread out:

I changed the definition of the problem so that the integrand fluctuates less and therefore appears less stiff to the solvers. Now it runs nicely with RK4() (and probably Tsit5() but each round of testing takes a long time so I won’t let perfect be the enemy of good and leave it at that).