Featured Post

We’re kicking off 2018 right, with another Meet Your Developers interview! This edition comes from Erik Postma, Manager of the Mathematical Software Group.

To catch up on previous interviews, search the “meet-your-developers” tag.

Without further ado…


  1. What do you do at Maplesoft?
    I’m the manager of the mathematical software group, a team of 7 mathematicians and computer scientists working on the mathematical algorithms in Maple (including myself). So my work comes in two flavours: I do the typical managerial things, involving meetings to plan new features and solve my team’s day to day problems, and in the remaining time I do my own development work.
  2. What did you study in school?
    I studied at Eindhoven University of Technology in the Netherlands. The first year, I took a combined program of mathematics and computer science; then for the rest of my undergrad, I studied mathematics. The program was called Applied Mathematics, but with the specialization I took it really wasn’t all that applied at all. Afterwards I continued in the PhD program at the same university, where my thesis was on a subject in abstract algebra (Lie algebras over finite fields).
  3. What area(s) of Maple are you currently focusing on in your development?
    I’ve spent quite a bit of time over the past two years making the facilities for working with units of measurement in Maple easier to use. There is a very powerful package for doing this that has been part of Maple for many years, but we keep hearing from our users it’s difficult to use. So I’ve worked on keeping the power of the package but making it easier to use.
  4. What’s the coolest feature of Maple that you’ve had a hand in developing?
    This was actually working on a problem in a part of the code that existed long before I started with Maplesoft. We have a very clever algorithm for drawing random numbers according to a custom, user-specified probability distribution. I wrote about it on MaplePrimes in a series of four blog posts, here. I’ve talked at various workshops and the like about this algorithm and how it is implemented in Maple.
  5. What do you like most about working at Maplesoft? How long have you worked here?
    I love working at the crossroads of mathematics and computer science; there aren’t many places in the world where you can do that as much as at Maplesoft. But the best thing is the people I work with: us mathematicians are all crazy in slightly different ways, and that makes for a very interesting working environment.
  6. Favourite hobby?
    Ultimate frisbee. I captain a mixed (i.e., coed) team called The Clockwork. (We play in orange jerseys – it references the book/movie A Clockwork Orange.) We play in a couple of local leagues, and some of the other members also work here. We don’t win much – but we work hard and have fun!
  7. What do you like on your pizza?
    Mushrooms. Mushrooms on everything!
  8. What’s your favourite movie?
    Probably Black Book, a dark movie about the Dutch resistance in the second world war from 2006, directed by Paul Verhoeven. I think what I like best about it is that it highlights the moral shades of grey in even so morally elevated a group as the resistance.
  9. What skill would you love to learn? Why?
    I’d love to learn to speak Russian! I’m trying, but I have a very hard time with it. It would allow me to communicate with my in-laws more easily; they speak Russian.
  10. Who’s your favourite mathematician?
    Oh, so many to choose from! I’m torn between:
  • Ada Lovelace (1815-1852), known as the first programmer.
  • Felix Klein (1849-1925), driving force behind a lot of research into geometries and their underlying symmetry groups.
  • Wilhelm Killing (1847-1923), a secondary school teacher who made big contributions to the theory of Lie algebras.

Or wait, can I choose my wife?

Featured Post

One case where an "expansion beyond all orders" may be needed is investigating the asymptotic behavior of the difference of two functions with coinciding dominant series.

We are interested in the asymptotic behavior of F(z) for large positive z:

h1 := proc (z) options operator, arrow; hypergeom([1/2], [5/4, 3/2, 7/4], z) end proc; h2 := proc (z) options operator, arrow; (3/4)*sqrt(2*Pi)*hypergeom([1/4], [3/4, 5/4, 3/2], z)/(GAMMA(3/4)*z^(1/4)) end proc; F := proc (z) options operator, arrow; h1(z)-h2(z)+(3/8)*sqrt(Pi)/sqrt(z) end proc

series does not succeed:

series(F(z), z = infinity, 20)



The reason is that the dominant terms containing the factor exp(3*z^(1/3)) in the two hypergeometric functions cancel exactly, and we have to look for the subdominant terms.

The order of the leading terms can be found from DETools:-formal_sol:

deq1 := FunctionAdvisor(DE, h1(z), f(z))[2, 1]

