rcorless

530 Reputation

10 Badges

4 years, 273 days

Social Networks and Content at Maplesoft.com

Editor-in-Chief of Maple Transactions (www.mapletransactions.org), longtime Maple user (1st use 1981, before Maple was even released). Most obscure piece of the library that I wrote? Probably `convert/MatrixPolynomialObject` which is called by LinearAlgebra[CompanionMatrix] to compute linearizations of matrix polynomials in several different bases. Do not look at the code. Seriously. Do not look. You have been warned.

MaplePrimes Activity


These are replies submitted by rcorless

@Joe Riel thanks for "thisproc".  I did use "procname" before but I was genuinely mistaken as to its meaning.  The keyword "thisproc" is correct.

However.  

g := proc(n) option remember; if n < 2 then 1; else procname(n - 1) + procname(n - 2); end if; end proc;

G := Compiler:-Compile(g);

CodeTools:-Usage( G(40) ); # Takes 1.36 seconds on my machine (!)

CodeTools:-Usage(g(40)); # Records 0ns .  Um, what?

Now, the first time G is invoked it might need something.  So if one progressively calls G(10), G(20), G(30), one sees the time creep up.  Then G(40) takes 1.3 seconds.

Or, don't bother.  It still takes about the same time, vastly longer than the uncompiled version.

Perhaps the lesson is, only compile more expensive programs?  There's some kind of overhead involved here?  Except that the larger the input, the longer it takes.  I suspect that the problem is "option remember" here, which is not being respected by the compiler.

This seems to be correct.  The "option remember" is being ignored, and the cost becomes exponential in n (as is well-known).  A non-recursive g works much better:

g := proc(n::nonnegint) local k, g0, g1, g; g0 := 1; g1 := 1; for k from 2 to n do g := g0 + g1; g0 := g1; g1 := g; end do; return g; end proc;

G := Compiler:-Compile(g);

Now G(91) gives the correct answer (no rounding error! Which seems remarkable because floats can't be involved, therefore, only integers) but G(92) does not give a negative answer but rather raises the overflow exception. 

Also the time taken is not detectable by CodeTools:-Usage; too fast.

Ok, thanks; I learned something.

-r

Aha! I used procname later, not thisproc. I shall investigate the differences, if any.

@Kitonum add( (-1)^(n+1)/(2*n-1), n=1..4) is *much* more intelligible than your (admittedly correct, just obscure) solution :)

"sum" does symbolic summation

"add" is (like seq) a command with special evaluation rules but is meant for adding together a fixed number of terms.

 

This is not an Ordinary differential equation, but a Functional differential equation instead.  Other examples include simple delay differential equations like y'(t) = a*y(t-1) and then you have to specify the "history" on an interval (say (0 <= x <= 1).  It gets complicated. More fun functional equations include things like y'(t) = a*y(t/2) and this one can be started from t=0 (doesn't need a history).  Yours is of a type I have never seen before.  The derivative at x=0 depends on the behaviour of the function at x=infinity!  If you specified the function on 0 < x <= 1 then it looks like it might be solvable on x>1, but I am not so sure. 

I'll think about this a little bit now---might try a Laplace transform, for instance---but yes as has already been observed, f(x) = 0 identically is a solution to this equation (but I have no idea if there are any others).

@Carl Love the word "height" was more commonly used for polynomials; there is a paper by Mahler for instance, but maybe the usage goes back to Littlewood.  The "height" of a polynomial is the infinity norm of the vector of coefficients (usually in the monomial basis); its "length" is the 1-norm.  In other words, the height of a polynomial was the largest absolute value of any coefficient.

Similarly, the height of a matrix is the largest absolute value of any entry.  If the entries are all integers, then there is a minimal possible height as well (and this matters: of course one can scale any matrix so its height is 1; but only at the cost of possibly introducing small entries).  The zero matrix has height 0 (and it's the only one).  A matrix like the Mandelbrot matrix of https://doi.org/10.5206/mt.v1i1.14037 has height 1, and this is pretty obviously the smallest interesting height.  It's astounding that the characteristic polynomial has height (we say "characteristic height", to distinguish from the height of the matrix) that is exponential in the degree, and is thus vastly ill-conditioned for evaluation and rootfinding.

 

The name "Bohemian" comes from the acronym "BOunded HEight Matrix of Integers" (BOHEMI) which is close enough.

 

intsolve was written in 1991, and fracdiff in 2004.  I'm kind of pleased that MapleMathMatt managed to force intsolve to work with fracdiff.  Nice work!

@Carl Love I once proposed that subject line as the title of our article on Gamma; both Jon and Judi recoiled in horror, even though it's a palindrome.

I had not known that this integral was expressible in terms of Psi.  Thank you---that derivation is cool, too; I like FormalPowerSeries.

@Kitonum I hadn't known about the "continuous" option; I am going to go look at it now.  This is something I really really should have known.

Okay, then---TIL about Algebra Tiles!  Totally going to show this to my daughter (who is a teacher).

Blink, blink!

@jeffreyrdavis75 I like the pictures too.  And, yes, one can do 3d versions (we haven't, so far, but there is a good reason to if we want to colour by condition number and by density both at once).  Have at it!

 

I do not know why my animation tool bar does not appear.  However, I can run the animation by using right-click and the context menu.  Fine.  Principle of conservation of recusancy, I guess.

 

Now I have a different problem; plots[animate] has replaced the old "insequence=true" option, and I find my handling of it to be hit-and-miss.  The following produces a static plot only, not an animation, the way I expected it to.  Any thoughts?

 

plots[animate](plots[pointplot], [[seq([n, sin(combinat[fibonacci](k + 1)*n/combinat[fibonacci](k))], n = 1000 .. 2000)], symbol = point, axes = boxed, labels = [n, typeset(sin(combinat[fibonacci](k + 1)*n/combinat[fibonacci](k)))], labeldirections = [horizontal, vertical]], k = [seq(j, j = 1 .. 12)]);

 

I'll try "with(plots)" now just in case that works, but I would expect this to work with the long forms.

@acer This earns you an acknowledgement in the paper :)  I am happy to use your real name there if you like, but if you want to keep the pseudonym I'm happy to do that, too!

@mmcdara I hope to get to a planned improvement in that package (our original version was a little better for derivatives).  Let's see what happens, now that I am retired :)

Here's one paper using it---there are others listed on David Jeffrey's home page.

https://www.uwo.ca/apmaths/faculty/jeffrey/pdfs/large.pdf

Thank you very much!  That explains everything.  Yes, circuits is what I meant (and not knowing any graph theory, I was unaware of the different definitions of "cycle" and "circuit").  When I converted to an undirected graph (to make CycleBasis work) via the UnderlyingGraph command, the pair of directed edges [12,4] and [4,12] collapsed to a single undirected edge; at that point no cycle remained.  But if three or more vertices are involved in a circuit (in my application), they have to collapse to a cycle.  Ok, that explains that puzzle very well.

In my application the directed graph comes from a deterministic dynamical system; there is only one edge leading from any given vertex, and so detection of a cycle in the underlying graph is automatically detection of a circuit in the original; but thank you for teaching me how to check it in a more general situation.  Your showing how to use ListTools there was also helpful---didn't know about Rotate!

1 2 3 Page 2 of 3