1484 Reputation

19 Badges

15 years, 180 days

Social Networks and Content at

Maple Application Center
I am the manager of the Mathematical Software Group, working mostly on the Maple library. I have been working at Maplesoft since 2007, mostly on the Statistics and Units packages and on contract work with industry users. My background is in abstract algebra, in which I completed a PhD at Eindhoven University of Technology. During my studies I was always searching for interesting questions at the crossroads of math and computer science. When I came to Canada in 2007, I had nothing but a work permit, some contacts at Maplesoft (whom I had met at a computer algebra conference a year earlier), and a plan to travel around beautiful Canada for a few months. Since Maplesoft is a company that solves more interesting math and computer science related questions than just about anywhere else, I was quite eager to join them, and after about three months, I could start.

MaplePrimes Activity

These are answers submitted by epostma

Hi Neill,

You've already found that you can do rhs(Eq1) = lhs(Eq1). Just a tiny bit shorter is this:

algsubs((rhs = lhs)(Eq1), Eq2);

This works because if you use something like f = g as a function applied to x, then this evaluates to f(x) = g(x).

As you already suggest in the title, you can also use algsubs. Whenever the thing you are replacing is not just a variable name, algsubs is more reliable than subs. So I'd recommend using that for this case. An alternative would be to roll your own version of algsubs using solve or isolate.

Hope this helps,

Erik Postma

Hi afeddersen,

Sorry - I seem to have a habit of answering direct questions with "here is something different you could do that might accomplish your goals in a better manner", and I'm afraid this is going to be such an answer.

The label feature in Maple is great, but it mostly works with direct user interaction: there's no easy access to programmatically retrieve the value of label (x.y) for example, if x and y are integers; instead you have to enter these numbers manually. Loops, of course, are very much programmatic constructs, and thus these two don't mesh very well.

One option would be to just store whatever value is printed out in a ?table or an ?rtable (such as a ?Vector). So if you have a loop that looks like this:

running_total := 0;
for i to 25 do
  running_total := running_total + i:
  i^2 - running_total;
end do;

to print the numbers i*(i-1)/2, you could instead decide to store these in a Vector:

results := Vector(25);
running_total := 0;
for i to 25 do
  running_total := running_total + i:
  results[i] := i^2 - running_total;
end do;

Now you can access the ith value as results[i] - with the added advantage over labels that it also works programmatically, not just when you type it directly. Does that help for what you're trying to do?

Erik Postma

Hi sund0002,

Glad that the plot worked. For import of text files as matrices, you can use this:

Data := ImportMatrix("/home/epostma/downloads/200-1.1.txt", skiplines=2, datatype=float[8]);
Data := Data[..-2];

The skiplines=2 option is to get rid of the first two lines of the file. Because there's an empty line at the end of the file, you'll need to get rid of the last row of the matrix as well, which happens in the second command.

The datatype=float[8] option is not necessary, but it makes subsequent operations on the matrix faster. This is especially useful if you have larger data sets.