diff(diff(diff(diff(f(z), z), z), z), z) = -(15/2)*(diff(diff(diff(f(z), z), z), z))/z-(195/16)*(diff(diff(f(z), z), z))/z^2+(1/32)*(32*z-105)*(diff(f(z), z))/z^3+(1/2)*f(z)/z^3


DETools:-formal_sol(deq1, f(z), z = infinity, order = 0)

[(1/z)^(1/2), -exp(-3/(-1/z)^(1/3))/z, -exp(3*(-1)^(1/3)/(-1/z)^(1/3))/z, -exp(-3*(-1)^(2/3)/(-1/z)^(1/3))/z]


As expected, one of the solutions (the third one for positive z) contains the exp(3*z^(1/3)) factor, the leading term being of the order exp(3*z^(1/3))/z.

Another, subdominant, solution is algebraic and, in fact, is a series containing only one term, as 1/z^(1/2) is an exact solution. It will turn out that the algebraic part in F(z) also cancels out.

Thus we have to look for the subsubdominant terms, which contain decaying exponentials. We will accomplish this by applying the steepest descent method to the integral representations of h1(z) and h2(z).

ms := convert([h1(z), h2(z)], MeijerG)

[(3/32)*Pi*2^(1/2)*MeijerG([[1/2], []], [[0], [-1/4, -1/2, -3/4]], -z), (3/32)*2^(1/2)*Pi*MeijerG([[3/4], []], [[0], [1/4, -1/4, -1/2]], -z)/z^(1/4)]


m2g := proc (m, y) local a, b, c, d; a, b := op(op(1, m)); c, d := op(op(2, m)); -((1/2)*I)*mul(`~`[GAMMA](`~`[`-`](1+y, a)))*mul(`~`[GAMMA](`~`[`-`](c, y)))*op(3, m)^y/(Pi*mul(`~`[GAMMA](`~`[`-`](b, y)))*mul(`~`[GAMMA](`~`[`-`](1+y, d)))) end proc

gs := applyrule(conditional(e::anything, _op(0, e) = MeijerG) = 'm2g(e, y)', ms)

[-((3/64)*I)*2^(1/2)*GAMMA(1/2+y)*GAMMA(-y)*(-z)^y/(GAMMA(5/4+y)*GAMMA(3/2+y)*GAMMA(7/4+y)), -((3/64)*I)*2^(1/2)*GAMMA(1/4+y)*GAMMA(-y)*(-z)^y/(z^(1/4)*GAMMA(3/4+y)*GAMMA(5/4+y)*GAMMA(3/2+y))]


gs[2] := combine(eval(gs[2], [1/z^(1/4) = exp(I*Pi*(1/4))/(-z)^(1/4), y = y+1/4]), power)



h1(z) and h2(z)are the integrals of gs[1] and of gs[2] over the same path, which is a loop encircling the poles ofGAMMA(-y) and of GAMMA(-1/4-y). Now a standard technique is to extend the integration contour far to the left, while still keeping both endpoints at "+infinity". Then the arguments of the gamma functions can be made large everywhere on the integration path, and the gamma functions can be replaced by their asymptotic approximations.

When moving the contour, we have to take into account the pole of the integrand at y = -1/2. The other poles of GAMMA(1/2+y) will be cancelled by the zeros of 1/GAMMA(3/2+y), which is why the algebraic part of the expansion will contain the single term of the order 1/z^(1/2).

This is the negative of the third term in F(z):

`assuming`([simplify((2*Pi*I)*residue(gs[1]-gs[2], y = -1/2))], [z > 0])



Expanding the gamma functions produces terms containing exp(-I*Pi*y) and exp(I*Pi*y)

`assuming`([simplify(convert(MultiSeries:-series((gs[1]-gs[2])/z^y, y = -infinity, 1), polynom))], [z > 0]); collect(convert(%, exp), exp)



As we shall see, those terms have saddle points y0(z) = exp(`&+-`((1/3)*(2*Pi*I)))*z^(1/3) located in the left half-plane and contribute exponentially small factors exp(3*y0(z)). The terms for which the saddle point would be located at y = z^(1/3) have cancelled out, thus cancelling the exponentially large contributions. Another possible way to achieve the same result was to write h1(z)-h2(z) as a single Meijer G-function -(3/32)*MeijerG([[1/2], []], [[-1/4, 0], [-3/4, -1/2]], z).

