What is the state of solvers of differential equations?

1,287 views
Skip to first unread message

jonas.a...@gmail.com

unread,
Jul 21, 2013, 9:37:48 PM7/21/13
to julia...@googlegroups.com
To gently introduce myself to using Julia, I'd like to solve differential equations of simple circuit models. What is the current state of affairs on this topic? Are there any such solvers at hand?

The closest I can find is an ODE package, seemingly not under heavy development, and a discussion from more than a year ago that indicates that, at least then, there existed no such facilities.

Tom Short

unread,
Jul 21, 2013, 10:07:12 PM7/21/13
to julia...@googlegroups.com
Try the Sundials.jl package:

https://github.com/tshort/Sundials.jl

For circuit models, you might also be interested in Sims.jl, a package that can model circuits somewhat like Modelica:

https://github.com/tshort/Sims.jl

It uses either Sundials or DASKR, a DASSL variant. No guarantees Sims will work right now. It needs updating and probably won't work with the latest Julia.



Tomas Lycken

unread,
Jul 21, 2013, 10:21:13 PM7/21/13
to julia...@googlegroups.com
The Sundials package is an excellent wrapper around the C library with the same name, and seems to be excellent (although I admit I haven't used it myself). If you need an ODE solver now, I definitely think this is the place to look.

I've been interested in doing work on the ODE package for quite some time, but as with so many things I'm lacking the time to spend on it. Even though Sundials exists, I do think there is merit in having a pure-Julia package for at least basic to semi-advanced ODE solvers. If you have suggestions on things that should be added or improved, or resources that could be useful (for example algorithm descriptions with licenses that are compatible with the MIT licence), do let me know. Knowing that any work I put in actually helps someone other than just myself is a great motivator =)

// Tomas

Steven G. Johnson

unread,
Jul 22, 2013, 12:52:55 AM7/22/13
to julia...@googlegroups.com

On Sunday, July 21, 2013 6:21:13 PM UTC-4, Tomas Lycken wrote:
I've been interested in doing work on the ODE package for quite some time, but as with so many things I'm lacking the time to spend on it. Even though Sundials exists, I do think there is merit in having a pure-Julia package for at least basic to semi-advanced ODE solvers. If you have suggestions on things that should be added or improved, or resources that could be useful (for example algorithm descriptions with licenses that are compatible with the MIT licence), do let me know. Knowing that any work I put in actually helps someone other than just myself is a great motivator =)

Yes, a pure Julia routine in Base would be nice, at least a solid implementation of ode45 (adaptive Runge-Kutta) or similar.

For example, this would allow one to integrate ODEs y' = f(t,y) where f returns any normed vector space (any type supporting norm, +, -, and multiplication by real scalars) (analogous to the quadgk integration routine we have now).

We can leave it to external packages for fancier features (parallelization, stiff equations, differential-algebraic equations, etcetera).

See also https://github.com/JuliaLang/julia/issues/75

--SGJ

Tomas Lycken

unread,
Jul 22, 2013, 8:16:50 AM7/22/13
to julia...@googlegroups.com
There are a couple of different ode45 versions as well as an ode23 and some other solvers included in ODE.jl, so I don't know if I agree that this should go into Base (having "using ODE" at the top of the file is simple enough...). There API is not set in stone, and there is some discussion in https://github.com/vtjnash/ODE.jl/issues/4 where you can chime in.

Also, I should note that I'm not the maintainer of this package. However, I seem to be one of quite few developers who want to work on it at the moment :P

// Tomas

Simon Frost

unread,
Jul 22, 2013, 10:23:55 PM7/22/13
to julia...@googlegroups.com
Dear Jonas,

I've played a bit with both Sims.jl and ODE.jl; the former has a friendlier syntax. For example, the 'standard' SIR epidemiological model in Sims looks like the following.

require("Sims")
function SIR(beta,gamma)
    S=Unknown(9999,"S")
    I=Unknown(1,"I")
    R=Unknown(0,"R")
    N=S+I+R
   {
    der(S)-(-beta*S*I)
    der(I)-(beta*S*I-gamma*I)
    der(R)-(gamma*I)
   }
end
(S,I,R)=sim(SIR(0.1/10000,0.05),1000)

Alas, the above doesn't currently work in the current build of Julia.

I've blogged about the ODE.jl interface; if I'm not mistaken, there isn't a way to pass parameters to the solver, so parameters have to be included in the state vector (http://phylodynamics.blogspot.co.uk/2013/07/differential-equation-modeling-with.html). It would be nice to have both a lower level (cf deSolve in R) and a higher level interface (cf simecol in R) to ODE solvers.

Best
Simon

Viral Shah

unread,
Jul 23, 2013, 6:10:16 AM7/23/13
to julia...@googlegroups.com
It would be great if you can file issues against the Sims package - packages often fall behind as the base moves fast. Hopefully, in a couple more releases, we will no longer need everyone to be using master.

-viral

Steven G. Johnson

unread,
Jul 23, 2013, 3:53:49 PM7/23/13
to julia...@googlegroups.com
On Monday, July 22, 2013 6:23:55 PM UTC-4, Simon Frost wrote:
I've blogged about the ODE.jl interface; if I'm not mistaken, there isn't a way to pass parameters to the solver, so parameters have to be included in the state vector (http://phylodynamics.blogspot.co.uk/2013/07/differential-equation-modeling-with.html). It would be nice to have both a lower level (cf deSolve in R) and a higher level interface (cf simecol in R) to ODE solvers.

Lexical scoping and anonymous functions are your friends here.  If you are solving x' = f(t,x, p), with p being some additional parameter, and you have the function f, you can just do:
         ode23((t,x) -> f(t,x,p), trange, x0)
in order to pass some additional parameteter p to f.

I'm less concerned at this point about small interface details than having a production-quality Julia implementation which supports arbitrary types (arbitrary Banach spaces).  It doesn't seem like we have that yet.  (The ODE.jl README explicitly says, "Some of these solvers are known to be produce incorrect results"!)

--SGJ

Tomas Lycken

unread,
Jul 23, 2013, 4:13:08 PM7/23/13
to julia...@googlegroups.com

I'm less concerned at this point about small interface details than having a production-quality Julia implementation which supports arbitrary types (arbitrary Banach spaces).  It doesn't seem like we have that yet.  (The ODE.jl README explicitly says, "Some of these solvers are known to be produce incorrect results"!)

If you look through the issues (there are only a handful of them) there is a problem with one of the solvers. The others should work fine - but one of the great lackings of the ODE package is that it doesn't have a good test suite. I've tried to do some research to find good test problems (non-trivial, yet with known solutions, and classified so that solvers can be tested with appropriate problems) but I'm stumbling and not getting very far. Ideally, I'd like to have at least one or two semi-difficult problems that are non-stiff, a couple of stiff problems and a couple of implicit ones. Any pointers to good resources are more than welcome.

As far as API decisions are concerned, there are other parameters as well that are nice to have - there's been discussion about how to specify tolerance arguments as well as the possibility to take a single step (mostly relevant for solvers with variable step lenght), but there might be other things that one wants to specify as well.

But the entire package is still, obviously, very much a work-in-(slow)-progress, so don't judge its potential after it's current state ;)

// Tomas

Tomas Lycken

unread,
Jul 23, 2013, 8:34:00 PM7/23/13
to julia...@googlegroups.com
Another question: what do you mean by supporting arbitrary Banach spaces? What would be required of the solver(s) to meet this requirement? Is it enough to be able to specify a function for the vector norm used to compare solution points, or does the solver need to do something else as well?

// T

Simon Frost

unread,
Jul 24, 2013, 1:06:44 PM7/24/13
to julia...@googlegroups.com
Dear Steven,

On Tuesday, July 23, 2013 4:53:49 PM UTC+1, Steven G. Johnson wrote:
Lexical scoping and anonymous functions are your friends here.  If you are solving x' = f(t,x, p), with p being some additional parameter, and you have the function f, you can just do:
         ode23((t,x) -> f(t,x,p), trange, x0)
in order to pass some additional parameteter p to f.

Thanks for the tip re: anonymous functions and lexical scoping; I'm still very much a Julia newbie. I'll try and hack a higher level interface similar to simecol, but as you point out, that doesn't fix the lower-level problems. Most, but not all, of the ODE systems I look at are rather well behaved. I'll have a go at trying to break the ODE solvers with some stiffer models.

Best
Simon
 

jonas.a...@gmail.com

unread,
Jul 26, 2013, 5:10:42 PM7/26/13
to julia...@googlegroups.com
Thanks for all the great replies! I didn't see them until now because I mistakenly thought I would be getting email notifications as they came in.

I now have some links and topics to check out.

Stefan Karpinski

unread,
Jul 26, 2013, 5:51:40 PM7/26/13
to julia...@googlegroups.com
You must not have your group settings setup to deliver emails. This can be adjusted.

jonas.a...@gmail.com

unread,
Jul 26, 2013, 6:17:57 PM7/26/13
to julia...@googlegroups.com
On Friday, July 26, 2013 7:51:40 PM UTC+2, Stefan Karpinski wrote:
You must not have your group settings setup to deliver emails. This can be adjusted.

 It actually works now all of a sudden. I must have switched it on unknowingly. It makes you wonder what else I do unknowingly!

Tomas Lycken

unread,
Aug 5, 2013, 10:30:51 AM8/5/13
to julia...@googlegroups.com
I'd like to revive this discussion - I submitted a pull request for the ODE.jl package a while ago which has received very little attention, so if more people than myself are interested in the state of the ODE.jl package it would be wonderful with some comments there.

The pull request doesn't nearly do all the things I'd like to do, but getting the API right is difficult, and just because I think something is a good idea doesn't mean it necessarily is...

// Tomas
Reply all
Reply to author
Forward
0 new messages