acer

32385 Reputation

29 Badges

19 years, 340 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Why are you trying to pass the plot function call around as if it were an object!? Do you have some special purpose where there really isn't a better programming methodology? (I can think of very few, and since you are new to Maple I suspect that your actual goal could be better done in another way altogether).

Having said that,

  eval('plot'(x, x = 0 .. a), a = 3);
  value(eval(%plot(x, x = 0 .. a), a = 3));
  ((e::uneval)->eval(e,a=3))( plot(x, x = 0 .. a) );

On this forum it is almost always better to inform us in detail of your actual motivating goal, rather than just inquiring about some barrier found during an attempt you devised for getting there.
 

You can customize and adjust this.

The idea is that you can write a convenient procedure to generate such plots in a reusable manner. And you could also augment the results futher.

restart;

F:=proc(ee,var,a,b,n::posint,
        {lineoptions::list:=[],
         ptoptions::list:=[]})
  local i,pts,xpts;
  uses plots;
  xpts:=[seq(a+(i-1)*(b-a)/(n-1),i=1..n)];
  try
    pts:=map(p->[p,eval(ee,var=p)],xpts);
  catch:
    pts:=map(p->[p,limit(ee,var=p)],xpts);
  end try;
  display(plot(ee,var=a..b,_rest),
          seq(plottools:-line([p[1],0],p,lineoptions[]),p=pts),
          pointplot(pts,ptoptions[]),
          ':-axis[1]'=[':-tickmarks'
                       =(u->u=Typesetting:-Typeset(u))~(xpts)]);
end proc:

F(sin(x),x,-2*Pi,2*Pi,7);

F(cos(x),x,-2*Pi,2*Pi,9,color=green,
  lineoptions=[color=blue,linestyle=dash],
  ptoptions=[color=red,symbol=solidcircle,symbolsize=10],
  size=[400,300]);

F(sin(x)/x,x,-2*Pi,2*Pi,9,color=blue,thickness=3,
  lineoptions=[color=blue,thickness=3],
  ptoptions=[color=blue,symbol=solidcircle,symbolsize=20],
  size=[500,200], labels=["",""]);

plots:-display(
  F(sin(x)/x,x,-2*Pi,2*Pi,13,color=blue,thickness=3,
    lineoptions=[color=blue,thickness=3],
    ptoptions=[color=blue,symbol=solidcircle,symbolsize=20],
    size=[600,300], labels=["",""], axes=box),
  plot(0,-2*Pi..2*Pi,color=black),
  plottools:-line([0,1],[0,-0.3],color=black),
  axis[2]=[tickmarks=[0,0.5,1.0]]
);

 

Download plotting_fun.mw

restart;

with(Units:-Standard):

a := 153350 * Unit(J);

153350*Units:-Unit(J)

#
# This doesn't work. (It's a mistaken attempt.)
#
convert(a,'units',Unit('kN')*Unit('m'));

(3067/20)*Units:-Unit(1000*J)

#
# This works
#
convert(a,'units',Unit('kN'*'m'));

(3067/20)*Units:-Unit(kN*m)

#
# The first attempt is a mistake, because the
# 3rd argument is itself combined, unitwise.
#
Unit('kN')*Unit('m');

1000*Units:-Unit(J)

#
# Alternatively, you can just use the synbols.
#
convert(a,'units','kN'*'m');

(3067/20)*Units:-Unit(kN*m)

It can also be done in 2D Input mode.

convert(a, 'units', 'kN'*'m')

(3067/20)*Units:-Unit(kN*m)

It can also be done sing the Units palette, by choosing
the blue "units" item. Hence the units appear in upright
Roman in Maple 2020 (or braced, when inlined on Mapleprimes).

convert(a, 'units', Unit('kN'*'m'))

(3067/20)*Units:-Unit(kN*m)

 

Download unitconvert.mw

Since NULL is ephemeral when used in ways that you intend (eg. being passed around, or returned, in an expression sequence) then it isn't a very good choice. Don't choose something that can be problematic for one of your basic usage patterns.

Using {} is not terrible, but it will fall down for a computation where {} is also a possible natural result and where you need to distinguish between the two situations. In that case you'd have to utilize something else instead, which means that your methodologies could be inconsistent. Note that you didn't cite only computations where the result denoted possible solutions -- so {} might make sense to denote the empty set, but is not great to denote general cases where something "could not be computed".

