acer

17418 Reputation

29 Badges

14 years, 138 days

On google groups. (comp.soft-sys.math.maple and sci.math.symbolic)

On stackoverflow.

On math.stackexchange.com.

MaplePrimes Activity


These are answers submitted by acer

In Maple syntax the equation of a plane could be entered as,

2*x+y+2*z-15=0

Or, in 2D Input mode you could put a space between a coefficient like 2 and the x, or 2 and z.

But 2x is not valid input.

You have omitted stating anything about the order being used ( eg, when calling Groebner[Basis] ). In fact you haven't said much about the kind of problem you're trying to resolve, which isn't terribly helpful.

If you read the Help page for topic Groebner,Basis_algorithms then you can read this bullet point in its Description section:

    method=default runs the default strategy, which is subject to change in the future.
    Currently this is method=direct for 'tdeg', 'wdeg', 'grlex' and 'matrix' orders and
    non-commutative problems, and method=convert for 'plex', 'lexdeg' and 'prod' orders
    in the commutative case.

And the bullet point for method=direct includes the following: (See also the bullet point for method=fgb)

    method=direct runs the fastest general purpose direct method.  This is currently
    the compiled F4 or Maple's F4 as a fallback.

Here is a link to a webpage related to the compiled F4 of J-C. Faugere. You may also see a wikipedia page for it, in which it describes parallelism in the linear algebra computation.

It's possible that you might also find the Help page for topic SolveTools,PolynomialSystem of interest. In the page's Description you can see that the option method=basis utilizes the Groebner[Basis] functionality while the option method=triade utilizes the RegularChains package's functionality. The RegularChains package incorporates an engine developed by the team led by Marc Moreno Maza.

What kind of additional, high level parallelization did you have in mind? Are you thinking of something very coarse (and as if you had resources to burn..), such as what, different orderings?

 

 

As member vv has cleverly demonstrated, in the case that the arc-length integral needs to be done numerically then solving an equivalent numeric differential equation will be far more efficient that most any naive numeric root-finding approach.

By "naive" I mean that there's no tracking of the last solutions or use of other leverage mechanisms. Only some easy memoization is done (since the root can get utilized in each entry of the Vector).

However, since fsolve was mentioned, here is one way that it may be used for the rootfinding part. This worked in both my Maple 2019.1 and Maple 2018.2 on 64bit Linux. It's still about 20 times slower than vv's dsolve/numeric approach, naturally.

kernelopts(version);

`Maple 2019.1, X86 64 LINUX, May 21 2019, Build ID 1399874`

restart;

r := t -> <cos(2*t), sin(3*t), 4*t>:

L := unapply('evalf[ceil(Digits/2)]'(Int(unapply(sqrt(add(diff(r(t),t)^~2)),t),
                                         0..tt)),tt,numeric):

f := unapply('fsolve'(s=L(t),t=0..infinity),
             s,numeric,proc_options=[remember]):

CodeTools:-Usage( plots:-spacecurve(r(f(s)), s=0..30, color=red) );

memory used=51.16MiB, alloc change=71.01MiB, cpu time=552.00ms, real time=538.00ms, gc time=48.35ms

 

Download spc.mw

This is one of about ten ways I got it to work, depending on whether I used a fixed or dynamic working precision, whether I used unapply or explicit procedure bodies, and how I guarded against premature evaluation, etc. Eg, another that worked in 2018.2 and 2019.1,

restart;
r := t -> <cos(2*t), sin(3*t), 4*t>:
L := subs(dummy=unapply(sqrt(add(diff(r(t),t)^~2)),t),
          proc(tt) option remember; evalf[10](Int(dummy,0..tt)); end proc):
f := proc(s) option remember; fsolve(L-evalf(SFloat(s)), 0..infinity); end proc:
CodeTools:-Usage( plots:-spacecurve(r('f'(s)), s=0..30, color=red) );

It may be slightly inconvenient if you have to know the range of before constructing the tickmarks.

I have rough code to do it more automatically, from another project, but for another day.It involves two stages: 1) raw plot, then 2) analyze data in generated plot and redisplay with forced tickmarks. I'm sure that you could write code for that too.

restart;

plot(x*0.0005e13+3.986e13, x=0..5, axes=box,
     axis[2]=[tickmarks=[seq(i=sprintf("%a",round(i)),
                             i=3.986e13..3.9885e13,0.0005e13)]]);

 

Download yticks_integers.mw

If you plot your objective over that domain you can see that the surface is quite flat outside a narrow band near the the lower end of the y_1 range.

The GlobalSolve solver splits its search across the whole domain, and doesn't appear to look harder near the boundary than elsewhere. So with default options it may not look search enough to discover the larger objective values in that narrow band.

You can either restrict the domain or adjust options to make it search more thoroughly.

restart;

#Digits:=16;

p_l := 10^(-15):
epsilon := 1.09*10^(-10):
p_B := 1.09*10^(-8):

n_A := 10^7:
k_A := 0:
n_B := 10^8:
k_B := 0:

l := (x, y) -> x^k_A*(1 - x)^(n_A - k_A)*y^k_B*(1 - y)^(n_B - k_B);

proc (x, y) options operator, arrow; x^k_A*(1-x)^(n_A-k_A)*y^k_B*(1-y)^(n_B-k_B) end proc

GlobalOptimization:-GlobalSolve(l(x_1, y_1), x_1 = p_l .. epsilon, y_1 = p_B .. 1,
                                maximize, initialpoint = [x_1 = 0, y_1 = 0]);

[-0., [x_1 = HFloat(1.09e-10), y_1 = HFloat(0.6335488702113811)]]

GlobalOptimization:-GlobalSolve(l(x_1, y_1), x_1 = p_l .. epsilon, y_1 = p_B .. 1,
                                maximize);

[-0., [x_1 = HFloat(1.09e-10), y_1 = HFloat(0.6335488702113811)]]

# The surface is very flat for most of its y_1 range.
# By restricting the y_1 range it has a better chance
# of finding a global optimum for the default options
# that control "how hard it looks".
#
GlobalOptimization:-GlobalSolve(l(x_1, y_1), x_1 = p_l .. epsilon, y_1 = p_B .. 1e-4,
                                maximize);

[.335974651460662466, [x_1 = HFloat(7.190020594155661e-11), y_1 = HFloat(1.0900005612664663e-8)]]

# Alternatively, adjust the options that control
# "how hard it looks".
#
GlobalOptimization:-GlobalSolve(l(x_1, y_1), x_1 = p_l .. epsilon, y_1 = p_B .. 1,
                                maximize, populationsize=900, evaluationlimit=900*50);

[.336216487516351359, [x_1 = HFloat(1.1240574542239072e-15), y_1 = HFloat(1.09e-8)]]

# This really doing next to no work, and "getting lucky".
# It stalls out because it figures that it cannot improve
# from the supplied starting point (too flat).
# It's returning the best value (of a very few that it tried).
#
Optimization:-Maximize(l(x_1, y_1), x_1 = p_l .. epsilon, y_1 = p_B .. 1,
                       initialpoint = [x_1 = 0, y_1 = 0]);

Warning, no iterations performed as initial point satisfies first-order conditions

[.335850214825617777, [x_1 = HFloat(1.09e-10), y_1 = HFloat(1.09e-8)]]

# Same as above, using its own generated initial point.
#
Optimization:-Maximize(l(x_1, y_1), x_1 = p_l .. epsilon, y_1 = p_B .. 1);

Warning, no iterations performed as initial point satisfies first-order conditions

[-0., [x_1 = HFloat(1.09e-10), y_1 = HFloat(0.50000000545)]]

# For a restricted y_1 range it actually does more
# significant work, ascertaining the local optimum.
#
Optimization:-Maximize(l(x_1, y_1), x_1 = p_l .. epsilon, y_1 = p_B .. 1e-7);

[.335850214825617777, [x_1 = HFloat(1.09e-10), y_1 = HFloat(1.0899999999999998e-8)]]

# The plots illustrate the situation.
#
#plot3d(l(x_1, y_1), x_1 = p_l .. epsilon, y_1 = p_B .. 1);

#plot3d(l(x_1, y_1), x_1 = p_l .. epsilon, y_1 = p_B .. 1e-7);

 

Download got_example.mw

By setting infolevel[GlobalOptimization]:=6  you can see information about the values its using. The evaluationlimit is divided by the populationsize to determine the number of iterations-per-population. So when I increase the populationsize option I usually increase the evaluationlimit as well.

Some parts of the Maple product are enhanced (internationalized) by various language packs.

This and this page indicate that German is currently not one of the language packs.

And when I open the Maple 2019.1 GUI's preferences using Tools->Options->General from the menubar, German is not one of the choices for the Language drop-menu.

ps. Maplesoft's web-pages are partially available in French and German. So, slightly ironically, one can also read it here and here.

Reduce the time-cost of evaluation of your objective and constraints. It's not possible to give precise details because you have not bothered to provide a complete example.

Reduce the ranges, if possible. Reduce the dimension of the problem, if possible.

You could try something with the Explore command, to give you an idea of how the parameters affect the polar plot.

Eg,

Explore(plots:-polarplot(a+b*sin(x), x=0..T, title=evalf[3](a+b*sin(x)),
                         gridlines=true, coordinateview=[-3..3,0..2*Pi]),
        parameters=[a=-1.0..1.0, b=-2.0..2.0, T=0..2*Pi],
        initialvalues=[T=2*Pi]);

polarplot_explored.mw

Suppose that Maple is installed at a base location which is also the value of the Linux environment variable MAPLE. In that case the command-line interface can be launched on Linux in a terminal by issuing the command,

  $MAPLE/bin/maple

where that is the name of the bash script which launches the cmaple interface.

That script accepts several options, most of which can be listed by passing the -h option.

What hardware does your device have? I mean, which CPU chipset? Which OS distribution and version? Tell us how it goes.

 

There is no version of Maple which can read and write a proper video file format (eg. mpeg-2, etc).

 

Here is one way.  Note that in Maple 2018 the double-braces are not rendered around Unit calls.

restart

with(Units:-Simple)

with(InertForm)

m__fuel := 10*Unit('kg'/'s')

10*Units:-Unit(kg/s)

h__steam := 100*Unit('kJ'/'kg')

100*Units:-Unit(kJ/kg)

P__steam := `%*`(m__fuel, h__steam)

`%*`(10*Units:-Unit(kg/s), 100*Units:-Unit(kJ/kg))

Display(P__steam, inert = false)

0, "%1 is not a command in the %2 package", _Hold, Typesetting

Value(P__steam)

1000000*Units:-Unit(W)

 

Download inertform.mw

A list is not a mutable structure. When evaluation of a list entails a change in the evaluated value of an entry then a copy is made (and the address changes).

A Vector (or Matrix, Array, ie. rtable) is a mutable structure. In order for in-place semantics to work for rtables via procedure calls they are not copied under eval. Thus a procedure receives the very same rtable, when it is passed in as argument to the procedure call and evaluated up front. And so the desire for in-place semantics is one reason that rtables are not copied under eval.

Another reason is efficiency: the term "rtable" stands for "rectangular table", where hardware numeric entries (or the references to symbolic DAG entries) are stored contiguously in memory. It's a great benefit that passing the data of a float[8] rtable to an external routine only needs to pass the reference to its beginning. Full copying under normal evaluation rules for procedure calls would be a performance damper.

The functionality to force an evaluation throughout all entries of an rtable is provided by the rtable_eval command. Below you can see that command producing a copy.

restart;

g := r -> k*r;

proc (r) options operator, arrow; k*r end proc

g1 := g([x,y]);

k*[x, y]

g2 := g(<x,y>);

Vector(2, {(1) = k*x, (2) = k*y})

k := 1;

1

g1;

[x, y]

g2;

Vector(2, {(1) = k*x, (2) = k*y})

lprint(eval(g1,1));

k*[x, y]

lprint(eval(g1));

[x, y]

lprint(eval(g2,1));

Vector[column](2,{1 = k*x, 2 = k*y},datatype = anything,storage = rectangular,
order = Fortran_order,shape = [])

lprint(eval(g2));

Vector[column](2,{1 = k*x, 2 = k*y},datatype = anything,storage = rectangular,

order = Fortran_order,shape = [])

lprint(rtable_eval(g2));

Vector[column](2,{1 = x, 2 = y},datatype = anything,storage = rectangular,order
= Fortran_order,shape = [])

H1:=eval(g1);
H1[1]:=2;
g1;

[x, y]

[2, y]

[x, y]

H2:=eval(g2);
H2[1]:=55;
g2;

Vector(2, {(1) = k*x, (2) = k*y})

H2[1] := 55

Vector[column](%id = 18446884262579183854)

H3:=rtable_eval(g2);
H3[1]:=77;
g2;

Vector(2, {(1) = 55, (2) = y})

H3[1] := 77

Vector[column](%id = 18446884262579183854)

 

Download list_rtable.mw

note: When doing such examinations it can sometimes be useful to utilize lprint and eval(...,1), so as not to be misled by evaluation occuring only during printing.

Access to individual rtable entries, via indexing, does entail full evaluation. This is natural, for programmatic purpose.

restart;

V := <k^2, sin(k)>;

Vector(2, {(1) = k^2, (2) = sin(k)})

k := 1.2;

1.2

V;

Vector(2, {(1) = k^2, (2) = sin(k)})

lprint(eval(V));

Vector[column](2,{1 = k^2, 2 = sin(k)},datatype = anything,storage =

rectangular,order = Fortran_order,shape = [])

V[1];
V[2];

1.44

.9320390860

Download rtableentryeval.mw

While assignment to an entry of a list "works" (up to the 100th entry, by default), it is not true in-place semantics. The underlying list is actually replaced by a full copy whenever this occurs. This is thus a very inefficient programming practice, and should be avoided.

Here's a little more related to list evaluation,

restart;

f := proc(x) addressof(x); end proc:

L:=[1,2,k];

[1, 2, k]

addressof(L);

18446883949991687742

addressof(eval(L,1));

18446883949991687742

addressof(eval(L));

18446883949991687742

f(L);

18446883949991687742

k := 1;

1

addressof(L); # !!

18446883949981400014

addressof(eval(L,1));  # !!

18446883949991687742

addressof(eval(L));

18446883949981400014

f(L);

18446883949981400014

Download list_eval.mw

You could start by reading the Programming Guide, (in particular chapter 9, Object Oriented Programming) and the help pages related to objects and modules.

See also the help pages for Topics:
  object
  object,create
  object,methods

The most up-to-date version of the Programming Guide is included in the latest release of Maple. A (possibly older) version of chapter 9 is online here.

plot3d(piecewise(x^2+y^2<=1,0,undefined), y=-1..1, x=-1..1,
       lightmodel=none,
       image=cat(kernelopts(homedir),"/mapleprimes/","penny.png"));

For this example we can discern the greatest argument in a call to F, and isolate for that up front.

restart;

term := (r+1)*(r+2)*(r+3)*(r+4)*F(r+4)*(k-r+1)*F(k-r+1):

H := isolate(sum(term, r=0..k-1)
             +eval(term, r=k)=c, F(k+4));

F(k+4) = (c-(sum((r+1)*(r+2)*(r+3)*(r+4)*F(r+4)*(k-r+1)*F(k-r+1), r = 0 .. k-1)))/((k+1)*(k+2)*(k+3)*(k+4)*F(1))

Q := unapply(H, k):

Q(0);

F(4) = (1/24)*c/F(1)

Q(1);

F(5) = (1/120)*(c-48*F(4)*F(2))/F(1)

Q(2);

F(6) = (1/360)*(c-72*F(4)*F(3)-240*F(5)*F(2))/F(1)

 

Download gen.mw

First 9 10 11 12 13 14 15 Last Page 11 of 201