Computing with rational functions
Tags: RKFUN, rational functions Name: Mario Berljafa and Stefan Güttel File: example_rkfun.m (PDF version) Date: 27/05/2015
Contents
Introduction
This toolbox comes with an implementation of a class called rkfun, which is the fundamental data type to represent and work with rational functions. Objects of this class are produced by the rkfit function (described in [2,3]), which in its simplest use case attemts to find a rational function such that
where are square matrices and
is a vector of compatible sizes. For example, let us consider
,
, and
. We know that
is a minimizer for the above problem. Let us try to find it via RKFIT:
N = 100; A = gallery('tridiag',N); I = speye(N); F = (A - 2*I)^2*inv(A^2 - I)*inv(A - 4*I); b = eye(N,1); xi = inf*ones(1,5); % initial poles at infinity maxit = 5; tol = 1e-12; [xi,ratfun,misfit] = rkfit(F,A,b,xi,maxit,tol,'real');
As we started with initial poles (at infinity), RKFIT will search for a rational function
of type
. When more than 1 iteration is performed and the tolerance tol is chosen sufficiently large, RKFIT will try to reduce the type of the rational function while still maintaining a relative misfit below the tolerance, i.e.,
. Indeed, a reduction has taken place and the type has been reduced to
, as can be seen from the following output:
ratfun
ratfun = RKFUN object of type (2,3). Real Hessenberg pencil (H,K) of size 4-by-3. coeffs = [ -0.0869, -0.394, -1.33, -1.5 ]
We can now perform various operations on the ratfun object, all implemented as methods of the class rkfun. To see a list of all methods just type help rkfun. We will now discuss some methods in more detail.
Evaluating an rkfun
We can easily evaluate at any point (or many points) in the complex plane. The following command will evaluate
and
simultaneously:
format long
ratfun([2 ; 3+1i])
ans = 0.000000000000034 + 0.000000000000000i 0.011764705882320 - 0.152941176470882i
We can also evaluate as a matrix function, i.e., computing
for matrices
and
, by using two input arguments. For example, by setting
we effectively compute the full matrix function
:
M = [ 3 , 1 ; 0 , 3 ]; B = eye(2); R = ratfun(M,B)
R = -0.125000000000272 -0.281250000000544 0 -0.125000000000272
Plotting
As can be evaluated at any point in the complex plane, it is straightforward to produce plots of this function. For example, here is a contour plot of
over the complex region
:
figure(1)
[X,Y] = meshgrid(-2:.01:5,-1:.01:1); Z = X + 1i*Y;
R = ratfun(Z);
contourf(X,Y,log10(abs(R)),-4:.5:2)
colormap hsv, colorbar

Another command ezplot can be used to get a quick idea of how ratfun looks over an interval on the real axis, in this case, :
figure(2)
ezplot(ratfun,[-2,5])
ylim([-5,5])
grid on

Pole- and root-finding
From the above plot we guess that has poles at
and
and a root at
, which is to be expected from the definition of
. The two commands poles and roots do exactly what their names suggest:
pls = poles(ratfun) rts = roots(ratfun)
pls = 0.999999999999999 -1.000000000000605 4.000000000000003 rts = 1.999999547063492 2.000000452935510
As expected from s type (2,3) rational function, there are two roots and three poles. Note that the pole at is identified with slightly less accuracy than the poles at
and
. This is because the point
is outside the spectral interval of
and hence sampled less accurately. Also the double root at
is identified up to an accuracy of
only. This is not surprising as the function is flat nearby multiple roots. However, the backward error of the roots is small:
ratfun(rts)
ans = 1.0e-15 * 0.138777878078145 0.138777878078145
Basic arithmetic operations and differentiation
We have implemented some very basic operations on rkfun's, namely, the multiplication by a scalar and the addition of a scalar. The result of such operations is again an rkfun object. For example, the following command computes points where
:
z = roots(2*ratfun - pi) ratfun(z)
z = -0.405042315586584 0.857694376021584 ans = 1.570796326794896 1.570796326794897
We have currently not implemented the addition and multiplication of two rkfun's, though this is doable in principle. However, we can already differentiate a rational function using the diff command. The following will find all local extrema of by computing the roots of
:
extrema = roots(diff(ratfun))
extrema = 0.407592768316864 + 0.000000000000000i 2.000000007213609 + 0.000000000000000i 2.796203600878076 + 2.627131616351852i 2.796203600878076 - 2.627131616351852i
There are two real extrema which we can add to the above plot of :
figure(2), hold on plot(extrema(1:2),ratfun(extrema(1:2)),'ro')

