1494 Reputation

19 Badges

15 years, 315 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

I'm not sure why you would want to construct a loop, specifically. The easiest way is to just use solve (or fsolve for floating point answers).

m := <x, x^2; 3*x, 10*x>:
det := LinearAlgebra:-Determinant(m);
solve(det = 0, x);
fsolve(det, x);

Hope this helps,

Erik Postma

I'm guessing that 'Fit' just returns unevaluated? The problem is probably with your "with". I see there's a space between "with" and the opening parenthesis. If you enter this in 2D math mode, Maple will interpret this as multiplying the variable "with" by the variable "Statistics". Get rid of the space and it should work. (It did for me; I got .850632001360003054+.500576419548484901e-2*t-.247858474516350479e-5*t^2+.\

By the way, as a shortcut, you should be able to call

PolynomialFit(3, Temperatures, BondEnergy, t);

without constructing the polynomial explicitly.

Hope this helps,

Erik Postma

Well... if you want to represent the numerator and denominator exactly, there's not much you can do. The coefficients are coprime, so there's no common factor to factor out. However, you could divide by the minimal coefficient, then approximate numerically:

# This is just your last line, assigning it to ratfun:
ratfun := CurveFitting[RationalInterpolation](XX, YY, x, degrees = [top_order, bottm_order]);
coefficients := [coeffs(expand(numer(ratfun))), coeffs(expand(denom(ratfun)))]:
# (you could verify that the gcd is one now: )
floatratfun := evalf(expand(numer(ratfun)) / min(coefficients)) / evalf(expand(denom(ratfun)) / min(coefficients));

Or you could even take the coefficients divided by the minimal coefficient (or otherwise scaled) and take truncated continued fractions. That would be a worse approximation than floating point numbers (probably), but it would leave you with an exact rational function.

Hope this helps,

Erik Postma

Hi HerClau,

Beyond ?factor , ?expand , ?collect , and ?normal (which I all used in the reply above), ?combine is also a standard tool for this sort of expression manipulation. (The links above are all to the online help pages.) I currently cannot think of any other commands that fall into the same category.

Hope this helps,

Erik Postma

I define eq1, eq2, eq3 as above. The first step in the paper is to substitute the latter two into the first:

eval(eq1, [eq2, eq3]);

Then we isolate m[ct]:

eq4 := isolate(%, m[ct]);

This leads to a slightly different expression than what you obtained (equivalent of course, but different still); we can get (essentially) the same thing by setting

expr5 := map(factor @ normal, collect(rhs(eq4), {`&Deltaw`, m[cr]}));

We can get at the term containing m[cr] like this:

term6, term7 := selectremove(has, expr5, m[cr]);

and do the rewritings that result in Ec. 3 as follows:

term8 := m[cr] * (1 + factor(normal(term6/m[cr] - 1)));

so that we can reproduce Ec. 3 like this:

eq9 := lhs(eq4) = term7 + term8;

Note some of the minus signs are chosen differently; I assume that's OK.
Now more happens to what we called term8: there's approximations. The way I would do this is to first do the mathematical operation (substitute rho[0] = rho[a] = 0 into the denominators) and then massage the result to look the way you'd like it to look. The operation "substitute into denominators only" is of course sufficiently rare that it's not built in to maple, but we can perform it using ?evalindets . So the first step is then:

term10 := evalindets(term8, anything ^ negative, eval, [rho[0] = 0, rho[a] = 0]);

In order to change the factor (rho[r] - rho[t])/rho[r]/rho[t] into 1/rho[t] - 1/rho[r], I would like to use expand; but we would need to apply it only to that subexpression. This is somewhat tricky.

fraction11 := term10 / m[cr] - 1;
fraction12 := expand(normal(fraction11 / (rho[a] - rho[0]))) * (rho[a] - rho[0]);

Now fraction12 corresponds to your variable C, and our reasoning shows that eq9 is equivalent to:

eq13 := lhs(eq9) = term7 + m[cr] * (1 + 'fraction12');

I'm not entirely sure how you make "term7" reduce to just `&Deltaw`, but I hope that this helps you along for a bit. Feel free to ask for more clarification!


Erik Postma

Let's see - I can more or less follow what you're doing in the paper (I'm not "checking every step", but I mean I see what you're doing, more or less), but I'm not sure what you're asking here - do you want to reproduce the result from the paper in Maple?

Another question: you specify three "input" equations,

eq1 := m[t]*(1-rho[a]/rho[t]) = m[r]*(1-rho[a]/rho[r])+`&Deltaw`;
eq2 := m[t] = m[ct]*(1-rho[0]/rho[c])/(1-rho[0]/rho[t]);
eq3 := m[r] = m[cr]*(1-rho[0]/rho[c])/(1-rho[0]/rho[r]);

The quantity m[a] does not appear here, but you say you need to get an equation involving m[a] (viz., m[a] = m[a]*(1+(rho[a]-rho[0])*(1/rho[t]-1/rho[r]))+`&Deltaw`, which you write using a variable C - do I understand this correctly?). Again I'm not sure what exactly you mean - do you want to derive this equation from eq1, eq2, eq3? That will not succeed, since m[a] does not appear in them.


Erik Postma

What's the "data type" of the index and of the superscript? The "s > 1" suggests that the subscript is a real number (maybe even integer?). Can you represent the superscripts (the generators) in a meaningful way? (E.g. as vectors over some field?)

Also: it looks like the multiplication that you wrote down is commutative, not anticommutative. (Unless the + you use in the sub- or superscript is intended to indicate a noncommutative operation?)


Erik Postma

Hi Jeremy,

I think the easiest solution may be the following.

First create a plot where lon and lat are regular Cartesian axes. (I'm using a silly quantity here, but you get the idea.)

p1 := plots[densityplot](lat * lon^2, lon = -Pi .. Pi, lat = -Pi/2 .. Pi/2):

You could also use a contour plot (with plots[contourplot], probably using filledregions=true). Then define the mapping from Cartesian to Aitoff:

mapping := (lon, lat) -> [2 * sqrt(2) * cos(lat) * sin(lon/2) / sqrt(1 + cos(lat) * cos(lon/2)),
                          sqrt(2) * sin(lat) / sqrt(1 + cos(lat) * cos(lon/2))];

And now the big trick: you can define a mapping that transforms a plot based on a mapping that transforms the coordinates. You can do this using plottools[transform], as follows:

plotmapping := plottools[transform](mapping):
plots[display](plotmapping(p1), scaling=constrained);

Hope this helps,

Erik Postma

Hi Magdalena,

I realize your question has already been solved, but for the record here's two more alternatives. You didn't specify if w is a list or some other data type; I'll first assume it's a list somewhat like this:

h := 5;
w := [seq(<i, i + 1, i + 2>, i=1..h)];

Then two very direct ways to make the matrix are:

W1 := `<|>`(op(w));
W2 := Matrix(w);

The funny looking `<|>` function from the first alternative is the same that you use if you type <w[1] | w[2] | w[3]>. The second alternative has the advantage that you can, if you need to, specify such things as datatype and storage (see ?Matrix ). If w is not a list (or you don't want to use all elements of w), you can still use the same techniques:

W3 := `<|>`(seq(w[i], i=1..h));
W4 := Matrix([seq(w[i], i=1..h)]);


Erik Postma

You can use ?mtaylor. Now mtaylor doesn't work with functions with a multidimensional codomain, but you can use map to get around that. The notation you use for expressing F in terms of f1 and f2 is a bit ambiguous, but the following should work for both Vectors and lists:

F := [f1(x, y), f2(x, y)];
map(mtaylor, F, [x = alpha, y = beta], 3);

This will cut off the taylor series at the third order term (that's what the 3 is for).

Hope this helps,

Erik Postma

I'll enter this in our database of issues. Thanks!

Erik Postma

Hi tombfh,

I can see four issues with the example code you posted. First: you type exp^(k*cos(x-mu)). That should be exp(k*cos(x - mu)); exp is the exponentiation function, not the base of the natural logarithm itself (that would be exp(1)).

Second: you use the variables x and y to represent both the data that you have for these variables, and for the abstract variables themselves. That doesn't work. I used xd and yd for the data corresponding to x and y in the example below.

Third: I don't think your plot command works. If the stats[fit] command would have worked (more about that in a minute), then you would have gotten back an equation in x and y. You type eq(x, y), which doesn't work -- you just need eq, it's already got x and y in it. Or rather, you need the right hand side of eq: that is the expression you want to plot. (Maple automatically assumes you'll give the expression for the y-value in terms of the x-value.) Then as the second argument you supply the range for the x-variable. The example below works.

Fourth: I see you try to use the stats[fit] command, which is deprecated. (Also, it can't be called that way without first doing with(stats).) (And finally, stats[fit] can't deal with expressions that are nonlinear in the parameters.) The recommended package to use now is Statistics. There, you can use the following:

(**) xd := <0,2,3>;
(**) yd := <5,1,0.6>;
(**) eq := a*exp(k*cos(x-mu))/(2*Pi*BesselJ(0,k));
(**) result := Statistics[Fit](eq, xd, yd, x);
(**) p1 := plots:-pointplot([seq([xd[i], yd[i]], i=1..3)]);
(**) p2 := plot(result, x=min(xd) .. max(xd));
(**) plots:-display(p1, p2);

You will note that in this case, there's an exact fit. That's possible now because the number of data points is the same as the number of parameters (both are three). If you add, say, the point (4, 0.5), that doesn't happen any more.

Hope this helps,

Erik Postma

Just call cmaple (in your directory) and give the .mpl file as a command line argument. Here's what I just tried in my Windows 32 VM:

C:\Program Files\Maple 13\>copy bla.mpl con
1 + 1;
(x+1)^2 - (x^2 + 2*x + 1);
        1 file(s) copied.

C:\Program Files\Maple 13\>cmaple bla.mpl
    |\^/|     Maple 13 (IBM INTEL NT)
._|\|   |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2009
 \  MAPLE  /  All rights reserved. Maple is a trademark of
 <____ ____>  Waterloo Maple Inc.
      |       Type ? for help.
> 1 + 1;

> (x+1)^2 - (x^2 + 2*x + 1);
                                   2    2
                            (x + 1)  - x  - 2 x - 1

> expand(%);

> quit
memory used=0.7MB, alloc=0.8MB, time=0.01

C:\Program Files\Maple 13\>

Here I only typed the things at the command prompt; Maple reads the file bla.mpl and executes it and quits all by itself. (Is there a less-archaic way of obtaining the functionality of the Unix-program cat in windows than copying to the con pseudofile?)

Hope this helps,

Erik Postma

Here's what I did. First, let's define the problem ODE and ICs (we usually say ICs if they are defined at a single point, as in this case):

ode := diff(t(x), x, x)+2*x*(diff(t(x), x))
+1/2+1/2*erf(x))*exp(-1/t(x)) = 0;
ics := t(-2) = c + (c-1)*erf(-2), D(t)(-2) = 2*(c-1)/exp(4)/sqrt(Pi);

Note that I used fractions to represent fractional numbers, not floating point numbers. This is generally a very good idea. Now we compute the solution procedure, for arbitrary c. Note we can only obtain numerical values and plots once we tell this procedure what the value for c is.

solution := dsolve({ode, ics}, numeric, t(x), parameters=[c]);

For c = 7.5, we get a singularity around -1.68:

plots[odeplot](solution, -2 .. 1);

So the first task is to find a value of c for which there is no singularity between -2 and 0. We can do that as follows:

findsing := proc(cvalue)
global solution, x;
  solution('parameters' = [cvalue]);
  end try;
  return eval(x, solution('last'));
end proc;
plot(findsing, 5 .. 10);

This shows that the singularity moves up as we increase c, and it moves up more as we increase c more. That gives some hope that it might move away as we increase c enough. Experimenting with the range, we see that the singularity moves past zero between c=190 and c=200. (I tried plot(findsing, 0..100) and plot(findsing, 0..1000), and then plot(findsing, 150 .. 250) to narrow it down.)

So now we are interested in finding a c-value where t(0)=8. Let's write a procedure that turns a c-value into the value t(0).

zerovalue := proc(cvalue)
global solution, t, x;
  solution('parameters' = [cvalue]);
    return 0;
  end try;
  return eval(t(x), solution(0));
end proc;

We can "solve" this graphically by issuing the command plot(zerovalue - 8, 195 .. 200); and see that it's about 198.8, or be a bit more sophisticated and say

fsolve(zerovalue - 8, 195 .. 200);

to obtain 198.8831747. Note that the curve does not seem to have its maximum at x=0, though; and I'm not sure I understand your claim that there should be a solution that is a 'bell curve'; do you mean that there should be a solution of the form t(x) = c1 * exp((x - c2)^2 / c3) ? Because that seems to be falsified by odetest.

Hope this helps,

Erik Postma

There's a very easy way to get a similar result as in Mathematica:

plot(piecewise(x < 1, x, x <= 2, 2, 3), x = 0 .. 3, discont = true);

The discont = true option makes Maple find the discontinuities and put symbols where the value at that point is. I don't think there's an easy option to get different symbols at the other value, but you could do that using discont itself. That would require some programming, though.

Hope this helps,

Erik Postma

First 7 8 9 10 11 Page 9 of 11