We write the first term above in the form g(y)*exp(f(y)):

f := proc (z, y) options operator, arrow; -3*y*(ln(-y)-1)-I*Pi*y+y*ln(z) end proc

g := proc (y) options operator, arrow; (-3/64-(3/64)*I)*sqrt(-1/y)/(sqrt(Pi)*y^3) end proc

diff(f(z, y), y)



For this to become zero, we need argument(-y) = -(1/3)*Pi, and thus y = exp((1/3)*(2*I)*Pi)*z^(1/3). We can visualize the paths where the imaginary part of f(z, y) stays constant. The path of the steepest descent is the one that goes through the saddle point in the direction exp(I*Pi*(1/3)); the blue color indicates smaller values of the real part of f(z, y):

y0 := proc (z) options operator, arrow; exp(((2/3)*I)*Pi)*z^(1/3) end proc

(proc () local z; z := 2; plots:-display(plots:-contourplot(Re(f(z, u+I*v)), u = -5 .. 5, v = -5 .. 5, contours = ([seq])(Re(f(z, y0(z)))+i, i = -30 .. 6, 6), filledregions, coloring = [blue, red], grid = [100, 100]), plots:-implicitplot(Im(f(z, u+I*v)-f(z, y0(z))), u = -5 .. 5, v = -5 .. 5, gridrefine = 5, color = green), plot([cos((1/3)*Pi)*xi+Re(y0(z)), sin((1/3)*Pi)*xi+Im(y0(z)), xi = -3 .. 3], linestyle = dot, color = white), axes = boxed) end proc)()


The real part of f(z, y) has a maximum along this path at y0(z).

`assuming`([(`@`(`@`(simplify, evalc), series))(f(z, y0(z)+exp(I*Pi*(1/3))*xi), xi = 0, 3)], [z > 0]); quad := convert(%, polynom)



Now we can compute the lead asymptotic term contributed by the saddle point y0(z):

lt1 := `assuming`([(`@`(simplify, evalc))(g(y0(z))*exp(I*Pi*(1/3))*(int(exp(quad), xi = -infinity .. infinity)))], [z > 0])



We repeat the same procedure for the second term of the integrand.

f := proc (z, y) options operator, arrow; -3*y*(ln(-y)-1)+I*Pi*y+y*ln(z) end proc

g := proc (y) options operator, arrow; (3/64-(3/64)*I)*sqrt(-1/y)/(sqrt(Pi)*y^3) end proc

diff(f(z, y), y)



y0 := proc (z) options operator, arrow; exp(-((2/3)*I)*Pi)*z^(1/3) end proc

The direction should be chosen as exp((1/3)*(2*I)*Pi) to be consistent with the direction of the integration contour, which goes from the lower to the upper half-plane.

lterm := proc (gy, fy, eq, dir) options operator, arrow; (eval(gy*exp(fy), eq))*dir*sqrt(-2*Pi/((eval(diff(fy, `$`(y, 2)), eq))*dir^2)) end proc

lt2 := `assuming`([(`@`(simplify, evalc))(lterm(g(y), f(z, y), y = y0(z), exp((1/3)*(2*I)*Pi)))], [z > 0])



Combining the two results yields the leading term of F(z). The next terms can be obtained by expanding gs[1] and gs[2] to higher orders.

Fasympt := unapply(simplify(lt1+lt2), z)

proc (z) options operator, arrow; (1/32)*exp(-(3/2)*z^(1/3))*3^(1/2)*2^(1/2)*(cos((3/2)*z^(1/3)*3^(1/2))+sin((3/2)*z^(1/3)*3^(1/2)))/z end proc


(proc () Digits := 50; plot(`~`[`*`](exp((3/2)*z^(1/3)), [F(z), Fasympt(z)]), z = 1000 .. 10000, linestyle = [solid, dot], thickness = [1, 5], axes = frame) end proc)()


Download steep.mw

How do I rearrange an expression?

Maple asked by pluto 10 Yesterday

animate a set of vectors

Maple asked by EngM 10 Today

long time calculation

Maple asked by arianbig 10 February 15

maple 6 download

Maple 6 asked by Lali_miani... 15 February 14

convert command

Maple 12 asked by tsunamiBTP... 253 February 14