epostma

1499 Reputation

19 Badges

16 years, 126 days
Maplesoft

Social Networks and Content at Maplesoft.com

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

This is because when you enter lsol, Maple will evaluate lsol and show the result, rather than just showing what is assigned to lsol. If you want to see just what is assigned to lsol, you need to enter eval(lsol, 1). This will show the unevaluated call to inttrans/laplace. So this aspect of your question is expected.

It is, however, a bug that the unevaluated call is there in the first place. It can be triggered with the simpler call:

inttrans:-laplace({f(t)=0, g(t)=0}, t, s);

Thank you for reporting this issue.

Hi tomtomo,

Maple has two different facilities for optimization: one that works symbolically (the minimize and maximize commands), and one that works numerically (the Optimization package). Only the numerical optimization package understands units.

Indeed, if you replace minimize with Optimization:-Minimize, and do a little surgery to the line after it to get the parameter assignments to work correctly, you get an answer. You need to change the assign command on the next line to just "assign(%[2])" to make it work with the simpler output from Optimization:-Minimize. Moreover, you also need to change the visualization commands a bit to work with units. However, once you've done this, you'll see that the fit is not really great. I think this is a case where you aren't really guaranteed a good fit without a global rather than local optimization solver. This does not come with Maple (though there are third party solutions, and indeed even a Maplesoft toolbox, that do that).

That said, as I was writing this up, I noticed a mistake earlier on. You use the convert/radians command, together with Unit(degrees). They do not work together, and indeed there should be at least a warning to that effect in both the help page and the command itself. This ends up causing *two* conversions from degrees to radians, which is not what we need. Instead, you can just use '90*Unit(degrees)' there.

I reworked that section of your worksheet a bit, which ended up with the following -- and note that the fit in this case ends up being very good. To make it self-contained, I'll include the definition of the 'data' Matrix. Some of the things I do manually below happen automatically when you load the Units:-Simple package; I would recommend that, actually, but figured I'd show how to do it without that.

data := Matrix([
[81*Unit(knot), 1600*Unit(ft), 0, 89*Unit(knot), 179*Unit(degrees), 179*Unit(degC)],
[81*Unit(knot), 1600* Unit(ft), 0, 77*Unit(knot), 303*Unit(degrees), 303* Unit(degC)],
[81*Unit(knot), 1600*Unit(ft), 0, 95*Unit (knot), 72*Unit(degrees), 72*Unit(degC)],
[81*Unit(knot), 1600*Unit(ft), 0, 98*Unit(knot), 95*Unit(degrees), 95* Unit(degC)]]);

F := a*x + b*y + x^2 + y^2 + c;
XY := Matrix(4, 2);
for i to 4 do
X := evalf(data(i, 4)*cos(-(data(i, 5) - 90*Unit(degrees))));
Y := evalf(data(i, 4)*sin(-(data(i, 5) - 90*Unit(degrees))));
XY[i, 1] := X;
XY[i, 2] := Y;
end do;

Optimization:-Minimize(add(eval(F, {x = XY[n, 1], y = XY[n, 2]})^2, n = 1 .. 4));
# result: approx.  [1941 * Unit(m^4/s^4), [a = -10.5*Unit(m/s), b = 2.8*Unit(m/s), c = -1971*Unit(m^2/s^2)]]

assign(%[2]);

# Original: TAS := evalf(sqrt(c)); # this shows the unit 'sqrt(m^2/s^2)'; to simplify that, use combine(..., units) - and you don't need the evalf here, it's automatic.
TAS := combine(sqrt(c), units);

xC := 0.5*a;
yC := 0.5*b;

wind_spd := sqrt(xC^2 + yC^2); # result: about 5.43*Unit(m/s)
# Original: wind_dir := arctan(yC/xC); # because Maple doesn't know where the negative sign comes from after the division, you don't find the correct quadrant -- instead, use 2-argument arctan:
wind_dir := arctan(yC, xC); # result: 2.88 (implied unit: radians)