By the way, note that you now have a Matrix, not an Array. This should make the last step of our previous conversation unnecessary (but it won't hurt either): the conversion from Array to Vector.

Hope this helps,

Erik Postma

Hi Morten,

The easiest way of doing what you ask is to run map(assign, A, 0). However, whenever variables have a mathematical meaning (as they probably do if they occur in a matrix), there is often a better solution than just assigning a value to them. For example, if you want to know the value of an expression where the variables have that value, you can use a set of equations that provide values for the variables:
:= map(`=`, indets(A, name), 0);
and then evaluate the expression at that set of values:
    result := eval((a + b + 2)^4, values);
This makes it much easier to also evaluate the expression at other values of your variables: if you had values assigned to a and b, then result := (a + b + 2)^4 just assigns the value 16 to result, with no way of retrieving the original expression in terms of a and b.

Hope this helps,

Erik Postma

The best solution for this sort of thing is generally ?plots/surfdata , except that requires that you already have data for a uniform grid of points, and your data is not in a uniform grid. So you'll need to fit your data to a surface first. And frankly, the arrangement of data that you show here is not going to be enough to get a very useful fit in between the data points that are given. But maybe you have more data that would just be too voluminous to post? Let's have a look at the options at any rate.

The usual suspects for constructing such a fit are to be found in the ?CurveFitting package, but here most procedures either again require data for a uniform grid of points, or they accept only one-dimensional input data whereas you have 2D input. Another option would be the ?Statistics/Regression subpackage. That provides two functions that could be useful here: ?Statistics/LinearFit and ?Statistics/NonlinearFit. Both can fit a nonlinear expression to data; the first command requires that that expression is linear in the parameters that are to be fitted but provides globally optimal results, the second command is more widely applicable but provides only locally optimal results. The ?Statistics/Fit command provides a way to let Maple choose which of the two commands is appropriate.

A final option, one that could be considered nonparametric, is to use a kernel density fit. Again, Maple has a facility for that (?Statistics/KernelDensity), but it only allows for 1D input / 1D output data. But this would maybe be useful to construct ourselves. Essentially what want to do then, is: given x- and y-coordinates, take a weighted average of the c-values corresponding to our input points, where the weight for ci depends on the vector (x - ai, y - bi). We'll do that as follows:

a := <.19375, .19494, .19581, .19632, .19646, .19624, .19565, .19472, .19349, .19198, .19024,  .18833>;
b := <-0.1517e-1, -0.1515e-1, -0.1512e-1, -0.1509e-1, -0.1506e-1, -0.1503e-1, -0.15e-1, -0.1497e-1, -0.1495e-1, -0.1492e-1, -0.1491e-1, -0.149e-1>;
c := <0.112346e-10, 0.111963e-10, 0.111544e-10, 0.111101e-10, 0.110649e-10, 0.110201e-10, 0.10977e-10, 0.10937e-10, 0.109011e-10, 0.108705e-10, 0.108458e-10, 0.108277e-10>;

kdFit := proc(xdata, ydata, zdata)
local xradius, yradius;
  xradius := (max - min)(xdata);
  yradius := (max - min)(ydata);
  return proc(x, y)
  local d, values;
    values := zip((xx, yy) -> exp(- 10 * (xx^2 + yy^2)), (xdata -~ x) / xradius, (ydata -~ y) / yradius);
    return (values . zdata) / add(d, d in values);
  end proc;
end proc;
ourFit := kdFit(a, b, c);
p1 := plots[contourplot3d](ourFit, min(a) .. max(a), min(b) .. max(b), axes=boxed, style=surfacecontour);
p2 := plots[pointplot3d](a, b, c, color=black, symbolsize=25);
plots[display](p1, p2);

Instead of the 3D contour plot that I plotted here, you can also create a regular, 'flat', contour plot by using plots[contourplot] - in which case you're probably less interested in the points, so you just need p1.

Hope this helps,

Erik Postma

Edit: had a typo in the code - fixed now.

It's not quite what you're asking for since it's also runtime behaviour, but I'd say it is somewhat in the spirit of what you're asking since it might allow you to catch a mistake before it actually leads to a wrong result. There are some underused safety-checking features that you can turn on with the command kernelopts(assertlevel = 2) (or by starting Maple with the -A 2 command line option): type assertions on local variablesassignments, and procedure return values. This makes Maple check (at runtime) the types of certain values, if you provide the type that they need to have. For example:

p := proc();
local a :: posint;
  a := -3/2;
  return a;
end proc;

q := proc()
local a;
  a :: posint := -3/2;
  return a;
end proc;

r := proc() :: posint;
local a;
  a := -3/2;
  return a;
end proc;

Each of these will throw an error when called; which admittedly is not as good as a compile time check (or rather, in Maple that would be an evaluation-(or simplification-)of-the-procedure type check), but it is at least earlier than when three procedures down the road you use the returned value to index a list.

Of course, parameter type checking (including delimiting the argument sequence by $ if you don't expect additional arguments) falls into the same category.

Hope that this slightly off-topic response still helps,

Erik Postma

Here's what I get if I copy and paste the 2D output from Maple:

This is an image! Note I have extended typesetting turned off in Maple, so that's presumably why BesselJ shows up as BesselJ instead of as J_0. I'm not sure exactly how copy and paste works in general, but I think the "providing" program (Maple, in this case) can offer a number of different formats for the copied content, and apparently the editing widget selects the image format. It's quite a wide image though; about double the width of the expression itself. So less suitable for in running text () because of the abundant amount of whitespace.

And here is what I get if I type or paste into the popup window:

BesselJ(0, x)+GAMMA(Psi)-sqrt(abs(omega))

It doesn't hang for me, but it didn't show anything useful in the little preview thingy until I also typed something into the box (a space, which I subsequently removed again). (That's probably because I pasted using the middle mouse button; when I retried using Ctrl-V it did update immediately.)

Note that the first image doesn't seem to have any useful "alt" text; the second one does. I concur with Alec and Acer that the second image looks horrible; the first is quite a bit better.

By the way, this is on Firefox 3.5.9 in Ubuntu 9.10.

Erik Postma

As in many other questions (and posts?), the first thing I am inclined to ask is: why do you want this, and could you do without it?

Maybe efficiency is of no concern, or the number of lists is always known beforehand to be very small. But if that is not the case, then Maple's lists are inefficient for building up incrementally; moreover, there is really no way to get around that if you need to keep all the intermediate results c[1], c[2], ... . But I find it hard to imagine a scenario where that would be absolutely necessary. For example, if you need access to only some of these partial lists after you have run the loop, then it might be better to just construct those partial lists on the fly: wherever in your program you had c[i], use seq(op(a[j]), j = 1 .. i) locally there. This is still not very efficient, but if you need only a few of these lists, it might be good enough.

If you need to do something to many (or all) partial lists, then the best solution might be to just store your data differently. For a few releases now, we have had automatically growing Array data structures (see the very last paragraph of ?rtable_indexing), which are very convenient to use; if you store all the entries of a[1], a[2], ..., in one big Array and use a separate Array (with datatype = integer[4] for bonus points) to store the last index of the elements belonging conceptually to a[1], a[2], etc., you can get much better performance. So the lists a[1] = [1, 2, 3], a[2] = [4, 5], a[3] = [], a[4] = [9, 4] would be represented with A1 := Array([1, 2, 3, 4, 5, 9, 4]) and A2 := Array([3, 5, 5, 7], datatype = integer[4]). You would then have access to the equivalent of c[i] by taking A1[1 .. A2[i]], or even better, ArrayTools[Alias](A1, [1 .. A2[i]]) (which would not involve any copying at all), or just A1 and knowing that you need to look only up to index A2[i]. Of course, all of this would result in an Array rather than a list, so you'd need to adapt the algorithm that operates on c[i].

Hope this helps,

Erik Postma

The easiest way to do that is by using a different command than Histogram: the  ?Statistics/TallyInto  command, which Histogram uses internally anyway:

a := [1, 3, 4, 6, 7, 8, 9, 4, 5, 9, 8];
minmax := min(a) - 1/2 .. max(a) + 1/2;
frequencies := map(rhs, Statistics[TallyInto](a, minmax, bins = (rhs - lhs)(minmax)));

This is much more reliable than obtaining the values from the Histogram itself, because the exact plot data structure that is used could be different in the next Maple release. However, if you need to do it, it can be done as follows (no guarantees about the future, of course):

valuesFromHistogram := proc(plotStructure :: specfunc(anything, PLOT), $)
local polygons, quadrangles;
  polygons := indets(plotStructure, 'specfunc'('anything', 'POLYGONS'));
  if nops(polygons) <> 1 then
    error "expected 1 POLYGONS structure, found %1", nops(polygons);
  end if;
  quadrangles := select(type, [op(polygons[1])], 'list'(['numeric', 'numeric']));
  return map(max, map[3](map2, op, 2, quadrangles));
end proc;


Hope this helps,

Erik Postma

Hi Nick,

One issue with what you posted; this may be just a typo in copying things to the site. You write step = 0.1, which is just an equation, not an assignment; so Maple evaluates that equation but doesn't do anything with it. That means that in the UpdatePosition procedure, Maple calculates expressions involving step as a symbolic variable (since it doesn't have a value assigned to it). That's relatively expensive.

If you replace the equation with the assignment step := 0.1; you already get a speedup; for me, the loop takes 0.00135 seconds (measured by running it 1000 times). But using the Compiler[Compile] command to compile a (slightly rewritten) version into C and run that, I got it down to 0.000039 seconds (measured running it 10^4 times). Here's what I did:

pos := Matrix(3, 241, datatype=float[8]); # Need datatype=float[8] for the compiler to be able to deal with it.
vel := Matrix(3, 241, datatype=float[8]);

step := 0.1;

pos(1, 1) := 2; pos(2, 1) := 0; pos(3, 1) := 5;  #arbitrary values
vel(1, 1) := 13; vel(2, 1) := 0; vel(3, 1) := 9;

UpdateAllPositions := proc (maxi :: posint, pos :: Matrix, vel :: Matrix, step :: float[8])
local i :: posint;
  for i from 2 to maxi do
pos[1, i] := pos[1, i-1]+vel[1, i]*step; #x position
pos[2, i] := pos[2, i-1]+vel[2, i]*step; #y position
pos[3, i] := pos[3, i-1]+vel[3, i]*step  #z position
end do:
end proc;

UpdateAllPositionsC := Compiler:-Compile(UpdateAllPositions):

ttt := time():
to 10000 do
UpdateAllPositionsC(241, pos, vel, step)
end do:
(time() - ttt)/10000;

By the way, since the velocities for i > 0 are all zero, there is an even faster solution - just copy the first column to all other columns (probably quickest if you use ArrayTools). But I guess that would be considered cheating, right?

A third option would be to use LinearAlgebra to do these operations; anyone want to try that?

Hope this helps,

Erik Postma

Edit: I meant to add a fourth option; if what you are modelling is essentially a discretized differential equation (as suggested by the position and velocity names for the matrices), possibly with discrete switching, it might be worth investigating the use of ?dsolve/numeric. It has had a very sophisticated facility for dealing with such discrete events, explained in ?dsolve/numeric/Events, and provides a very natural, accurate, and efficient way to deal with such stuff.

I tried a few random (positive) values for the other parameters, and for a couple of them, there was a unique real solution for i, but not with these values:

y = 2
a = 9
b = 1/10
c = 1/10
d = 5
g = 3/4
h = 1/8

Here, there are still a number of solutions (Maple finds 51 of them), but they're all complex, and I guess that would not be valid in the context of neural networks. Here's a plot of the right hand side of the equation minus y:

plot(-2+3/4/(1+exp(-9-i/10))+1/8/(1+exp(-1/10-5*i)), i=-infinity..infinity)

And indeed, we can see that if y > g + h, then substituting anything positive for the exponentials is not going to lead to a solution. So this seems to indicate that the function is not invertible -- unless this violates some constraint on the values that the parameters should have?


There are two degrees of freedom involved here, expressible through three quantities: the number of bars, the width of each bar, and their product: the width of the range of data that is covered. If you use the ?Histogram command from the Statistics package, then you can either set these explicitly or select a rule of thumb that Maple should use to select these parameters; the default is one of the rules. In more detail:

  • The range option is set to deduce by default, which includes all given data points. You can select a smaller (or indeed larger) range as well.
  • The bincount option selects the number of "bins" the data is collected into - in other words, the number of bars. If this option is given, then the next option is ignored:
  • The binwidth option selects either the absolute width of a bar (specified as a positive number), or the bin width calculation method as one of the keywords Sturges, Scott, or FreedmanDiaconis. The default is Sturges.

So if you want to have a fixed bar width, just specify the binwidth option and don't specify the bincount option.

Hope this helps,

Erik Postma

If I do 'plot'(sin(t^2), t=`0`*`..`*`3`, labels=[sqrt(time), amplitude]), then I get both labels:

plot(sin(t^2), t=0..3, labels=[sqrt(time),amplitude])


Hope this helps,

Erik Postma

Hi Chewyraver,

I haven't tried to manipulate your equation for very long, but I think it's not unlikely that there is no closed form solution for i. So that means that you will have to make do with the RootOf expression.

As you can see on the ?RootOf help page, RootOf is a general notation for "the solution of this equation". It may sound a bit like cheating to return such a general notation from solve ("What is the solution to this equation?" "The solution to this equation is the solution to this equation."), but in many cases there just doesn't exist a combination of functions that have a name and have been studied in the mathematical literature that is an inverse to the expression you give solve. This may well be a case like that. And in fact, even if there is a function with that name, that only means that some of the properties of that function are known. Take for example the equation z*exp(z)=y; the solution for z in terms of y is given by Lambert's W function, as z=LambertW(y). From a purely algebraic point of view, you could say that this is the definition of Lambert's W function, and the solution would be nothing more than a tautology. In some sense, the same is true for the square root, which could be defined as (one branch of) the inverse function of the square. So knowing to call the solution LambertW(y) is only useful because there might be other algorithms that can take advantage of that fact - for example, we know the derivative of LambertW, and we can evaluate it numerically.

However, most Maple commands can handle these RootOf expressions fairly well; in some sense, they are just another specially named function. You can differentiate them, and generally you can evaluate them numerically as well - provided you substitute values for the other parameters, as would be true for the LambertW function as well. You write that you don't know how to handle them; the natural question is then: what do you want to do with them? If you let us know, we might be able to help you.

Hope this helps,

Erik Postma

P.S.: you know that to enter an exponential, you need to type exp(...) or use the button on the expression palette, right? If you just type "e ^ x", then you get the variable e raised to the x'th power, but e is just another variable, not the base of the natural logarithm.

Hi vuishl,

Restart is probably not the best solution, since it only works at the top level (see ?restart) - that is, not inside a loop. But there's a framework that you can use to run these sorts of things much more efficiently, if the objective is to plot: ?dsolve/numeric/Events, which uses the numeric ODE solver instead of the exact one.

I'm trying to reverse engineer what you're doing with the RootFinding:-NextZero etc., and I think this is more or less what happens: you want to follow a solution until x(t) is plus or minus one, then negate the derivative of the x(t). Is that correct? Then a more efficient way to generate such a series of plots starts as follows, in a fresh session:

ODE := diff(diff(x(t), t), t)+.2*(diff(x(t), t))-x(t) = aa*sin(t);
ICs := x(0)=mm, D(x)(0)=nn;
solution := dsolve([ODE, ICs], numeric, parameters=[aa, mm, nn],
events=[[x(t)^2-1, diff(x(t), t)=-diff(x(t), t)]], range=0..100):

Now solution is a procedure that contains everything necessary for simulation of your system, including the "bouncing back" of the derivative. The parameters aa, mm, and nn are placeholders that you need to substitute values for before simulation starts (I'll show how in a second). Then you can request the values at a particular point in time, say t = 5, as follows:

solution(parameters = [aa = 1.5, mm = 0.999, nn = -2.5]);

This will return a list of equations saying that (approximately) t = 5, x(t) = -0.56, and diff(x(t), t) = -1.09. But much more interesting is a routine that will do the phase plot for you:

plots:-odeplot(solution, [x(t), diff(x(t), t)], t=0..100, refine = 1);

This we can use in a loop.

for a from 1.0 to 2.0 by 0.4 do
  for m from 0.990 to 1 by 0.004 do
    for n from -3.0 to -2.0 by 0.4 do
      solution(parameters = [aa = a, mm = m, nn = n]);
      plotsetup(ps, plotoutput = sprintf("c:/tmp/partA[%a]M[%a]N[%a].eps", a, m, n), plotoptions = `color,portrait,noborder,width=16cm,height=12cm`);
      plots:-odeplot(solution, [x(t), diff(x(t), t)], t=0..100, refine=1);
    end do;
  end do;
end do:

As you see, there is no call to dsolve inside the loop - that's the beauty of using parameters! All that happens is that the simulation is run and the plot is drawn. Note that the name of the loop variable is different from the name of the parameter - that's always a good idea.

I would try to see if you can run this loop on your older system first. If that doesn't work, we can still break it up into smaller pieces, but I'd think that this is going to be much more efficient than what you had before.

One note: your original version always counted to 200 "bounces", whereas I bounded the simulation by letting it always run to t=100 irrespective of the number of bounces. I think we can make the simulation always run to 200 bounces, if you'd strongly prefer that; it would be a bit more complicated though. And we'd need a reasonable upper bound on how long that's going to take. Let us know if that is necessary.

Hope this helps,

Erik Postma

5 6 7 8 9 10 11 Page 7 of 11