search

SciML / DifferentialEquations.jl

Multi-language suite for high-performance solvers of differential equations and scientific machine learning (SciML) components. Ordinary differential equations (ODEs), stochastic differential equations (SDEs), delay differential equations (DDEs), differential-algebraic equations (DAEs), and more in Julia.
https://docs.sciml.ai/DiffEqDocs/stable/
Other
2.85k stars 226 forks source link

Differential Matrix Equations and Parareal #225

Closed maximilianbehr closed 6 years ago

maximilianbehr commented 6 years ago

Hi, i am a Phd Student at the Max Planck Institute Magdeburg. I am working on numerical methods for large scale differential Riccati equations. I want to use DifferentialEquations.jl to solve a medium scale differential Riccati equation,

,

where are dense matrices of medium size, lets say .

My first question: Is it possible to use DifferentialEquations.jl to solve equations like that ?

I am familiar with Julia, due to multiple dispatching and abstraction i suspected that this is possible.

If you apply numerical method to the differential Riccati equation you have to solve several algebraic Lyapunov equations or Riccati equations. I have special solvers for that. How can i use my solvers together with your package.

My second question: Can i specifiy a solver for the algebraic systems which have to be solved during integration.

Even for medium size Differential Riccati Equations it is quite time consuming to solve such systems, due to the fact that an ODE of size have to be solved. I started together with a colleague to implement Parareal in Julia. We have a basic implementation but it is not ready for publishing yet.

My last question: Have you already worked on adding Parareal to the DifferentialEquations.jl package?

My goal would be that our Parareal code works smoothly together with DifferentialEquations.jl.

Thanks in advance.

ChrisRackauckas commented 6 years ago

Hey,

Sorry it took awhile to get to this.

Is it possible to use DifferentialEquations.jl to solve equations like that ?

Yes, you just define the ODE on a matrix. I recently showed something like that very similar in a blog post:

http://www.stochasticlifestyle.com/solving-systems-stochastic-pdes-using-gpus-julia/

Additionally, the first tutorial on ODEs discusses the the dependent variable's type is defined by the initial condition, so if you want to solve an ODE on matrices just pass a matrix:

http://docs.juliadiffeq.org/latest/tutorials/ode_example.html#Example-3:-Using-Other-Types-for-Systems-of-Equations-1

I suggest reading the full ODE tutorial since it's a good overview of standard usage.

How can i use my solvers together with your package.

See this page on specifying linear solvers

http://docs.juliadiffeq.org/latest/features/linear_nonlinear.html

The implicit solvers for which that can be used is discussed on the ODE solvers page:

http://docs.juliadiffeq.org/latest/solvers/ode_solve.html#Extra-Options-1

Can i specifiy a solver for the algebraic systems which have to be solved during integration.

Do you have a differential-algebraic equation instead of an ODE? You might want to specify it in the fully implicit form then and use a Krylov method from Sundials. The specification of a DAE is shown here:

http://docs.juliadiffeq.org/latest/tutorials/dae_example.html

Or you can use a singular mass matrix with stiff ODE solvers in order to include algebraic equations, in which case you can add your own linear solver for the whole system which then allows you to take control of this as well.

Even for medium size Differential Riccati Equations it is quite time consuming to solve such systems, due to the fact that an ODE of size have to be solved.

Try to BDF + Krylov methods from Sundials. 200^2 really isn't that bad and should only take like 10's of seconds. Make sure you write everything in a non-allocating way. You may want to get some help optimizing your ODEs by positing on the Julia Discourse

https://discourse.julialang.org/

I go through how to optimize away allocations in that blog post as well.

Have you already worked on adding Parareal to the DifferentialEquations.jl package?

Not yet. There's two ways it can be done. One way is to make an algorithm in OrdinaryDiffEq.jl which wraps the perform_step! calls of another algorithm, and takes control of the steps to piece them together. The advantage of this is that it can be done in a way that's still compatible with things like event handling, but this implementation would take some know-how.

The other way to implement it is to use the integrator interface.

http://docs.juliadiffeq.org/latest/basics/integrator.html

You can directly control how a method steps (this is exactly how the internals are working BTW). Using that, you can make multiple independent integrators and piece them together to a parareal algorithm. I think this would actually be a good first step to then create the one that's internal to OrdinaryDiffEq.jl.

The nice thing about this setup is that it then works over any timestepping method of OrdinaryDiffEq.jl, which instantly gives a whole lot of choices. But you can also define it the more traditional way of just defining the update equation. It's quite easy to add algorithms to OrdinaryDiffEq.jl since you just need to define the method that says "given a value at t, compute the value at t+dt" and the rest is done for you. The devdocs on this are here:

http://devdocs.juliadiffeq.org/latest/contributing/adding_algorithms.html

I am quite interested in getting some parallel algorithms up and running. Last summer as part of a Google Summer of Code (students get paid :) ) we built solvers which work by training a neural net in parallel, but by the end we were able to show it's just not a great method.

https://julialang.org/blog/2017/10/gsoc-NeuralNetDiffEq

Parareal is probably a better choice, so I'm ready for round two, especially if any student wants to take this on with us!

I'm closing this because there isn't anything actionable here, but feel free to keep the discussion going. Or, for more Julia usage questions, you may want to post to the Discourse or use our chat channel ( https://github.com/JuliaDiffEq/DifferentialEquations.jl ). Hopefully this helps!