# Need to evaluate the sines and cosines; this would happen automatically with Units:-Simple
XY_combined := map(combine, XY, units);
A := plots:-pointplot(XY_combined);

# implicitplot doesn't know how to handle units, so we need to strip them out and deal with them manually
F_no_units := eval(F, Units:-Unit = 1);
B := plots:-implicitplot(F_no_units, x=-180..180, y=-180..180, color=green);

plots:-display(A, B);


Hope this helps!

Erik Postma
Manager, mathematical software group
Maplesoft.

It looks like the DataTable component doesn't always store attributes of function calls inside it. When I do:

op(2, newtable["M"]); # shows Unit(kN * m)
[attributes(op(2, newtable["M"]))]; # shows an empty list
[attributes(Unit(kN * m))]; # shows: [Units:-UnitStruct(1000, newton, SI)*Units:-UnitStruct(1, metre, SI), inert

For newtable["F"], the correct value is shown. This is a bug in the code dealing with data tables which I will enter into our bug tracking system -- thank you for finding it!

The attributes are what's used in the Units package for doing the actual manipulation of units, including the arithmetic implemented in Units:-Simple. So that's running into the absence of this attribute.

As a workaround, if you do eval(newtable["M"]), that re-evaluates the Unit function call, which re-computes the correct attribute.

Hi everyone,

This is indeed a bug. Maple was not careful enough to avoid overflow when trying to evaluate the probability mass function in evalhf mode. It needs to be evaluated in regular evalf instead, or use the cleverer routine that uses the external code implementation.

I'm currently testing a fix, which, if everything goes as expected, should make it into the next major version of Maple. Thanks Nikol for finding it!

Erik Postma
Manager, mathematical software group

I agree; the documentation is not precise enough here. I'll change it to say that the element is inserted after position n, which is less likely to be misunderstood.

I just submitted a change into our code repository that fixes this issue. Unless we find a major unexpected problem caused by this fix, the issue should disappear in any future major versions of Maple.

I also ran Carl's sheet, and it shows a couple of single asterisks (as you would expect) and two double asterisks that disappear on rerunning that case. I think that can be considered passing the test.

Finally, just out of interest, I traced which distributions' custom samplers (potentially) depend on the normal distribution one (other than just reduction in a special case); it's quite a list. I figured some of you might be interested:

  • the gamma, log normal, student's T, and non-central versions of chi square and student's T directly call the normal distribution sampler;
  • the beta, chi square, and error distribution call the gamma distribution sampler;
  • and the F ratio, inverse gaussian, and non-central F ratio and beta distributions call the chi square distribution sampler.

Hi mmcdara,

I'm afraid I can't reproduce your finding that the Student[Statistics] package does this better than the top level Statistics package; in Maple 2018, if I do:

with(Student[Statistics]):
R := NegativeBinomialRandomVariable(1.965, 0.003);
Quantile(R, 0.98);

I get no answer, just like with the top level Statistics package. In fact, that is to be expected: if you run the command showstat(Student[Statistics][Quantile]), you can see that it directly calls the Quantile command in the top level Statistics package. (That's how the Student[Statistics] package is set up in general - it's meant to be a slightly simpler user interface to the same functionality. Among the few bits of functionality that are exclusive to Student[Statistics] are plots associated with quantities (activated with the output option) and the TestsGuide command).

Maybe more importantly, as I mention in my response to Teep's post, I have a fix to this bug in the works.

Erik Postma
Manager, mathematical software group.

Hi Teep,

Thank you for this report; I saw it today and I have a fix in progress. It looks like the biggest problem is that Maple uses code here that is meant for continuous distributions instead of discrete distributions, for doing root finding, trying to find the value where the CDF crosses the value 0.98. The code fairly quickly finds that at 1919.998, the CDF is 0.9799689, and at 1920.42, the CDF is 0.9800204, but then gets stuck trying to find an epsilon-sized interval in between that contains it. It then tries an alternate method that is a very bad idea.

With this fix, Maple will use a discrete interpolation method for discrete distributions. We've had the code for this since the Statistics package was first written, but it wasn't used in this branch of code. The fix should make it into Maple 2019, unless something unexpected happens.

 

Erik Postma
Manager, mathematical software group.

Thanks for this report, _Maxim_. This issue appears to be fixed in our internal development version. I can't promise any particular future version the fix will show up in, but it should make it out eventually!


Thanks for this report, _Maxim_ - it looks like this issue is fixed in Maple 2018.


As a (partial?) workaround, there is

dsolve(diff(f(z), z) - f(z)/z = 0, f(z), series);

which returns  _C1*z*(1 + O(z^6)) with the default setting for Order.

Unfortunately, this is a known limitation in the current versions of Maple. There are various stages of evaluation of a module definition, and at the time the option option object(ParentObject) is being processed, the ParentObject cannot yet be recognized as a valid object. I am not sure we can fix this in the near future, and all the potential workarounds that I'm aware of have significant downsides.

One such workaround is of course to make the objects top-level entities, instead of living inside a module. In this case, you lose the name space separation that a module provides.

Another would be to define the objects after defining the package:

module Package()
  option package;
  export ParentObject, ChildObject;
end module:

unprotect('Package:-ParentObject');
Package:-ParentObject := module()
  option object;
end module:
protect('Package:-ParentObject');

unprotect('Package:-ChildObject');
Package:-ChildObject := module()
  option object(Package:-ParentObject);
end module:
protect('Package:-ChildObject');

This works well if there is nothing else inside the package, but:

  1. It is clumsy with the protect / unprotect pairs.
  2. If there are other members of Package, then you don't have lexical access to them. In particular, locals of Package are invisible (unless you use kernelopts(opaquemodules=false), which has its own problems).
  3. As a consequence, it requires that all classes that you define are exports of Package. Otherwise, you might for example want to have a common parent class of a number of derived classes be a local that's invisible to the outside world.

Erik Postma
Manager, mathematical software group.


Hi tbfl1919,

 

Let's first investigate what exactly is going on.

The functions V__tot and similar all take a unitless quantity t as their argument, as you can see by evaluating, for example, V__tot(1.3). In the integral, however, you specify the lower and upper bounds for t, the argument of V__tot, as T__1 and T__2, which are 0 and 1 ms. (I should note here that 0 is a special case - it can only be unitless in Maple for technical reasons - but the int command interprets the range as expressed in units of time.) So the int command tries to substitute time values for t, which causes the sin command inside V__tot to fail, as it does on the command line if you evaluate V__tot(1.3*Unit(ms)). Maple decides that that's an expression it can't integrate and merely takes the unit, V2, out of the integral. So we have now evaluated the integral to some unevaluated integral times the unit V2. This is then divided by T__2 - T__1, a duration, to obtain an unevaluated integral with the unit V2/s. The square root of this then has the unit V/s1/2, which is equivalent to the unit you see. It could be argued that this is correct, because if we could specify the integrand better, the integral would come out to be a duration (it is a unitless quantity, integrated over a period of time, so the dimension is time), its square root therefore has the unit s1/2, and multiplying that with the unit we see gives the unit m2 kg / s3 / A, which is the same as V.

 

However, this is all of course incredibly clunky and it doesn't lead to Maple giving you the answer you're looking for. So let's talk about solutions. The first option is to make the integral strip off the time unit and add it back on manually. I think the appropriate unit is ms, correct? So we can do that by specifying the lower and upper bounds as T__1/Unit(ms) and T__2/Unit(ms) and multiplying the integral by Unit(ms). This appears to lead to a correct result (of sqrt(2)/2 V). But this is, in some sense, cheating: we just remove the units manually and then put them back on. It's possible to do better.

 

The second option is to make V__tot accept an expression in terms of time, rather than a unitless quantity. For this we need to do two things:

  1. Redefine V__tot_AC(t) := A * sin(omega * t);
  2. Be very careful to never submit a unitless quantity into V__tot or V__tot_AC.

The latter issue appears, for example, when we try to plot these functions: the only way I can think to do this is to change the command to plot(V__tot(t * Unit(ms)), t = 0 .. 5) (and similar for V__tot_AC). It also appears when we want to do the integration: just entering V__tot(t) already causes an error. However, we can prevent the error by using unevaluation quotes (single forward quotes) around V__tot(t): int('V__tot(t)'2, t = T__1 .. T__2) evaluates no problem. This leads to the same result of sqrt(2)/2 V.

 

You don't write what version of Maple you are using; if it is version 2017, you can use another trick to get around some of the problems here. In 2017 we introduced the Units[Simple] package as an alternative to Units[Standard]. The most important difference between the two is the following: the Units[Standard] package assumes that every unassigned variable, such as t, represents a unitless quantity. The Units[Simple] package does not. This restriction may have played a role in your decision to define V__tot_AC(t) the way you did, where you included the unit for t in the definition, rather than it just being part of the value t that the user submits. Indeed, with the Units[Standard] package, you cannot evaluate sin(omega * t), if t is unassigned and omega has a unit of frequency, because the argument to the sine function is then assumed to have a unit of frequency, which is illegal. On the other hand, if you load the Units[Simple] package, it will see the unassigned variable t and not assume anything about it. (It does assume that every unassigned variable has a single well-defined dimension within one expression; if you evaluate (t + 5 kg) * (t + 5 ms), then Maple will complain: the first factor implies that t must be a mass and the second that t must be a duration.) So compared with the Units[Standard] solution, you don't need to worry about what arguments you submit to V__tot so much: if the argument can be a duration for some assignment of units to your variables, then Maple will accept it. So item 2. above is not so important anymore. In particular, you get the correct answer for V__tot_rms by only removing the unit ms from the definition of V__tot_AC, no quoting needed anywhere.

 

Finally, for the upcoming version of Maple (which is still some time away!) we have overhauled plotting with units somewhat; if you use my suggested definition of V__tot and you use Units[Simple] so that you don't need to worry about evaluating V__tot on symbolic arguments, you will be able to get the correct plot just by writing plot(V__tot(t), t = 0*Unit(ms) .. 5*Unit(ms)) or, if you prefer, plot(V__tot, 0*Unit(ms) .. 5*Unit(ms)). We feel this is a more natural way to use units in plots.

If you replace the last loop with:

e2X := 0;
for a from 1 to 1 do
  for alpha from 0 to k-1 do
    for b from 0 to k-1 do
      for beta from 0 to k-1 do
        f := 2*S*(d[beta+1,alpha+1]*W(beta))*L(alpha)*d[a+1,b+1]*L(b);
        e2X := e2X - f;
      end do;
    end do;
  end do;
end do:
e2 := int(e2X, z = -h/2..h/2);

then you save about half of the running time.

Kitonum has a good suggestion in their answer. Here is an alternative: the typematch command. You pass it the expression to be tested, and a pattern. This pattern consists of type names with a double colon and a name. For example, to match a product of two sums, each with two terms, you could match to the expression (a::algebraic + b::algebraic) * (c::algebraic + d::algebraic). The command returns true if and only if the expression matches the pattern, and in that case it assigns the matching values to the variables occurring in the pattern; in this case to a,b,c,d. See the documentation here: https://www.maplesoft.com/support/help/Maple/view.aspx?path=typematch. One thing that's explained there is how to deal with evaluation problems if the variables in your pattern are assigned values, e.g. if you've used the same pattern successfully before. (It involves unevaluation quotes.)

 

If you don't care for the matching values, just for the yes-or-no result, you can use the type command. It takes a structured type as its second argument, which is basically like the second argument of typematch without the variables. A tricky thing here is that Maple automatically simplifies something like algebraic+algebraic to 2*algebraic, and you don't want that to happen with your pattern.

 

In either case, you may need to study structured types a bit too come up with the right pattern/type. They are explained (rather succinctly) on these help pages: http://www.maplesoft.com/support/help/Maple/view.aspx?path=type/structure and http://www.maplesoft.com/support/help/Maple/view.aspx?path=type. Let us know if you have trouble coming up with the right type to check!

1 2 3 4 5 6 7 Last Page 1 of 11