Bug in implementation of FSAL in oderkf. #15
Conversation
|
Sure, go for it. (Not to make excuses, but the algorithms I added are all taken with a PowerPoint presentation of lecture notes as the only real reference; the existence of such a bug is unsurprising.) At any rate, I can hit the button, but if you're looking for a review you might want to ping someone who knows the subject area better. |
|
@pao : It's really not your fault; the octave implementation already had this bug (and I only figured it out by directly comparing The fix is IMO straight-forward, but let's see if @tlycken or someone else has any objections? |
|
It looks straightforward and not obviously wrong. I just don't feel too comfortable in this code anymore. |
|
You mentioned that the Octave implementation had this bug; does it still? We should probably let them know. http://savannah.gnu.org/bugs/?group=octave |
|
Well, the implementation at http://users.powernet.co.uk/kienzle/octave/matcompat/scripts/ode_v1.11/ode45.m (version v1.11) has this bug. I don't know if it ever was an official part of Octave? However, the author (Marc Compere) has a new version (v1.15) here https://sites.google.com/site/comperem/home/ode_solvers where the error has been corrected, so I guess the Octaveians (?) should know about it. On an unrelated issue: Do you know how this licensing business works? Since |
|
Yeah, I had completely forgotten where this came from; this was a refactoring of Regarding the licensing question, I'm pretty sure (but IANAL) we're bound by GPL unless we were to replace the lot with a from-scratch implementation based on the literature. Or I believe we could also ask Octave to allow us to distribute our derived version of Octave's implementation of ode45 under a more permissive license. |
| @@ -197,25 +197,13 @@ function oderkf{T}(F::Function, tspan::AbstractVector, x0::AbstractVector{T}, a, | |||
| xout = x.' # first output solution | |||
|
|
|||
| k = zeros(eltype(x), length(c), length(x)) | |||
|
|
|||
| k[1,:] = F(t,x) # first stage | |||
|
|
|||
There was a problem hiding this comment.
Whitespace error - there's a new tab at the end of this line.
|
Except for the whitespace errors I noted above, I don't see any problems with this. (And I don't know if we care enough about whitespace errors to do anything about them - just thought it'd be better to make a concious decision not to care, than to just forget about it...) However, I don't know how I suddenly became the expert on adaptive-step RK methods... :P |
You know, it is the usual story: one day you are sitting in your couch and enjoying life, next day AP calls and asks about your expert opinion on the world crisis in ... |
|
@pao: I think this is good to merge ... |
Bug in implementation of FSAL in oderkf.
Looking into eventually using
oderkfalso for the Bogacki–Shampine method (akaode23), I found that the "First Same As Last " (FSAL) property of the Dormand-Prince pair is not correctly implemented. The error was already present in the "source" ode45, but has been corrected.Essentially,
k[1,:] = k[end,:]should only be done if the integration step was accepted. I also took the liberty to remove an unnecessary check ofdelta==0.(cc: @pao)