The syntax and "feel" of these computations is inspired by the Chebfun system [4], which represents polynomials via Chebyshev interpolants and allows for many more operations to be performed than our rkfun implementation. Here we are representing rational functions through their coefficients in a discrete-orthogonal rational function basis. Working with rational functions poses some challenges not encountered with polynomials. For example, the indefinite integral of a rational function is not necessarily a rational function but may contain logarithmic terms.
Multiple precision computations
Objects of class rkfun can be converted to MATLAB's Variable Precision Arithmetic (VPA) as follows:
vpa(ratfun)
ans = RKFUN object of type (2,3). Real Hessenberg pencil (H,K) of size 4-by-3. Variable precision arithmetic (VPA) activated. coeffs = [ -0.0869, -0.394, -1.33, -1.5 ]
Alternatively, we can also use the Advanpix Multiple Precision (MP) toolbox [1], which is typically more efficient and reliable than VPA. We recommend the use of this toolbox in particular for high-precision root-finding of rkfun's and conversion to partial fraction form, as MATLAB's VPA does not support generalized eigenvalue problems:
ratfun = mp(ratfun)
ratfun = RKFUN object of type (2,3). Real Hessenberg pencil (H,K) of size 4-by-3. Multiple precision arithmetic (ADVANPIX) activated. coeffs = [ -0.0869, -0.394, -1.33, -1.5 ]
When evaluating a multiple precision ratfun, the result will be returned in multiple precision:
ratfun([2 ; 3+1i])
ans = Columns 1 through 1 3.434024222017004731285023676831756e-14 + 0i 0.01176470588232020213010063005597186 - 0.1529411764708824495013595783793991i
It is important to understand that, although the evaluation of ratfun is now done in multiple precision, the computation of ratfun using the rkfit command has been performed in standard double precision. rkfit does not currently support the computation of rkfun's in multiple precision. Here are the roots of ratfun computed in multiple precision:
roots(mp(ratfun))
ans = 1.999999546081609784329105502644233 2.000000453917395201475935899049275
Conversion to partial fraction form
It is often convenient to convert a rational function into its partial fraction form
The residue command of our toolbox performs such a conversion. Currently, this only works when the poles of
are distinct and
is not of superdiagonal type (i.e., there is no linear term in
). As the conversion to partial fraction form can be an ill-conditioned transformation, we recommend to use residue in conjunction with the multiple precision feature. Here are the poles
and residues
(
), as well as the absolute term
, of the function
defined above:
[alpha,xi,alpha0,cnd] = residue(mp(ratfun)); double([xi , alpha]) double(alpha0)
ans = -1.000000000000605 0.900000000001050 0.999999999999999 -0.166666666666673 4.000000000000003 0.266666666667165 ans = 1.025954395909300e-16
The absolute term is close to zero as is of subdiagonal type. The output cnd of residue corresponds to the condition number of the transformation to partial fraction form. In this case of a low-order rational function with well separated poles the condition number is actually quite moderate:
cnd
cnd = 56.150570384550242
Conversion to quotient and continued fraction form
Our toolbox also implements the conversion of a rkfun to quotient form with two polynomials
and
given in the monomial basis. As with the conversion to partial fraction form, we recommend performing this transformation in multiple precision arithmetic due to potential ill-conditioning. Here we convert
into the
form and evaluate it at
using MATLAB's polyval command:
[p,q] = poly(mp(ratfun)); polyval(p,2)./polyval(q,2)
ans = 3.434024222017004236112745256053473e-14
A rkfun can also be converted into continued fraction form
1 r(z) = -------------------------------------------------------- + absterm hh(1)*z + 1 ---------------------------------------------- h(1) + 1 --------------------------------------- hh(2)*z +...+ 1 ------------------------- h(m-1) + 1 ---------------- hh(m)*z + 1/h(m)
as follows:
[h,hh,absterm,cnd] = contfrac(mp(ratfun));
That's it for this tutorial. Note that more methods will be added over time and we'd be happy to receive any feedback or bug reports. For more details about the internal representation of rkfun, see [3].
The following command merely creates a nice thumbnail.
figure(1), hold on, plot(NaN) % create thumbnail

References
[1] Advanpix LLC., Multiprecision Computing Toolbox for MATLAB, ver 3.8.3.8882, Tokyo, Japan, 2015. http://www.advanpix.com/.
[2] M. Berljafa and S. Güttel. Generalized rational Krylov decompositions with an application to rational approximation, accepted for publication in SIAM J. Matrix Anal. Appl., 2015. MIMS EPrint 2014.59 (http://eprints.ma.man.ac.uk/2278/), Manchester Institute for Mathematical Sciences, The University of Manchester, UK, 2014.
[3] M. Berljafa and S. Güttel. The RKFIT algorithm for nonlinear rational approximation, MIMS EPrint 2015.38 (http://eprints.ma.man.ac.uk/2309/), Manchester Institute for Mathematical Sciences, The University of Manchester, UK, 2015.
[4] T. A. Driscoll, N. Hale, and L. N. Trefethen. Chebfun Guide, Pafnuty Publications, Oxford, 2014. http://www.chebfun.org