[Tutor] hi

Oscar Benjamin oscar.j.benjamin at gmail.com
Mon Aug 12 22:11:21 CEST 2013


On 12 August 2013 20:47, Vick <vick1975 at orange.mu> wrote:
>> -----Original Message-----
>> From: Oscar Benjamin [mailto:oscar.j.benjamin at gmail.com]
>> Sent: Monday, 12 August, 2013 20:53
>>
>> With a smaller stepsize my rk4 integrator can get an error of 1e-15
>> (1e-13%) in 20 milliseconds:
>>
>> $ ./ode.py
>> dt : 9.76562e-04  error : 2.83e-15
>
> [Vick] I have compared your rk4 code with mine and they are virtually the
> same. It is the standard rk4 method. Using only the rk4 method the error is
> 1e-4. However you seem to be getting different result due to the fact that
> you have good programming skills as evidenced by the other defs in your code
> and I really don't understand what the other defs are really doing. The
> solveto and solve defs are certainly doing much to get your result.
>
> Just strip your code of all the superpowers it has got and leave just the
> rk4 method to compute. How much error do you get then? ;)

It's because I'm using double precision and you're using mpmath. The
mpmath library slows your code down by several orders of magnitude for
little gain in accuracy (since the dopri coefficients are not
sufficiently accurate). There are several reasons that rk4 is so
popular and one is that it pretty much has as high an order as is
needed for working in double precision (in most cases).

As I said before mpmath is massive overkill for the requested level of
accuracy. On the other hand if you used a 10th (or higher) order
Adams-Bashforth with mpmath you might start to see the gains from high
order outweigh the costs from mpmath at *really* small error levels.
Of course you have to ask yourself how much time you're prepared to
spend coding up integrators to make the error smaller when you could
really just use scipy's odeint and get on with plotting your results!


Oscar


More information about the Tutor mailing list