Convergence acceleration / extrapolation methods for series and
sequences.
References:
Carl M. Bender & Steven A. Orszag, “Advanced Mathematical Methods for
Scientists and Engineers: Asymptotic Methods and Perturbation Theory”,
Springer 1999. (Shanks transformation: pp. 368-375, Richardson
extrapolation: pp. 375-377.)
modelparameters.sympy.series.acceleration.richardson(A, k, n, N)[source]¶
Calculate an approximation for lim k->oo A(k) using Richardson
extrapolation with the terms A(n), A(n+1), …, A(n+N+1).
Choosing N ~= 2*n often gives good results.
A simple example is to calculate exp(1) using the limit definition.
This limit converges slowly; n = 100 only produces two accurate
digits:
modelparameters.sympy.series.acceleration.shanks(A, k, n, m=1)[source]¶
Calculate an approximation for lim k->oo A(k) using the n-term Shanks
transformation S(A)(n). With m > 1, calculate the m-fold recursive
Shanks transformation S(S(…S(A)…))(n).
The Shanks transformation is useful for summing Taylor series that
converge slowly near a pole or singularity, e.g. for log(2):
Return a generator for consecutive Pade approximants for a series.
It can also be used for computing the rational generating function of a
series when possible, since the last approximant returned by the generator
will be the generating function (if any).
The input list can contain more complex expressions than integer or rational
numbers; symbols may also be involved in the computation. An example below
show how to compute the generating function of the whole Pascal triangle.
The generator can be asked to apply the sympy.simplify function on each
generated term, which will make the computation slower; however it may be
useful when symbols are involved in the expressions.
x0 (number, optional) – Point to perform series expansion about. Default is 0.
dir ({1, -1, '+', '-'}, optional) – If dir is 1 or ‘+’ the series is calculated from the right and
for -1 or ‘-’ the series is calculated from the left. For smooth
functions this flag will not alter the results. Default is 1.
hyper ({True, False}, optional) – Set hyper to False to skip the hypergeometric algorithm.
By default it is set to False.
order (int, optional) – Order of the derivative of f, Default is 4.
rational ({True, False}, optional) – Set rational to False to skip rational algorithm. By default it is set
to True.
full ({True, False}, optional) – Set full to True to increase the range of rational algorithm.
See rational_algorithm() for details. By default it is set to
False.
Returns the formal series expansion of f around x=x0
with respect to x in the form of a FormalPowerSeries object.
Formal Power Series is represented using an explicit formula
computed using different algorithms.
See compute_fps() for the more details regarding the computation
of formula.
Parameters:
x (Symbol, optional) – If x is None and f is univariate, the univariate symbols will be
supplied, otherwise an error will be raised.
x0 (number, optional) – Point to perform series expansion about. Default is 0.
dir ({1, -1, '+', '-'}, optional) – If dir is 1 or ‘+’ the series is calculated from the right and
for -1 or ‘-’ the series is calculated from the left. For smooth
functions this flag will not alter the results. Default is 1.
hyper ({True, False}, optional) – Set hyper to False to skip the hypergeometric algorithm.
By default it is set to False.
order (int, optional) – Order of the derivative of f, Default is 4.
rational ({True, False}, optional) – Set rational to False to skip rational algorithm. By default it is set
to True.
full ({True, False}, optional) – Set full to True to increase the range of rational algorithm.
See rational_algorithm() for details. By default it is set to
False.
Rational algorithm for computing
formula of coefficients of Formal Power Series
of a function.
Applicable when f(x) or some derivative of f(x)
is a rational function in x.
rational_algorithm() uses apart() function for partial fraction
decomposition. apart() by default uses ‘undetermined coefficients
method’. By setting full=True, ‘Bronstein’s algorithm’ can be used
instead.
Looks for derivative of a function up to 4’th order (by default).
This can be overriden using order option.
By setting full=True, range of admissible functions to be solved using
rational_algorithm can be increased. This option should be used
carefully as it can signifcantly slow down the computation as doit is
performed on the RootSum object returned by the apart function.
Use full=False whenever possible.
Return \(\sigma\)-approximation of Fourier series with respect
to order n.
Sigma approximation adjusts a Fourier summation to eliminate the Gibbs
phenomenon which would otherwise occur at discontinuities.
A sigma-approximated summation for a Fourier series of a T-periodical
function can be written as
where \(a_0, a_k, b_k, k=1,\ldots,{m-1}\) are standard Fourier
series coefficients and
\(\operatorname{sinc} \Bigl( \frac{k}{m} \Bigr)\) is a Lanczos
\(\sigma\) factor (expressed in terms of normalized
\(\operatorname{sinc}\) function).
Parameters:
n (int) – Highest order of the terms taken into account in approximation.
Returns:
Sigma approximation of function expanded into Fourier series.
The behaviour of
sigma_approximation()
is different from truncate()
- it takes all nonzero terms of degree smaller than n, rather than
first n nonzero ones.
Implemented according to the PhD thesis
http://www.cybertester.com/data/gruntz.pdf, which contains very thorough
descriptions of the algorithm including many examples. We summarize here
the gist of it.
All functions are sorted according to how rapidly varying they are at
infinity using the following rules. Any two functions f and g can be
compared using the properties of L:
So we can divide all the functions into comparability classes (x and x^2
belong to one class, exp(x) and exp(-x) belong to some other class). In
principle, we could compare any two functions, but in our algorithm, we
don’t compare anything below the class 2~3~-5 (for example log(x) is
below this), so we set 2~3~-5 as the lowest comparability class.
Given the function f, we find the list of most rapidly varying (mrv set)
subexpressions of it. This list belongs to the same comparability class.
Let’s say it is {exp(x), exp(2x)}. Using the rule f ~ 1/f we find an
element “w” (either from the list or a new one) from the same
comparability class which goes to zero at infinity. In our example we
set w=exp(-x) (but we could also set w=exp(-2x) or w=exp(-3x) …). We
rewrite the mrv set using w, in our case {1/w, 1/w^2}, and substitute it
into f. Then we expand f into a series in w:
We need to recursively compute limits at several places of the algorithm, but
as is shown in the PhD thesis, it always finishes.
Important functions from the implementation:
compare(a, b, x) compares “a” and “b” by computing the limit L.
mrv(e, x) returns list of most rapidly varying (mrv) subexpressions of “e”
rewrite(e, Omega, x, wsym) rewrites “e” in terms of w
leadterm(f, x) returns the lowest power term in the series of f
mrv_leadterm(e, x) returns the lead term (c0, e0) for e
limitinf(e, x) computes lim e (for x->oo)
limit(e, z, z0) computes any limit by converting it to the case x->oo
All the functions are really simple and straightforward except
rewrite(), which is the most difficult/complex part of the algorithm.
When the algorithm fails, the bugs are usually in the series expansion
(i.e. in SymPy) or in rewrite.
This code is almost exact rewrite of the Maple code inside the Gruntz
thesis.
Because the gruntz algorithm is highly recursive, it’s difficult to
figure out what went wrong inside a debugger. Instead, turn on nice
debug prints by defining the environment variable SYMPY_DEBUG. For
example:
Stores (expr, dummy) pairs, and how to rewrite expr-s.
The gruntz algorithm needs to rewrite certain expressions in term of a new
variable w. We cannot use subs, because it is just too smart for us. For
example:
So we do it the hard way and keep track of all the things we potentially
want to substitute by dummy variables. Consider the expression:
exp(x-exp(-x))+exp(x)+x.
The mrv set is {exp(x), exp(-x), exp(x - exp(-x))}.
We introduce corresponding dummy variables d1, d2, d3 and rewrite:
d3+d1+x.
This class first of all keeps track of the mapping expr->variable, i.e.
will at this stage be a dictionary:
{exp(x):d1,exp(-x):d2,exp(x-exp(-x)):d3}.
[It turns out to be more convenient this way round.]
But sometimes expressions in the mrv set have other expressions from the
mrv set as subexpressions, and we need to keep track of that as well. In
this case, d3 is really exp(x - d2), so rewrites at this stage is:
{d3:exp(x-d2)}.
The function rewrite uses all this information to correctly rewrite our
expression in terms of w. In this case w can be choosen to be exp(-x),
i.e. d2. The correct rewriting then is:
Compute the limit of e(z) at the point z0 using the Gruntz algorithm.
z0 can be any expression, including oo and -oo.
For dir=”+” (default) it calculates the limit from the right
(z->z0+) and for dir=”-” the limit from the left (z->z0-). For infinite z0
(oo or -oo), the dir argument doesn’t matter.
This algorithm is fully described in the module docstring in the gruntz.py
file. It relies heavily on the series expansion. Most frequently, gruntz()
is only used if the faster limit() function (which uses heuristics) fails.
Computes the maximum of two sets of expressions f and g, which
are in the same comparability class, i.e. mrv_max1() compares (two elements of)
f and g and returns the set, which is in the higher comparability class
of the union of both, if they have the same order of variation.
Also returns exps, with the appropriate substitutions made.
Computes the maximum of two sets of expressions f and g, which
are in the same comparability class, i.e. max() compares (two elements of)
f and g and returns either (f, expsf) [if f is larger], (g, expsg)
[if g is larger] or (union, expsboth) [if f, g are of the same class].
The result of this function is currently undefined if e changes sign
arbitarily often for arbitrarily large x (e.g. sin(x)).
Note that this returns zero only if e is constantly zero
for x sufficiently large. [If e is constant, of course, this is just
the same thing as the sign of e.]
Takes as input a polynomial expression and the variable used to construct
it and returns the difference between function’s value when the input is
incremented to 1 and the original function value. If you want an increment
other than one supply it as a third argument.
Takes as input a Sum instance and returns the difference between the sum
with the upper index incremented by 1 and the original sum. For example,
if S(n) is a sum, then finite_diff_kauers will return S(n + 1) - S(n).
Return from the atoms of self those which are free symbols.
For most expressions, all symbols are free symbols. For some classes
this is not true. e.g. Integrals use Symbols for the dummy variables
which are bound variables, so Integral has a method to return all
symbols except those. Derivative keeps track of symbols with respect
to which it will perform a derivative; those are
bound variables, too, so it has its own free_symbols method.
Any other method that uses bound variables should implement a
free_symbols method.
For dir=”+” (default) it calculates the limit from the right
(z->z0+) and for dir=”-” the limit from the left (z->z0-). For infinite
z0 (oo or -oo), the dir argument is determined from the direction
of the infinity (i.e., dir=”-” for oo).
First we try some heuristics for easy and frequent cases like “x”, “1/x”,
“x**2” and similar, so that it’s fast. For all other cases, we use the
Gruntz algorithm (see the gruntz() function).
trials (int, optional) – The algorithm is highly recursive. trials is a safeguard from
infinite recursion incase limit is not easily computed by the
algorithm. Try increasing trials if the algorithm returns None.
Terms (Admissible) –
================ –
functions (The terms should be built from rational) –
sums (indefinite) –
:param :
:param and indefinite products over an indeterminate n. A term is admissible:
:param if the scope of all product quantifiers are asymptotically positive.:
:param Every admissible term is asymptoticically monotonous.:
The order of a function characterizes the function based on the limiting
behavior of the function as it goes to some limit. Only taking the limit
point to be a number is currently supported. This is expressed in
big O notation [1]_.
The formal definition for the order of a function g(x) about a point a
is such that g(x) = O(f(x)) as x rightarrow a if and only if for any
delta > 0 there exists a M > 0 such that |g(x)| leq M|f(x)| for
|x-a| < delta. This is equivalent to lim_{x rightarrow a}
sup |g(x)/f(x)| < infty.
Let’s illustrate it on the following example by taking the expansion of
sin(x) about 0:
\[\sin(x) = x - x^3/3! + O(x^5)\]
where in this case O(x^5) = x^5/5! - x^7/7! + cdots. By the definition
of O, for any delta > 0 there is an M such that:
As it is usually used, the order of a function can be intuitively thought
of representing all terms of powers greater than the one specified. For
example, O(x^3) corresponds to any terms proportional to x^3,
x^4,ldots and any higher power. For a polynomial, this leaves terms
proportional to x^2, x and constants.
Return True if expr belongs to Order(self.expr, *self.variables).
Return False if self belongs to expr.
Return None if the inclusion relation cannot be determined
(e.g. when self and expr have different symbols).
Return from the atoms of self those which are free symbols.
For most expressions, all symbols are free symbols. For some classes
this is not true. e.g. Integrals use Symbols for the dummy variables
which are bound variables, so Integral has a method to return all
symbols except those. Derivative keeps track of symbols with respect
to which it will perform a derivative; those are
bound variables, too, so it has its own free_symbols method.
Any other method that uses bound variables should implement a
free_symbols method.
Finds the shortest linear recurrence that satisfies the first n
terms of sequence of order leq n/2 if possible.
If d is specified, find shortest linear recurrence of order
leq min(d, n/2) if possible.
Returns list of coefficients [b(1),b(2),...] corresponding to the
recurrence relation x(n)=b(1)*x(n-1)+b(2)*x(n-2)+...
Returns [] if no recurrence is found.
If gfvar is specified, also returns ordinary generating function as a
function of gfvar.
Return from the atoms of self those which are free symbols.
For most expressions, all symbols are free symbols. For some classes
this is not true. e.g. Integrals use Symbols for the dummy variables
which are bound variables, so Integral has a method to return all
symbols except those. Derivative keeps track of symbols with respect
to which it will perform a derivative; those are
bound variables, too, so it has its own free_symbols method.
Any other method that uses bound variables should implement a
free_symbols method.
The order of a function characterizes the function based on the limiting
behavior of the function as it goes to some limit. Only taking the limit
point to be a number is currently supported. This is expressed in
big O notation [1]_.
The formal definition for the order of a function g(x) about a point a
is such that g(x) = O(f(x)) as x rightarrow a if and only if for any
delta > 0 there exists a M > 0 such that |g(x)| leq M|f(x)| for
|x-a| < delta. This is equivalent to lim_{x rightarrow a}
sup |g(x)/f(x)| < infty.
Let’s illustrate it on the following example by taking the expansion of
sin(x) about 0:
\[\sin(x) = x - x^3/3! + O(x^5)\]
where in this case O(x^5) = x^5/5! - x^7/7! + cdots. By the definition
of O, for any delta > 0 there is an M such that:
As it is usually used, the order of a function can be intuitively thought
of representing all terms of powers greater than the one specified. For
example, O(x^3) corresponds to any terms proportional to x^3,
x^4,ldots and any higher power. For a polynomial, this leaves terms
proportional to x^2, x and constants.
Return True if expr belongs to Order(self.expr, *self.variables).
Return False if self belongs to expr.
Return None if the inclusion relation cannot be determined
(e.g. when self and expr have different symbols).
Return from the atoms of self those which are free symbols.
For most expressions, all symbols are free symbols. For some classes
this is not true. e.g. Integrals use Symbols for the dummy variables
which are bound variables, so Integral has a method to return all
symbols except those. Derivative keeps track of symbols with respect
to which it will perform a derivative; those are
bound variables, too, so it has its own free_symbols method.
Any other method that uses bound variables should implement a
free_symbols method.
Returns the formal series expansion of f around x=x0
with respect to x in the form of a FormalPowerSeries object.
Formal Power Series is represented using an explicit formula
computed using different algorithms.
See compute_fps() for the more details regarding the computation
of formula.
Parameters:
x (Symbol, optional) – If x is None and f is univariate, the univariate symbols will be
supplied, otherwise an error will be raised.
x0 (number, optional) – Point to perform series expansion about. Default is 0.
dir ({1, -1, '+', '-'}, optional) – If dir is 1 or ‘+’ the series is calculated from the right and
for -1 or ‘-’ the series is calculated from the left. For smooth
functions this flag will not alter the results. Default is 1.
hyper ({True, False}, optional) – Set hyper to False to skip the hypergeometric algorithm.
By default it is set to False.
order (int, optional) – Order of the derivative of f, Default is 4.
rational ({True, False}, optional) – Set rational to False to skip rational algorithm. By default it is set
to True.
full ({True, False}, optional) – Set full to True to increase the range of rational algorithm.
See rational_algorithm() for details. By default it is set to
False.
Compute the limit of e(z) at the point z0 using the Gruntz algorithm.
z0 can be any expression, including oo and -oo.
For dir=”+” (default) it calculates the limit from the right
(z->z0+) and for dir=”-” the limit from the left (z->z0-). For infinite z0
(oo or -oo), the dir argument doesn’t matter.
This algorithm is fully described in the module docstring in the gruntz.py
file. It relies heavily on the series expansion. Most frequently, gruntz()
is only used if the faster limit() function (which uses heuristics) fails.
For dir=”+” (default) it calculates the limit from the right
(z->z0+) and for dir=”-” the limit from the left (z->z0-). For infinite
z0 (oo or -oo), the dir argument is determined from the direction
of the infinity (i.e., dir=”-” for oo).
First we try some heuristics for easy and frequent cases like “x”, “1/x”,
“x**2” and similar, so that it’s fast. For all other cases, we use the
Gruntz algorithm (see the gruntz() function).
trials (int, optional) – The algorithm is highly recursive. trials is a safeguard from
infinite recursion incase limit is not easily computed by the
algorithm. Try increasing trials if the algorithm returns None.
Terms (Admissible) –
================ –
functions (The terms should be built from rational) –
sums (indefinite) –
:param :
:param and indefinite products over an indeterminate n. A term is admissible:
:param if the scope of all product quantifiers are asymptotically positive.:
:param Every admissible term is asymptoticically monotonous.: