Exponential Solvers

Shifted Exponentials

The shifted exponentials, usually called φ-functions, are defined by the following series:

\[φ_{\ell}(z) = ∑_{k=0}^{∞} \frac{z^k}{(\ell+k)!} \qquad \ell \in\NN\]

Notice in particular that

\[φ_0(z) = \ee^z\]

We observe the following recursion relation

\[φ_{\ell}(z) = \frac{1}{\ell!} + z ∑_{k=1}^{∞}\frac{x^{k-1}}{(\ell+k)!} = \frac{1}{\ell!} + z φ_{\ell+1}\]

This allows to prove the useful identity, valid for \(\ell≥1\):

\[φ_{\ell}(z) = \int_0^1 \ee^{z (1-x)} \frac{x^{\ell-1}}{(\ell - 1)!} \dd x\]

Indeed, by integration by parts, the recursion relation is the same as above, and for \(\ell=1\) one has

\[φ_1(z) = \frac{\ee^z-1}{z}\]

Padé Approximations

One computes Padé approximations of the form

\[φ_{\ell}(z) = \frac{N(z)}{D(z)}\]

\[d^{\ell} := \frac{d!}{(2d-\ell)!}\]

\[D_j^\ell := d^{\ell} (-1)^j \frac{ (2d + \ell -j)!}{j! (d-j)!}\]

So one may start with:

\[D_0^0 = 1\]

and compute for \(0≤j<d\):

\[D_{j+1}^0 = \frac{-(d-j)}{(2d -j)(j+1)} D_{j}^{0}\]

then, for \(0≤ \ell < \ell_{\max}\):

\[D_{j}^{\ell+1} = (2d-\ell)(2d+\ell+1-j) D_{j}^{\ell}\]

\[C_j := \frac{1}{j!}\]

\[N^{\ell} = C^{\ell} \star D^{\ell}\]

phi_pade – Phi Padé

Computation of φ functions using Padé approximations and Scaling and Squaring.

Some formulae are taken from the Expint documentation, of Håvard Berland, Bård Skaflestad and Will Wright.

class odelab.phi_pade.Phi(k, d=6)

Main class to compute the \(φ_l\) functions. The simplest way to define those functions is by the formula:

\[φ_{\ell}(x) = ∑_{k=0}^{∞} \frac{x^k}{(\ell+k)!}\]

Usage is as follows:

phi = Phi(k,d)
result = phi(M)

where M is a square array. The variable result is a list of all the values of \(φ_{k}(M)\) for \(0≤k≤l\).

compute_Pade()

Compute the Padé approximations of order \(d\) of \(φ_l\), for \(0 ≤ l ≤ k\).

The goal is to produce an approximation of the form:

\[φ_{\ell} = \frac{N}{D}\]

where \(N\) and \(D\) are polynomials. The formula for \(D\) is first computed recursively using the following recursion relations:

\[\begin{split}D_0^0 = 1\\ D_{j+1}^0 = \frac{-(d-j)}{(2d -j)(j+1)} D_{j}^0\\ D_{j}^{\ell+1} = (2d-\ell)(2d+\ell+1-j) D_{j}^{\ell}\end{split}\]

Then, considering the inverse factorial series:

\[C_j := \frac{1}{j!}\]

The numerator \(N\) is now computed by:

\[N = D*C\]

eval_pade(z, s=None)

Evaluate \(φ_l(z)\) using the Padé approximation.

phi_square(l)

Formula for squaring phi_l from existing phi_k for k≤l, taken from the Expint documentation.

class odelab.phi_pade.Polynomial(coeffs)

Polynomial class used in the Padé approximation. Usage is as follows:

p = Polynomial([4.,2.,1.,...]) # polynomial 4 + 2x + x**2 + ...
Z = Polynomial.exponents(z, 3) # compute z**2 and z**3
p(Z) # value of the polynomial p(z)

The evaluation at a matrix value is computed using the Paterson and Stockmeyer method (see [Golub] § 11.2.4). The polynomial is split into chunks of size s.

[Golub]

Golub, G.H. and van Loan, C.F., Matrix Computations, 3rd ed. isbn:9780801854149

classmethod exponents(z, s=1)

Compute the first s+1 exponents of z.

Returns :[I,z,z**2,...,z**s]

class odelab.phi_pade.RationalFraction(numerator, denominator)

Super simple Rational Function class, used to compute N(Z)/D(Z) for polynomials N and D.