Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Already on GitHub? Sign in to your account

Stiff ODE solver + small fix #615

Merged
merged 2 commits into from Mar 20, 2012

Conversation

Projects
None yet
3 participants
Member

pao commented Mar 20, 2012

Finally figured out the error in the Rosenbrock solver after carefully auditing vs. Numerical Recipes. The Jacobian matrix estimator is not industrial-strength but seems to get the job done. Suggestions for a good way to handle providing an explicit function for the Jacobian matrix are welcome.

I also cleaned up a small fix I made earlier in ode23() since the promotion behavior in min/max has changed to make it unneeded.

Member

pao commented Mar 20, 2012

I'd also like @JonasRauch to take a look at this if he gets a chance. (If I understand right, mentioning him should summon him here, I think that's part of the GitHub "magic Markdown.")

Contributor

JonasRauch commented Mar 20, 2012

Yes, it did send me an email :-)

Looks great to me. The one thing I would change is the interface for the jacobian function. I see two ways:

  1. Have it exposed (like it is now). Then it would make sense to have
    function jacobian(F::Function, x::AbstractVector)
    instead of
    function jacobian(F::Function, t::Number, x::AbstractVector)
    for the same reason we don't have parameters for F in the ODE functions. You can just call jacobian(x->F(t,x), x). If the function is exposed, the interface should be generic and not include the t parameter that is specific to the ODE problem.
  2. Leave the interface the way it is now and move the definition of jacobian into the body of the integrator, making it a private function.

I would vote for 2. atm, because it is a rather crude implementation of a jacobian estimation with a fixed step difference quotient. For the long run, a flexible, robust jacobian estimate should be on the TODO list for numerical methods anyway. Once that is there, this one can be removed.

Member

pao commented Mar 20, 2012

It might be best to do both, so if we are able to take an analytic form of the function we have a common interface. In any case, please do not merge until I'm able to make these changes later today. Thanks, Jonas!

Contributor

JonasRauch commented Mar 20, 2012

That's good, too. But: I think if you have

function ode4s(...)
   function jacobian(F, x)=...
   ...
end

and then another function

`````` jacobian(F,x)=...```
outside, the inner one will always be called. So you could not use an external jacobian inside if it has the same name and signature.

You could leave it to the user to provide a Jacobian function G(t, x) in addition to F:
function oderosenbrock{T}(F::Function, G:::Function, tspan::AbstractVector, x0::AbstractVector{T}, gamma, a, b, c)
and implement the "standard" ODE interface as

    jacobian(...)=...
    oderosenbrock(F, (t,x)->jacobian(F, t, x), x0, gamma, a, b, c)
end
Member

pao commented Mar 20, 2012

I suspect that's where I'll end up, yes. Multiple dispatch FTW.

StefanKarpinski added a commit that referenced this pull request Mar 20, 2012

Merge pull request #615 from pao/topic/ode
Stiff ODE solver + small fix

@StefanKarpinski StefanKarpinski merged commit e862fa1 into JuliaLang:master Mar 20, 2012

Owner

StefanKarpinski commented Mar 20, 2012

GitHub needs to make buttons that do very different things harder to confuse :-/

Member

pao commented Mar 20, 2012

I wanted to fix this before merging...oh well, I'll just do another pull request.

Owner

StefanKarpinski commented Mar 20, 2012

Oh, sorry. I was going to ask, but didn't see any harm in merging it since it only touches ode.jl. Plus, this holds your feet to the fire to fix it ;-)

Owner

StefanKarpinski commented Mar 20, 2012

Not that the fire appears to be necessary at all...

Member

pao commented Mar 21, 2012

It burns, it burns! Working it now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment