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+.\
466299938140965520e-9*t^3.)

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
Maplesoft.

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: )
igcd(op(coefficients));
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
Maplesoft.

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
Maplesoft.

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!

Best,

Erik Postma
Maplesoft.

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.

Best,

Erik Postma
Maplesoft.

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?)

Best,

Erik Postma
Maplesoft.

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):
p1;

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
Maplesoft.

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)]);

Best,

Erik Postma
Maplesoft.

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
Maplesoft.

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

Erik Postma
Maplesoft.

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
Maplesoft.

Just call cmaple (in your bin.win 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\bin.win>copy bla.mpl con
1 + 1;
(x+1)^2 - (x^2 + 2*x + 1);
expand(%);
        1 file(s) copied.

C:\Program Files\Maple 13\bin.win>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;
                                       2

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

> expand(%);
                                       0

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

C:\Program Files\Maple 13\bin.win>

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
Maplesoft.

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))
+140*(1/10)+((1-t(x))*(1/14)+1/2-1/2*erf(x))*((1-t(x))*(1/14)
+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:

solution(parameters=[7.5]);
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]);
  try
    solution(5);
  catch:
  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]);
  try
    solution(0);
  catch:
    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
Maplesoft.

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
Maplesoft

"unable to delimit strings / identifiers" sounds like there is a problem with your use of double quotes (for strings) or backquotes (for identifiers). This could be apparent from the single line that throws the error.

You write "RETURN". Are you using the old RETURN function? That has been deprecated for a very long time; we recommend use of the "return" statement now. See ?return.  This may be completely unrelated to your error message, though.

Hope this helps,

Erik Postma
Maplesoft.