[Tutor] hi
Oscar Benjamin
oscar.j.benjamin at gmail.com
Mon Aug 12 18:42:05 CEST 2013
On 12 August 2013 16:32, Vick <vick1975 at orange.mu> wrote:
>> From: Oscar Benjamin [mailto:oscar.j.benjamin at gmail.com]
>
>> The rational approximations are intended to
>> be good enough for use with double precision but are insufficient for the
> 30
>> dps setting that you are using. The reason given by the authors for
> publishing
>> approximate coefficients is that (in 1981) it was deemed computationally
>> intractable to compute the exact rational coefficients. However that may
> not
>> still be the case with modern computing power. I doubt that I'll attempt
> this
>> any time soon but if you're able (and bothered) to do that I'd love to
> know
>> what the coefficients would be.
>>
> [Vick] Yeah they are rational coefficients and I got them from the google
> book satellites orbits which I believe I posted the links for it previously.
>
> The problem is how are the coefficients derived? I've never seen an
> algorithm to compute the coefficients of an integrator method. Anyway it
> would surely involve complicated coding to start an algorithm to compute
> them. If I can get it somewhere, I wouldn't mind doing the DOPRI8
> computations.
Are you able to access this paper (link at the bottom)?
P.J. Prince, J.R. Dormand, High order embedded Runge-Kutta
formulae, Journal of Computational and Applied Mathematics, Volume 7,
Issue 1, March 1981, Pages 67-75, ISSN 0377-0427,
http://dx.doi.org/10.1016/0771-050X(81)90010-3.
(http://www.sciencedirect.com/science/article/pii/0771050X81900103)
I believe this is the original source for the coefficients you are
using. They describe constructing a system of linear rational
equations for the coefficients of the method. I think that the
equations must be satisfied for any 13-stage explicit RK method that
has an 8th-order with embedded 7th-order method. Solving the equations
leads to a solution space with 10 degrees of freedom. You could use
e.g. sympy for this step of the calculation.
They then use these 10 degrees of freedom to try and satisfy some
additional conditions that they describe in the paper such as
minimising the coefficient of the leading error term. It is this last
step that they claim was too computationally expensive, so they did it
in 24 d.p numerical precision and then simplified the results to
fractions that are accurate to 18 d.p. suitable for use in double
precision computation. I haven't studied the equations in detail to
see why they might be expensive to solve. If it is possible you would
probably want to use the fractions module for the arithmetic at this
step.
As a simpler example, imagine constructing the 4th order RK method
which has Butcher table:
0 |
1/2 | 1/2
1/2 | 0 1/2
1 | 0 0 1
----|----------------
| 1/6 1/3 1/3 1/6
Essentially there are a number of constraints that must be satisfied
for the method to be sensible and to have the appropriate order:
1) The top-left entry is always zero.
2) The sum of any row in the upper-right part is equal to number at
the left of the row (3 equations)
3) The coefficients in the Taylor expansion for the first four terms
in the error must be zero (4 equations)
4) Possibly more ...
This gives us 8 equations for the 15 coefficients:
a1 |
a2 | b21
a3 | b31 b32
a4 | b41 b42 b43
---|-----------------
| c1 c2 c3 c4
If the equations have full rank then their solution has 7 degrees of
freedom. More equations can easily be constructed by imposing
additional conditions on how you want the resulting table to look.
Standard rk4 is one way to satisfy these conditions. Another way is
the 3/8 rule:
0 |
1/3 | 1/3
2/3 | -1/3 1
1 | 1 -1 1
----|--------------------
| 1/8 3/8 3/8 1/8
Which you prefer is partly arbitrary. The second will usually have
smaller error terms but the first is very efficient to implement and
code because it has a few zeros in it. Essentially the authors in the
Dormand and Prince paper have chosen some complex additional
conditions for their 7/8 table and somehow tried to find the
coefficients that satisfy them.
Much as I would love to see those coefficients, for your own sake it
would be a lot easier to just use an Adams-Bashforth integrator.
Here's a stripped down script that can compute the exact coefficients
for any arbitrary order:
#!/usr/bin/env python
#
# Compute coefficients for Adams-Moulton integrator.
#
from math import factorial
from operator import neg
from itertools import dropwhile
from fractions import Fraction
class Polynomial(tuple):
def __new__(typ, coefficients):
coeffs = map(Fraction, coefficients)
coeffs = dropwhile(lambda c: c==0, coeffs)
return tuple.__new__(typ, coeffs)
@property
def order(self):
return len(self) - 1
def evaluate(self, x, denomlimit=None):
total, = self[:1] or [0]
for c in self[1:]:
total = total * x + c
if denomlimit is not None:
total = total.limit_denominator(denomlimit)
return total
def integrate(self, a, b):
coeffs = [c / (self.order - n + 1) for n, c in enumerate(self)]
anti = Polynomial(coeffs + [0])
return anti.evaluate(b) - anti.evaluate(a)
def __mul__(p1, p2):
'''
>>> from eventphys.poly import Polynomial
>>> Polynomial(['2', '-1/2']) * Polynomial([3, 2, '1/2'])
Polynomial(['6', '5/2', '0', '-1/4'])
'''
p3 = [0] * (len(p1) + len(p2) - 1)
for n1, c1 in enumerate(p1):
for n2, c2 in enumerate(p2):
p3[n1 + n2] += c1 * c2
return Polynomial(p3)
def ab_coeffs(nsteps):
b = []
N = nsteps
for j in range(N):
p = Polynomial([1])
for i in range(N):
if i != j:
p *= Polynomial([1, i])
a = Fraction((-1)**j, factorial(j) * factorial(N-j-1))
b.insert(0, a * p.integrate(0, 1))
return b
def main(nsteps):
coeffs = ab_coeffs(int(nsteps))
print('Coefficients for Adams-Moulton with %s steps' % nsteps)
for c in coeffs:
print(c)
if __name__ == "__main__":
import sys
main(*sys.argv[1:])
Example output:
$ ./am.py 5
Coefficients for Adams-Moulton with 5 steps
251/720
-637/360
109/30
-1387/360
1901/720
Compare with
http://en.wikipedia.org/wiki/Linear_multistep_method#Adams.E2.80.93Bashforth_methods
Oscar
More information about the Tutor
mailing list