## 1394 Reputation

13 years, 180 days
Maplesoft

## Social Networks and Content at Maplesoft.com

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.

## Some commands...

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.

## Ok - let's try...

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.

## Two questions...

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.

## Questions...

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.

## A possible solution...

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.

## Alternatives...

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.

## mtaylor...

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.

## It looks like Mean(A^2/(1+A)) is the cul...

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

Erik Postma
Maplesoft.

## Four issues...

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

## Partial solution...

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.

## the discont option...

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

## error message...

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

## The other option:...

In the interest of completeness, I believe the other option would be as follows:

floatProcessor := proc(x)
local mantissa, exponent;
mantissa, exponent := op(x);
if abs(mantissa) < 10 and exponent = 0 then
return Float(mantissa * 10, -1);
else
return x;
end if;
end proc;
a := x - 3.;

subsindets(a, float, floatProcessor);

This returns x - 3.0. Hope this helps!

Erik Postma
Maplesoft

Hi Jonas,

For the first, I'm not sure why you set f(x) + f'(x) equal to zero. I would approach it as follows. If we take a tangent at x0, then the value of that tangent at x is f(x0) + (x - x0) * f'(x0). This should be equal to -3 for x=-1: so we need to solve f(x0) + (-1 - x0)*f'(x0) = -3. This will give you two solutions, which both lead to a correct answer.

Also when you obtain 1/2*x^2 + 3*x + 5, and try to solve that, I think you obtain an incorrect answer. Indeed, if you plot that expression, you'll see that it does not have any roots (which confirms the solutions -3 + I, -3 - I). The formula you use for solving the quadric doesn't seem to be correct: the determinant should be "b^2 - 4*a*c", which is in your case 3^2-4*1/2*5 = 9-10=-1.

Then on to your question about regression. The ExponentialFit command that you show gives the same result that you give for f(x), although it's written slightly differently. I could not find a very easy way to change the exp(0.043*x) into 1.04^x, by the way; I'm sure some of the experts in the forum can come up with something though. If desperate, you can always do something like this:

Q := ExponentialFit(X, Y, x);
expr := indets(Q, specfunc(anything, exp))[1];
newQ := subs(expr = exp(op([1, 1], expr))^op([1, 2], expr), Q);

For the value of R^2: what's your definition, exactly? I tried several things, but could not get the number 0.9596 that you suggest. This is what I would think gives the correct value (taking some cues from Wikipedia's page about the coefficient of determination):

QFunction := unapply(Q, x);
QValues := QFunction~(X);
option1 := Correlation(Y, QValues)^2; # returns 0.9831 (approx)
option2 := Correlation(ln~(Y), ln~(QValues))^2; # returns 0.9963 (approx)

Here I liberally use Maple 13's elementwise syntax with the tilde operator.

Hope this helps,

Erik Postma
Maplesoft.

 First 7 8 9 10 11 Page 9 of 11
﻿