On Tue, Oct 2, 2012 at 5:12 PM, Ondřej Čertík <[email protected]> wrote:
[...]
> For example the difficult part of I_{9/2}(x) is:
>
> In [1]: r = (105/x**4 + 45/x**2 + 1)*sinh(x) - (105/x**3 + 10/x)*cosh(x)
[...]
> So after this painful analysis, we have found:
>
> [0, 0.4] use [5]
> [4, oo] use [1]
>
>
> Now we can expand the function around x=1, and repeat the analysis
> above until we cover the whole real axis. So first of all, this should
> be automated.
> But then the series expansion is *not* the best approximation, because
> the series is very (too much) accurate around x=0, and barely accurate
> around x=0.4.
> A better approximation is to use a so called Chebyshev approximation,
> which gives uniform accuracy (the end result is that less terms are
> needed).
>
> Finally, even better than just using one series, it's better to use a
> rational function, which is a fraction of two polynomials. This seems
> to be the most effective.
I've used MiniMaxApproximation in Mathematica for the function
I_{4+1/2}(x) and using trial and error
I divided the real axis into the intervals (0, 0.2), (0.2, 1.7), (1.7,
4), (4, 10), (10, 20), (20, oo),
in the first I use the series from sympy, second, third and fourth
uses rational approximation
from Mathematica, fifth uses the exact analytic formula using
sinh/cosh, and the last one
uses the fact that sinh(x)/exp(x) = 1/2 for x from (20, oo).
The Python code and graphs are here:
https://gist.github.com/3837394
You can see from the error plot, that the rational approximation is
accurate to machine precision
(compared to the "exact" answer from mpmath), pretty much uniformly for any "x".
A possible improvement is to use the Horner's rule to evaluate the
polynomials, but otherwise
I think this is pretty much as fast as it can be, the only way to make
it faster is to put there
more intervals and lower the polynomial order.
So this is the solution that I am looking for. It'd be nice to be able
to do the same thing using sympy or mpmath.
Ondrej
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
