ImplicitMidpoint and other implicit methods: option for fixed point iteration as nonlinear solver?

[AnderMuruaUria]

I'd like to solve a system of ODEs with the algorithm ImplicitMidpoint() implemented with fixed point iteration. I've looked into the documentation of DifferentialEquations.jl, and have found no simple way to do it. Should I explicitly define a nonlinear solver ImplicitMidpoint(nlsolve=?) that does the trick?

[ChrisRackauckas]

Should I explicitly define a nonlinear solver ImplicitMidpoint(nlsolve=?) that does the trick?

Yes. If you have a problem setup as G(x)=0, then f(x) = G(x) + x is the fixed point problem, so applying a fixed point solver like Picard or Anderson acceleration to that will do it. So the way it's setup now you'd want to pass an nlsolve problem that does that. However, that is admittedly much too cumbersome and we need a better way to handle this.

Normally functional iteration is bad for stability and thus cannot be used for stiff equations. I assume this is because you want to use ImplicitMidpoint for the symplectic properties on a non-stiff equation? If so, that's a good use and we should make the option much easier for this (though of course, if the system is partitioned, make use of the symplectic RK methods).

[AnderMuruaUria]

I want to use ImplicitMidpoint because it is time-reversible and it exactly conserves quadratic invariants. I'm not interested in high order methods, but rather in a low order method with constant time-step that tries to capture the symmetries of the problem.

My problem has no partitioned structure, so that PRK symplectic methods cannot be used. Due to the particularities of the problem at hand, I have to use very small time-steps, which makes possible (and very convenient, because its Jacobian is very expensive dense matrix) to use a fixed point iteration.

Of course, I could implement the implicit midpoint rule with constant time-step myself, but I am very much interested in learning more about the powerful and versatile functionalities of DifferentialEquations.jl. Unfortunately, I'm not currently familiar enough with the way the algorithms in DifferentialEquations.jl deal with nonlinear solvers.

In principle, a simple way to effectively apply fixed point iteration in a framework of Newton-like iterations is to "approximate" the Jacobian of the right-hand side of the ODE by a zero matix. In other words, to replace the basic linear solve v -> (I - h c J)^{-1}v common in several implementation of implicit methods by just v -> v. Is there a simple way to do that for the implicit integrators in DifferentialEquations (in particular, for the algorithm ImplicitMidpoint)?

[ChrisRackauckas]

One way to "hack" it in there is to define a zero Jacobian:

http://docs.juliadiffeq.org/latest/features/performance_overloads.html#Declaring-Explicit-Jacobians-1

But I agree this is an important use-case so we should find a way to support this in the non-generic nlsolve algorithms. It wasn't in the current implementations but it would make sense to add Picard and Anderson accelerated iteration.

[AnderMuruaUria]

Thanks for your quick answer. Yes, defining a zero Jacobian will do the trick, but at the cost of some loss of efficiency, since it will be unnecessarily solving linear systems with identity coefficient matrix. Is that right?

On the other hand, I have implemented my problem as a system u'=f(u,p,t) where u is a Nx3 matrix. How should I define explicitly the zero Jacobian in that case?

[ChrisRackauckas]

On the other hand, I have implemented my problem as a system u'=f(u,p,t) where u is a Nx3 matrix. How should I define explicitly the zero Jacobian in that case?

If you're looking for efficiency you should be using the in-place form. This is described in the optimizing notebook:

http://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqTutorials.jl/blob/master/Introduction/OptimizingDiffEqCode.ipynb

But this question, along with

but at the cost of some loss of efficiency, since it will be unnecessarily solving linear systems with identity coefficient matrix. Is that right?

are being addressed in https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/220 . Jacobians for out of place methods and allowing user-chosen lazy Jacobian types (so that I-gamma*J = I, the UniformScaling operator instead of an array, which makes factorization free. That would be possible with that change). So there's three things involved here. I plan to do the Jacobian types within the next month, but iteration and out of place Jacobian usage will take a little bit more.

[AnderMuruaUria]

Thanks a lot. (By the way, I was already using the in-place form.)

[AnderMuruaUria]

I'm still trying to apply ImplicitMidpoint with fixed point iteration without success. I've tried to define a custom non-linear solver for the option nlsolve, but apparently the option nlsolve is not available for ImplicitMidpoint nor for Trapezoid. However, I've checked in the stable documentation and in the latest documentation (for the version of the documentation which contains the un-released features), and according to both versions of the documentation, the option nlsolve seems to be available, at least for Trapezoid. Am I wrong?

[ChrisRackauckas]

Yes sorry, that's what I meant by:

But I agree this is an important use-case so we should find a way to support this in the non-generic nlsolve algorithms.

Most of the algorithms don't allow generic nlsolve and so it's not supported yet. It'll take some work. Sorry that wasn't more clear.

[AnderMuruaUria]

Thanks a lot. You've done (and are doing) a great job with DifferentialEquations.jl!

I now realize that you already said that non-generic nlsolve algorithms aren't available. What has confused me is the fact that in (both in the stable and the newest version of) the documentation one can explicitly read that non-generic nlsolve are available, at least for the trapezoidal rule.

[ChrisRackauckas]

What has confused me is the fact that in (both in the stable and the newest version of) the documentation one can explicitly read that non-generic nlsolve are available, at least for the trapezoidal rule.

Yes, for specific algorithms its exists like GenericImplicitEuler and GenericTrapezoid. But those cases are rare.

[AnderMuruaUria]

I've finally succeeded in applying ImplicitMidpoint with fixed point iteration in an efficient way: prob = ODEProblem(f,u0,tspan,p); sol = solve(prob,ImplicitMidpoint(linsolve=linsolve!) with linsolve!(::Type{Val{:init}},f,x) = function _linsolve!(x,A,b,update_matrix=false) x .= b end In order to avoid redundant evaluations of the Jacobian, I set f(::Type{Val{:jac}},J,u,p,t) = nothing

[ChrisRackauckas]

Yeah, we have to make that smoother...

[ChrisRackauckas]

This PR is solving this issue: https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/pull/451

[ChrisRackauckas]

Added. Now the implicit methods can do nlsolve = NLFunctional(). Will get docs soon (@YingboMa )