numericalnim
numericalnim copied to clipboard
SciNim
add vern76 to ode integrators
Hello.
Another algorithm, this time it is Verners 7th order algorithm. Tests work fine, only the vector high precision one takes a little longer. I think that is just the high order of the algo, however you might want to take a look at it.
I modified the description for Verner's 6th order to include the "most efficient" description because on the publication Page there is also a "most robust" version. According to a paper by Verner the "most efficient" ones are to be preferred.
Also sorry for bothering you with all these pull requests. ^^
Bother?! 😍 Your pull requests are amazing!! Thanks to you NumericalNim has a much higher performance now than it did have just a week ago.
The vectors are really slow (I made them 🤷) so it is no big surprise. I might even when I think about it have a different way of measuring error for it compared to Arrayman er's Tensors 😅
The code looks good except that the order is set to 6.0 instead of 7.0 😊 that might have to do with the slowing. And for some reasong I haven't past the order variable to the step function 😅 so the 1/5 you see there in the adjustment of timestep is actually 1/order so it should be 1/7 in this case
That's the case for Vern65 as well 😅 I should have made it more clear. So it entirely my fault.
Great initiative to write even more clearly which algorithm it is 😉👍
Oh, I will fix these things when I am back home which should be in about 1.5-2h. :grinning:
I might look into the vectors implementation too, maybe we can speed them up a bit. :wink:
Great! No hurry at all 😊👍
I can't stop you, but I don't know how to optimize them much more. They are just for-loops iterating over the vector's elements 😅🤷
Nice! Do you see any difference in performance after the changes? 😄 (I'm not currently at a computer)
Not for the small test I do. But that is really just a test wether the algorithm works.
However Vern7 is about ten times slower than Vern6. 🤔 Will look into this.
I was thinking about leaving out the calculations when a_ij = 0.0 during calculation of the ks. I am not sure whether the compiler kicks those out or not?
The last commit is just a little cosmetic fix for the documentation comments of Vern76 and Vern65. They still referred to Dopri54 because I used that as a template.
Is Vern7 consistently slower than Vern6 or is it only for bigger tol? It does sound suspicious otherwise. 🤔 I'll do some testing as well tomorrow. Nothing is like finding what tiny detail is wrecking an algorithm 😂
It wouldn't hurt removing the 0.0's but I sure hope the compiler is smart enough to optimize them away (at least if a_ij is a const).
I can confirm that Vern7 is really slow.... Took it 16 second while Vern6 & Tsit5 took less than my system could measure 😅 so there is probably something wrong in the time step calculations because the error is there nothing wrong with
I did some more tests and Vern7 gets slower by increasing the timespan. It becomes slower by factors > 100 hundred compared to Vern6.
Decreasing the default dtMin back to 1e-8 makes the algorithms I tested (Dopri54, Tsit54, Vern65) run 10-20% faster for shorter calculations. For longer ones there I see no effect but we might consider changing the default.
Yeah we should change it back, the reason we changed it last time wasn't even a valid reason 😅 it was just a misunderstanding in the Tsit5.
If it gets slower for longer timespans it definitely feels like it has to do with the adaptive timestep as it plays a larger role if it's faulty on longer time spand
I did the same kind of test as I did before measuring how many steps each method took and what error it thought it had (error is what the method self thought the error was):
- Vern6:
Steps: 103 Total error: 0.0002288001802899231 Error/step: 2.221360973688574e-006 - Tsit54:
Steps: 103 Total error: 0.008341642556670692 Error/step: 8.09868209385504e-005 - Bs32:
Steps: 536 Total error: 0.3926441066888433 Error/step: 0.0007325449751657524 - Vern7:
Steps: 5428350 Total error: 40231.00641081012 Error/step: 0.007411277167244212
This confirms that the time step is chosen wrongly, it wouldn't take 10^4 times too many steps otherwise.... 😅 So next up is checking that all the coefficients are correct. It is really easy to copy the wrong numbers 👍😄
I just checked the coefficients twice, they are correct. Made that mistake when implementing Tsit54 and Vern6 but then it was obvious because of big errors. Vern7 does not create those so I think the coefficients are OK.
If you test to return yLow instead of yHigh, do you get acceptable errors if you force the minDt pretty high (1e-1) if you compare it to now when we return yHigh. The adaptive time step shouldn't play much of a difference then so we should be able to see if the methods are correct or not
Returning yLow increases the error significantly and also increases time needed. Setting minDt to high values did not change much however.
That's really odd as well... 😵 The time shouldn't change between the to as we have already done all the time step adjustments when we return yLow. But that the error significantly rises is concerning...
Tomorrow I will check the coefficients for the y's against JuliaDiffEq to see if I can see any difference
Just checked yLow once more, Dopri54 and Vern65 show no difference, tsit54 and Vern76 need double time in my test.
Ahhhh! 😵😱😮 This is so confusing. Can't see how it would affect the time to switch. But Tsit54 isn't strange as yLow isn't what we mean by yLow 😅 Do you have sny idea why it behaves so strange?
Not really. I just drop in replaced the coefficients of "most efficient" Vern7 with those of "most robust" Vern7, and now it works flawlessly! :flushed: So there is probably still a mistake in the coefficients.
Strange, but great that it works 😊👍 did you change all of the coefficients or was there common ones?
There are like 5-10 common ones, when I omit the 0.0s. I am currently checking the coeffs of "most efficient" once again, but rigorous this time.
I think we can go with the "most robust" if it works. Can be good to have an integrator that is robust 😄👍😊
I feel like the "most effective" is getting very personal 😂 can be a dangerous rabbit hole to step into