Related to issues #73 and #74. Would be useful to solve all of these at the same time.

I don't have access to a copy of Levine at the moment, will see if I
can get hold of one.

If I understand your interpretation correctly, the idea is to split
the delay up, and find Pade approximations to each and cascade them.

An alternative might be:

1. Find high order Pade approximation for exp(-s) (i.e., T=1)
2. Split this up into second-order-sections, as is done for high-order digital filters
3. Scale each second order section by T; one only then ever scales by T and T2, not Tn

If step 1 and 2 are done using high precision math
(fractions.Fraction, or probably better gmpy2.mpq), I think one might
be able get to absurdly high order, over a wide range of T.

Having said all of this, I'm curious to know what the practical or
academic applications of high order (n > 5?) Pade approximations are?

If I understand your interpretation correctly, the idea is to split
the delay up, and find Pade approximations to each and cascade them.

Well, that's what Levine says (I don't have access to it either, just via Google books).

But in the last few weeks I had some ideas along the same lines as what you mentioned:

An alternative might be:

1. Find high order Pade approximation for exp(-s) (i.e., T=1)
2. Split this up into second-order-sections, as is done for high-order digital filters
3. Scale each second order section by T; one only then ever scales by T and T2, not Tn

and I've got some really interesting results, however it may take a week or two for me to write up as I work full-time and have some life issues coming up very shortly.

ok, this took more than a week or two to write up, but here it is.

http://www.embeddedrelated.com/showarticle/927.php

The key is getting accurate values for poles and zeros without going through polynomial coefficients. Applications of high-order Pade approximations are time-domain simulation; if there's a good alternative that works well to simulate time delays in continuous-time systems, I'd like to know what they are, as I'm not aware of any, except that Simulink has some proprietary method of doing so.

Just skimmed your write-up, looks great. Looking forward to giving it a thorough read. (Had previously read your Chebyshev polynomial article, which was great.)

Regarding simulating delays: have you looked at delay differential equation solvers? There do appear to be some free/open source ones, e.g., http://www.unige.ch/~hairer/software.html (have never used them myself).

Sorry for jumping in suddenly: this is something I have been chewing in my mind for a while too.

As far as I know, Simulink doesn't pull off any fancy tricks other than the discretization via Thiran filters for the systems and numerical integration for the delay simulation. You can look it up via doc thiran for the Spectrum article reference for the discretization part which is basically Lagrange approximation instead of your Chebyschev argument.

I'm (very slowly) writing something similar to that kind of simulation substitute for lsim(2) but now I'm curious about that method you have referred to. Would you please direct me to that proprietary method reference?

By the way, your article was really fun!

Matlab has dde23 [1], which doesn't obviously use Thiran filters (interpolation is mentioned in the docs).

[2] says Simulink can handle variable delays; do Thiran filters allow for that?

[1] http://www.mathworks.com/help/matlab/ref/dde23.html
[2] https://books.google.co.za/books?id=6d7iAAAAQBAJ&lpg=PT183&ots=w4MLvRyJBO&dq=simulink%20dde23&pg=PT183#v=onepage&q=simulink%20dde23&f=false

dde23 doesn't solve stiff systems as far as I know. Also if you use ode23/45 it still doesn't have the problems given in the blogpost of @jason-s. So I don't see how it comes into play.

For time variable delays, it uses intersample integration but I'm not sure how they precondition the right hand side for it in the simulink.

Did you mean dde23 as the proprietary method above?