[Tutor] hi

Vick vick1975 at orange.mu
Mon Aug 12 17:11:21 CEST 2013


> -----Original Message-----
> From: Oscar Benjamin [mailto:oscar.j.benjamin at gmail.com]
> Sent: Monday, 12 August, 2013 16:08
> To: Vick
> 
> 
> Maybe. It still depends. If you use a higher-order method where it's just
not
> needed you can end up wasting a lot of CPU-cycles. The same goes for using
> multi-precision arithmetic when you could happily work in double
precision.

[Vick] Well I'm using the higher order method for n-body problem. The
problem has four 1st order ODEs to solve for one body; two trivial and two
very significant ones. So for n-body problem you would need (n x 4) 1st
order ODEs to solve.

The equation is set for 2-D coordinates x and y as follows for a 2-body
problem example:

dxdt(3) = vx1
dxdt(4) = vy1
dxdt(1) =  ((G*m2)/(((x1-x2)^2 + (y1-y2)^2)^(3/2)))*(x2-x1)
dxdt(2) = ((G*m2)/(((x1-x2)^2 + (y1-y2)^2)^(3/2)))*(y2-y1)

dxdt(5) = vx2
dxdt(6) = vy2
dxdt(7) = ((G * m1) / (((x1 - x2) ^ 2 + (y1 -y2) ^ 2) ^ (3 / 2))) * (x1 -
x2)
dxdt(8) = ((G * m1) / (((x1 - x2) ^ 2 + (y1 - y2) ^ 2) ^ (3 / 2))) * (y1 -
y2)

>The second is slow but achieves absurdly high levels of accuracy. It
> gives better accuracy than I got with your rungkdp8 (see below) 

[Vick] That's the whole point! I had to wait for about 40 minutes to get
very high accuracy which the method was not designed to obtain.
Whereas the DOPRI8 got the result with 1e-13 error in less than a second.
(Before when I said 2 or less seconds, it is obvious I didn't time them
really; but it just printed the result as soon as I pressed F5)	

> > In that amount of time, the memory was taxed a lot.
> 
> I don't think you really know what you're talking about on this subject.
The
> algorithm works with a constant, very small, amount of memory regardless
of
> stepsize. In fact the rungkdp8 method uses more memory, however, the
> difference is entirely insignificant for such a low-dimensional system. I
think
> you're confusing the fact that the way
> *you* have implemented the algorithm you store every intermediate step in
> the computation. My code only stores the desired output and the current
> value of the intermediate results.

[Vick] Well yeah, maybe there is a confusion on my use of computer memory.
It is obvious that the Runge Kutta 4th order method has only 4 steps and
that it is using very small numbers to compute and whereas the DOPRI8 has 13
steps in which it has to use large numbers to compute. Yes the DOPRI8 uses
more memory than the RK4 in that regard. However the computing speed and the
computer memory the code needs to operate on the algorithm is excessive with
your code as it has to run for 10^(2-6).

> > For the same
> > ODE problem, RK4 should produce an error of 1e-4 in a reasonable
> > amount of time (about 2 secs) and utilizing normal computer memory for
> > its calculations.
> 
> If you modify the check_errors() function in the first script I posted so
that it
> looks like this:
> Then you can run it and get the following:
> 
> $ time ./ode.py
> dt : 7.53390e+01  error : 3.28e-07
> 
> real    0m0.094s
> user    0m0.030s
> sys     0m0.046s
> 
> In other words it takes 90 milliseconds to run the script and get an error
less
> than 1e-6 (or 1e-4% if you prefer). Actually it really just takes 90
milliseconds
> to *initialise* the script. We can use timeit to check how long the actual
> computation takes (with the print line commented out):

[Vick] That's just it. The error margin of RK4 in that instance is just 1e-4
in as you say 90 msec but for about the same amount of time DOPRI8 gets
1e-13.

 
> It takes 168 microseconds to achieve the accuracy you requested, for your
> test problem, using my rk4 integrator. Scipy's odeint is probably faster
(I'll let
> you do the timing if you want to know but bear in mind that it takes a
> significant amount of time to *import* scipy so you should use timeit for
> benchmarking).

[Vick] What I meant was that I don't need to wait for 40,30...10,5 or even 1
minute to get the result. I am not actually measuring the real speed of each
computation by each method. With human perception in view, all I'm saying is
that for less than a second which is a reasonable "waiting" time for a human
computing a running animation of the n-body problem, the error gotten by
each method just demonstrates the superiority of the one over the other. It
has been tested anyway. The DOPRI8 is an 8th order method whereas the RK4 is
a 4th order. For the same amount of time, the latter can offer only 1e-4 of
error for that test problem whereas DOPRI8 offers 1e-13.

> > As I mentioned before I used this problem as a test for accuracy to
> > compare different methods. I use percentage change (x 100) to derive
> > the error for each method.
> 
> I don't understand why you insist on defying convention by using
> percentages. Really your errors should be so small that percent is a
> meaningless unit of measure. Fractional errors are useful because you can
> easily relate them to absolute errors without the pointless cognitive
burden
> of a factor of 100 getting in the way.

[Vick] I don't understand! Normally if a got 12 and b got 12.5, I want to
know by how much a is different from b. So the normal arithmetic if I wanted
the percentage would be: 12*100/12.5 = 96% so the error is 100%-96% =
4%.......doing so with my error calculation we take abs(12.5-12)/12.5 * 100
= 4 whereas with your error calculation it gives (12.5-12)/12.5 = 0.04 . The
percentage is 4% and the fraction of that error is 4/100 = 0.04

Both are valid, and anyway, it is really no big deal as both calculations
are distant by 2 digits only.


> > RK4 in this instance scores about 1e-4 accuracy. With this problem
> > I've tested RKF5, DPRK45, some 6 and 7 order RK methods and the last
> > one is the DOPRI8(7)13.
> 
> Are you testing all the methods with the same step-sizes? If so that's not
a
> fair comparison since they do different amounts of work within each step.
> 

[Vick] That's the whole point for one method to be better than the other in
that for the same amount of time and for the same test problem, in one
single step of their algorithms they are achieving different accuracies.

> 
> I would report a comparison of timings between this code and my decimal
> rk4 but I think it would be unfair as I think there may be something wrong
> with the way your's is implemented. Unless I've misunderstood something
> it's not achieving the proper benefits of the mpmath library since for
some
> reason the error bottoms out at 1e-18. I replaced the bottom of your
script
> with this:
> 
> 
> When I run that I get:
> 
> $ ./vick.py
> dt : 1.0e+00  error : 1.1e-09
> dt : 1.0e-01  error : 9.4e-18
> dt : 1.0e-02  error : 2.3e-18
> dt : 1.0e-03  error : 2.3e-18
> 
> Do you get the same (it could be a difference in the underlying library
etc.)?
> 

[Vick] yeah I get the same results.  The mpf is here to give me more digits
to compute. I think Python would normally give up to 6 digits. When I first
started out with Python I had trouble getting very long decimals. But when I
started to use mpmath, all my problems were solved. So I think it's a habit
to get mpmath involve in any code I make in python which involves scientific
or advanced math calculations. If you increase the dps level the number of
decimals the code will print will become greater.

Vick




More information about the Tutor mailing list