How about using FAIL instead? Note that your two cited example circumstances were, "...it could be something that could not be computed inside foo, or some solution that could not be found". The common phrasing there is "could not be".

For the cases of computing possible solutions, would you want to be able to distinguish between the case of provably no solutions and the case where it happens to not succeed in finding any solution?

[edit] I don't think that solutions involving returning uneval-quoted NULL is great. You'd need to avoid using whatever it got assigned to in an analogous situation. Why introduce any need for extra care, if you don't have to?

Here is one way to deal with that, which consists of the step,

   subsindets(eq[4], float,u->`if`(frac(u)=0.0,trunc(u), u));

(There are more complicated ways, some of which deal with other things that you might encounter and want, such as stripping off trailing zeroes. Let us know...)

restart

with(LinearAlgebra):

par := {a = 2.5, alpha = 2, k_r = .5, k_rc = .2, k_rq = .2, kappa = 0.1e-2, mu_k = .35, mu_s = .44, omega = .766620580157922, sigma_0 = 110, sigma_1 = 1.37, sigma_2 = 0.823e-1, x_s3 = -1.04003626422324936017819852700633621040584050594846801927800, zeta = 0.904504977553123318334601762827181333680096702957228781770315e-1}:

f := proc (v_r) options operator, arrow; mu_k+(mu_s-mu_k)*exp(-a*v_r) end proc:

``

g_exp1 := taylor(1/f(v_rv+x21), x21 = 0, 4):

for k to 4 do g_coeff[k] := taylor(subs(v_r = v_rv, coeff(g_exp1, x21, k-1)), epsilon = 0, 3) end do:

g0 := eval(subs(par, coeff(g_coeff[1], epsilon, 0))):

g1 := eval(subs(par, coeff(g_coeff[2], epsilon, 0))):

g2 := eval(subs(par, coeff(g_coeff[3], epsilon, 0))):

g3 := eval(subs(par, coeff(g_coeff[4], epsilon, 0))):

``

NULL

eq[1] := subs(par, diff(x[1](t), t)-x[2](t)):

eq[2] := subs(par, diff(x[2](t), t)+2*zeta*x[2](t)+x[1](t)+k_r*(x[1](t)-x[3](t))+2*kappa*(x[2](t)-x[4](t))+k_rq*(x[1](t)-x[3](t))^2+k_rc*(x[1](t)-x[3](t))^3)

diff(x[2](t), t)+.1829009955*x[2](t)+1.5*x[1](t)-.5*x[3](t)-0.2e-2*x[4](t)+.2*(x[1](t)-x[3](t))^2+.2*(x[1](t)-x[3](t))^3

eq[3] := subs(par, diff(x[3](t), t)-x[4](t))

diff(x[3](t), t)-x[4](t)

eq[4] := subs(par, diff(x[4](t), t)+2*kappa*alpha*(x[4](t)-x[2](t))+k_r*alpha*(x[3](t)-x[1](t))-k_rq*alpha*(x[3](t)-x[1](t))^2+k_rc*alpha*(x[3](t)-x[1](t))^3+alpha*(sigma_0*x[5](t)+sigma_1*v_r*(1-sigma_0*x[5](t)*(g_0+g_1*x[4](t)+g_2*x[4](t)^2+g_3*x[4](t)^3))+sigma_2*v_r)):

diff(x[4](t), t)+0.4e-2*x[4](t)-0.4e-2*x[2](t)+x[3](t)-x[1](t)-.4*(x[3](t)-x[1](t))^2+.4*(x[3](t)-x[1](t))^3+220*x[5](t)+2.74*v_r*(1-110*x[5](t)*(g_0+g_1*x[4](t)+g_2*x[4](t)^2+g_3*x[4](t)^3))+.1646*v_r

eq[5] := subs(par, diff(x[5](t), t)-v_r*(1-sigma_0*x[5](t)*(g_0+g_1*x[4](t)+g_2*x[4](t)^2+g_3*x[4](t)^3)))

diff(x[5](t), t)-v_r*(1-110*x[5](t)*(g_0+g_1*x[4](t)+g_2*x[4](t)^2+g_3*x[4](t)^3))

NULL

Download question_sunit_ac1.mw

Here is one way.

Let me know if you really prefer that multiple calls foo(ode,y,x) would produce the exact same result (in which case it could adjusted to be store and manage that, I think.)

restart;

foo := proc(ode::`=`,y,x::symbol)
  local t:=`tools/genglobal`(x);
  if t=x then t:=`tools/genglobal`(x); end if;
  return PDEtools:-dchange({x=t^2},ode,t)  
end proc:

ode:=diff(y(x),x)=sin(x);
foo(ode,y,x);
foo(ode,y,x);

diff(y(x), x) = sin(x)

(1/2)*(diff(y(x0), x0))/x0 = sin(x0^2)

(1/2)*(diff(y(x1), x1))/x1 = sin(x1^2)

ode:=diff(y(t),t)=sin(t);
foo(ode,y,t)

diff(y(t), t) = sin(t)

(1/2)*(diff(y(t0), t0))/t0 = sin(t0^2)

 

Download genglobal.mw

Be careful about terminology. You wrote, "I think, `evalb` is used to compare numbers and `is` used to compare expression, logic. Is this true?" That sentence is not precise, in Maple terminology, and that can lead to misunderstanding here.

Think instead in terms of technical Maple terminology, for example whether something is of type numeric. That includes floating-point numbers, integers, and fractions (rationals, with integers in numerator and denominator). That does not include exact representations such as sqrt(2). Yes, that denotes a "number" in the common mathematical sense, but it is not of type numeric to Maple. A related and useful class is things of type realcons, which would produce something of type numeric if evalf were applied.

The evalb command can compare equality and inequality relations of items that are both of type numeric. It can also test whether two structues are identically the same, by which I mean that they have the same address in memory.

So you cannot generally, reliably compare two expressions (of type algebraic but not numeric) with the evalb command unless you can put both into a canonical form, in which case they are uniquified and the canonical forms are identically the same.

Your example,
  f := x + 1 = (x^2 - 1)/(x - 1);
can be simplified (using simplify or normal, for example) so that lhs and rhs of that normal form are in fact the very same expression. That is why evalb(simplify(f)) returns true.

But not all expressions can be easily put by Maple into some canonical form. (For example, there are mathematically equivalent formulations for some expressions involving trig and special function calls for which a choice of canonical reformulation is not obvious.) That's one advantage of using is, as it has a variety of mathematical operations it can do for comparison.

Another significant advantage of the is command is that it can (often) utilize assumptions on names in symbolic expressions.

Here are some tips:

1) For mathematical testing of an equation (ie. A=B), it can sometimes help is by querying is(A-B=0) since (while zero-testing is hard) at least the target form is clear: zero.

2) In general the is command performs more strongly that simplify (or its friends, normal, radnormal, shake, evalc, rationalize, combine, and expand) alone. It is rare that it helps to make a call such as is(simplify(A-B)=0) , but examples do exist.

3) For inequality testing involving expressions of type realcons you may generally be better off invoking the is command, rather than applying evalf and risking a mistake due to close roundoff error. Like most "rules", this too has its exceptions.

4) If you are going to use the is comand then you might have to write your code so that it handles the case that is returns FAIL instead of true or false.

5) For verification of floating-point results from a computation against a target value you might look at the testfloat command, or the verify command with a float comparator.

6) For verification of data structures (lists, sets, Vectors, etc) you could also look at a call of the verify command, according to the type of the structures' members.

It's a large topic, and a textbook on Maple could be written while covering it comprehensively.

Background reading, Help pages for Topics:
   -  testfloat
   - verify
   - evalrshake
   - testeq
   - signum
   - type numeric, realcons, algebraic
   - evalc, radnormal, evala, normal

The catch clause is intended to catch whatever throws a catchable error with the given prefix1, and there could be more than one statement in the try clause that emits such.

One point of the finally clause is that its code is supposed to run if either an error is caught or no error occurs. If you don't want some code run when a error is caught then put it later in the try clause, instead -- eg. after your conditional.

In complicated scenarios you may find it helpful to set and use one or more flags, or to use a nested try..catch, etc.

Does the error message upon time-out always begin with "time expired"?

1 Managing uncatchable errors is a different story.

You can write your own wrapper around ColumnGraph, allowing you to retain your plots:-setoptions setting in your initialization file. You can make it ignore or respect the setting, and override it.

(A bug report will be submitted...)

restart;

kernelopts(version);

`Maple 2020.1, X86 64 LINUX, Jun 8 2020, Build ID 1474533`

plots:-setoptions(size=[250,250]);

with(Statistics):

#
# This version ignores setoptions choice for size,
# but allows you to force it.
#
SColumnGraph:=proc()
  uses plots; local oldsize, P;
  oldsize:=setoptions(':-size'):
  setoptions(':-size'=[':-default',':-default']);
  try
    P:=Statistics:-ColumnGraph(_passed);
  catch "":
    error;
  finally
    setoptions(':-size'=oldsize);
  end try;
  return P;
end proc:

T := [StringTools[CharacterFrequencies]("antidisestablishmentarianism")]:

SColumnGraph(T);

SColumnGraph(T, 'size'=plots:-setoptions(':-size'));

SColumnGraph(T, 'size'=[300,150]);

#
# This version respects setoptions choice for size,
# but allows you to disregard it.
#
RColumnGraph:=proc()
  uses plots; local oldsize, P, S;
  oldsize:=setoptions(':-size'):
  setoptions(':-size'=[':-default',':-default']);
  try
    S:=select(type,[_passed],identical(size)=anything);
    S:=`if`(nops(S)>0,S[-1],':-size'=oldsize);
    P:=Statistics:-ColumnGraph(_passed,S);
  catch "":
    error;
  finally
    setoptions(':-size'=oldsize);
  end try;
  return P;
end proc:

T := [StringTools[CharacterFrequencies]("antidisestablishmentarianism")]:

RColumnGraph(T);

RColumnGraph(T, 'size'=[':-default',':-default']);

RColumnGraph(T, 'size'=[300,150]);

 

Download ColumnGraph_Problem_ac.mw

Could you not also try using a personal Maple initialization file, so that you wouldn't have to change its contents when you upgrade your Maple release?

That is to say, from the initialization Help page,
    "C:\Users\userid\maple.ini"
where userid is your Windows id. (See also kernelopts(':-homedir') output, if relevant.)

Let me know if I did this incorrectly. This is not thread-safe (but I suspect that you could make a variant that is -- each thread getting its own parameter name and dummy container Vector).

[edit] I see that mmcdara has written about subscribers being a discrete uniform random variable over {0, ..., 16}, and getting the fast direct generation of the sample. But I won't delete this, since it might interest someone later, to compare varying the parameter value of a preformed RV against generating a sample from a fresh call to Binomial (or other distribution flavour) for each entry.

restart:

randomize():

with(Statistics):

n_draw      := 10000:
prior_rate  := Sample(Uniform(0, 1), n_draw):
n_trials    := 16:
dummyV      := Vector(1,datatype=float[8]):
RV          := RandomVariable(Binomial(n_trials,pat)):

subscribers := CodeTools:-Usage( map(proc(u) option hfloat;
                                       global pat;
                                       pat := u;
                                       Sample(RV,dummyV);
                                       dummyV[1];
                                     end proc,
                                     prior_rate) ):

memory used=110.80MiB, alloc change=33.00MiB, cpu time=753.00ms, real time=754.00ms, gc time=48.49ms

Histogram(subscribers, discrete=true, view=[-1..n_trials+1, default], thickness=10);

 

Download subscribers_ac1.mw

You could compute it using the geom3d package, either by command or from coordinates (after a projection). Or you can compute it with LinearAlgebra. See below.

If you wanted you could add labels to the points (vertices, midpoint, and projection to the plane).

restart;

with(geom3d):

point(A,1,1,1), point(B,1,-1,-1),
point(C,-1,1,-1), point(F,-1,-1,1):

gtetrahedron(T,[A,B,C,F]):

plane(P1,[A,B,C]), plane(P2,[A,B,F]):

FindAngle(P1,P2);
evalf(%);
evalf(%*180/Pi);

arccos(1/3)

1.230959417

70.52877935

projection(Q, F, P1):  # Q is the projection of F onto plane P1
midpoint(M, A, B):     # M is the midpoint between A and B

triangle(Tr1,[A,B,C]), triangle(Tr2,[A,B,F]):
segment(L_FQ,[F,Q]), segment(L_MQ,[M,Q]), segment(L_MF,[M,F]):

plots:-display(draw(Tr1,color="Orange"),draw(Tr2,color="Orange"),
               draw(T,transparency=0.9,color=gray),
               draw(Q,symbolsize=20,symbol=solidsphere),
               draw(F,symbolsize=20,symbol=solidsphere),
               draw(M,symbolsize=20,symbol=solidsphere),
               draw(L_FQ,linestyle=dash,thickness=3),
               draw(L_MQ,linestyle=dot,thickness=3),
               draw(L_MF,linestyle=dot,thickness=3),
               orientation=[5,-15,-20], size=[500,500],
               view=[-1..1,-1..1,-1..1], axes=box);

# Alternatively, compute it from the coordinates
# of F, Q, and M

cQ:=coordinates(Q);
cF:=coordinates(F);
cM:=coordinates(M);

[1/3, 1/3, -1/3]

[-1, -1, 1]

[1, 0, 0]

with(LinearAlgebra):

arccos( <cM-cF>^%T . <cM-cQ> / (Norm(<cM-cF>,2)*Norm(<cM-cQ>,2)) );

arccos(1/3)

# Or, without geom3d

vA:=<1,1,1>: vB:=<1,-1,-1>:
vC:=<-1,1,-1>: vF:=<-1,-1,1>:
X1:=CrossProduct(vA-vC, vA-vB):
X2:=CrossProduct(vA-vF, vA-vB):
arccos( X1 . X2 / (Norm(X1,2)*Norm(X2,2)) );

arccos(1/3)

 

Download regtetra_dihedral.mw

Firstly, using a common bracketed notation for pretty-printing the result, and then some alternatives.

restart;

MySumP := (j::integer) ->
  seq(%binomial((j+10),(i+10))*p^(j+i)*(1-p)^(j+i), i= 0..0):

foo := MySumP(10);

%binomial(20, 10)*p^10*(1-p)^10

foo;

%binomial(20, 10)*p^10*(1-p)^10

value(foo);

184756*p^10*(1-p)^10

# If you do not prefer the gray round brackets

InertForm:-Display(foo, ':-inert'=false);

%binomial(20, 10)*p^10*(1-p)^10

# optionally

`print/%binomial`:=proc() ':-C'(args); end proc:

 

foo;

%binomial(20, 10)*p^10*(1-p)^10

value(foo);

184756*p^10*(1-p)^10

`print/%binomial`:=proc(m,n) ':-C'[n]^m; end proc:

 

foo;

%binomial(20, 10)*p^10*(1-p)^10

value(foo);

184756*p^10*(1-p)^10

# Now with upright Roman "C"

`print/%binomial`:=proc(m,n)
    uses Typesetting;
    msubsup(mtext("C"),mn(convert(n,string)),
            mn(convert(m,string)));
end proc:

 

foo;

%binomial(20, 10)*p^10*(1-p)^10

value(foo);

184756*p^10*(1-p)^10

`print/%binomial`:=proc(m,n)
    uses Typesetting;
    mscripts(mi("C"),mn(convert(n,string)),
             none()$3,mn(convert(m,string)),none());
end proc:

 

foo;

%binomial(20, 10)*p^10*(1-p)^10

value(foo);

184756*p^10*(1-p)^10

# Now with upright Roman "C", and smaller numerals

`print/%binomial`:=proc(m,n)
    uses Typesetting;
    mscripts(mtext("C"),mn(convert(n,string),':-mathsize'=10),
             none()$3,mn(convert(m,string),':-mathsize'=10),none());
end proc:

 

foo;

%binomial(20, 10)*p^10*(1-p)^10

value(foo);

184756*p^10*(1-p)^10

NULL

Download inertbinomial.mw

If you don't want your PDF file interspersed with empty spaces due to forced page-breaks then you might also try using Array-plots for your 3D plots.

Eg,
   plots:-display(Array([plot3d(sin(x))]));
instead of,
   plot3d(sin(x));

That will put a Table border around each such plots, though you can toggle off the visibility of its borders using right-click menu choice Table->Properties just before exporting.

A:=[1,2,3]:
B:=[7,8,9]:

zip(f,A,B);

            [f(1, 7), f(2, 8), f(3, 9)]

First 114 115 116 117 118 119 120 Last Page 116 of 336