acer

32333 Reputation

29 Badges

19 years, 318 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

The result coming out of symbolic int (from the meijerg method) is incorrect.

You can very quickly produce a plot as follows, however, using numeric integration. Note the use of inert (captalized) Int below, instead of active int.

Here I get good speed by forcing a numeric integration method which avoids (here unnecessary) expensive discontinuity checking.

plot(Int(exp(-(x^2-a)^2),x=0..infinity,method=_d01amc),
     a=-2..5);

Other quick approaches here include specifying a decent method other than _d01amc (eg. _Dexp) or utilizing operator form for the integrand within the Int call (slightly tricky since the value of parameter a has to be shoehorned in). These are not quite as quick as the approach given above.

plot(Int(exp(-(x^2-a)^2),x=0..infinity,method=_Dexp),
     a=-2..5);

plot(a->Int(unapply(exp(-(x^2-a)^2),x),0..infinity),
     -2..5);

Using Int (or int with its numeric option) alone for this example is considerably slower slower than the approaches mentioned above. These next approaches will also produce the better plot, but they are not nearly as fast (because of hybrid symbolic-numeric checking), taking about 60sec versus less than 1sec for the first approach above. The slowness is not due to the accuracy tolerance.

plot(Int(exp(-(x^2-a)^2),x=0..infinity),a=-2..5);

plot(a->int(exp(-(x^2-a)^2),x=0..infinity,numeric),-2..5);

I will submit a bug report against the symbolic int command, for the incorrect result coming out of the meijerg method.

If you read the Help page for Topic proc you can see a description of the fact that local variables only get 1-level evaluation.

You code does a rather unusual and roundabout trick of creating expressions with references to y, and then temporarily assigning specific values to y (eg. y=0, y=b), and then unassigning y. It does this while constructing the expressions that get assigned to EQ1,EQ2,EQ3, and EQ4.

When you do all that at the top-level then the expressions assigned to EQ1,etc get fully evaluated, and the references to y are replaced by the values temporarily assigned to y, and no references to y remain.

But when you do that within a procedure (where y, GenFunc, EQHelp, and EQHelp2 are local variables -- whether declared explicitly or not) the references to local y only get 1-level evaluation, and are not resolved to the temporary assigned values.

There are several ways around this:

A) You could wrap in a simple call to eval, each of the four constructions. For example, instead of,
          EQ1 := GenFunc = 0;
you could have,
          EQ1 := eval(GenFunc = 0);
while retaining you temp assignments to y, and unassignments of y (so that the diff calls were still valid). This is workable, but retains what I consider to be a very awkward way of doing something that should be simple. The fact that you have to repeatedly unassign y is a hint that is suboptimal. This is not the best direction in which to move. Test_GetDet_MaplePrimes_assign_eval_unassign.mw

B) You could declare EQHelp, EQHelp2, and GenFunc as a global at the top of the wrapping GetDet procedure, and keep all the other mechanisms as they are. That works, but retains all the awkwardness. It's actually more awkward and convoluted. (For other examples you may get lost trying to figure out what to declare global. Also, unnecessary global declarations are evil and wicked.) This is the wrong direction in which to move. Test_GetDet_MaplePrimes_globaldecl.mw

C) You could get rid of all the need for repeated assigning and unassigning to y, by instead using the so-called 2-argument calling sequence of the eval command to evaluate GenFunc, EQHelp, and EQHelp2 at the specific values of y. This is much more usual and convenient way to utilize a value for a variable without having to unassign it for subsequence symbolic use (the ensuing diff calls). This is the right way to accomplish this part of the task (IMO), regardless of whether inside a wrapping procedure or all done at the top-level. Test_GetDet_MaplePrimes_ac.mw

You keep submitting duplicates of earlier question threads. You spam this site with them. It is inconsiderate and unhelpful. I flag them as duplicates, and then delete them. I have given you messages about this -- as private messages and comments on earlier Questions.

Put your followup queries onto one of your related earlier Questions, instead.

If you have other kinds of example then you should show them.

restart;

expr := 4.22647099771398*exp(-0.0104163888375266*t)
        - 4.22647099771398*exp(-2.00033505979518*t);

4.22647099771398*exp(-0.104163888375266e-1*t)-4.22647099771398*exp(-2.00033505979518*t)

indets(expr,specfunc(anything,exp));

{exp(-2.00033505979518*t), exp(-0.104163888375266e-1*t)}

# note the order of the terms (we don't know
# whether that matters to you.
map(op,indets(expr,specfunc(anything,exp)));

{-2.00033505979518*t, -0.104163888375266e-1*t}

# This preserves the "ordering" of the terms in the original sum.
[seq([op(map(op,indets(term,specfunc(anything,exp))))], term=expr)];

[[-0.104163888375266e-1*t], [-2.00033505979518*t]]

# This preserves the "ordering" of the terms in the original sum.
# Here the operands of the exp calls are not separated, so
# only use this if you know that each term has only one such
# call.
[seq(op(map(op,indets(term,specfunc(anything,exp)))), term=expr)];

[-0.104163888375266e-1*t, -2.00033505979518*t]

 

Download indets_exp.mw

Some instructor has a sense of humour.

(I made no effort for efficiency.)

Guessing at how num2text might work was not difficult.

restart;

num2text := proc(n::posint) local C, i;
  C := convert(n,'base',100);
  cat(seq("abcdefghijklmnopqrstuvwxyz"[C[nops(C)-i+1]],
          i=1..nops(C)));
end proc:

eqs := {z =  4*x + 19*y + 2515211725275120,
        z = 52*x + 27*y + 2515210613496048,
        z = 36*x + 65*y + 2515210981587340,
        z =  6*x + 60*y + 2515211676449260}:

Q := eval([x,y,z], solve(eqs,{x,y,z}));

[23151811, 61518, 2515211819051206]

map(num2text, Q)[];

"work", "for", "yourself"

cat(num2text(Q[1])," ",num2text(Q[2])," ",num2text(Q[3]));
 

"work for yourself"

#cat(num2text(Q[1]), seq([" ",num2text(Q[i])][],i=2..nops(Q)));

 

Download hwQuestion.mw

If you have a set of the form,

   {var1=..., lambda22=..., var2=..., var3=...}

then you can evaluate any formula at that "point" (set of equations of the form name=value).

In particular you can evaluate the simple formula of just lambda22 at that point.

This is especially useful when the equations are in a set, because the mechanism works regardless of the order of the equations in the set.

restart;

with(LinearAlgebra):

kappa:=30:
tensF:=Matrix([[lambda11,0,0],[0,lambda22,0],[0,0,lambda22]]):
tensFbar:=Matrix([[lambda11bar,0,0],[0,1/sqrt(lambda11bar),0],[0,0,1/sqrt(lambda11bar)]]):
tensBbar:=Multiply(tensFbar,Transpose(tensFbar)):

tensPK:=Matrix([[pi,0,0],[0,0,0],[0,0,0]]):
tensPK1:=Multiply(Transpose(MatrixInverse(tensF)),(2/Determinant(tensF))*a*(tensBbar-(1/3)*Trace(tensBbar)*Matrix(3,shape=identity))+kappa*(Determinant(tensF)-1)*Matrix(3,shape=identity)):
eq1:=tensPK[1,1]=tensPK1[1,1]:
eq2:=tensPK[2,2]=tensPK1[2,2]:
sys:={eq1,eq2}:

sol1:=[solve( sys, {lambda11bar,lambda22},explicit)]:
nops(sol1);

6

eval(lambda22,sol1[1]);

(1/30)*10^(1/2)*(lambda11*(pi*lambda11+90))^(1/2)/lambda11

eval(lambda22,sol1[2]);

-(1/30)*10^(1/2)*(lambda11*(pi*lambda11+90))^(1/2)/lambda11

seq([i,eval(lambda22,sol1[i])], i=1..nops(sol1));

[1, (1/30)*10^(1/2)*(lambda11*(pi*lambda11+90))^(1/2)/lambda11], [2, -(1/30)*10^(1/2)*(lambda11*(pi*lambda11+90))^(1/2)/lambda11], [3, (1/30)*10^(1/2)*(lambda11*(pi*lambda11+90))^(1/2)/lambda11], [4, -(1/30)*10^(1/2)*(lambda11*(pi*lambda11+90))^(1/2)/lambda11], [5, (1/30)*10^(1/2)*(lambda11*(pi*lambda11+90))^(1/2)/lambda11], [6, -(1/30)*10^(1/2)*(lambda11*(pi*lambda11+90))^(1/2)/lambda11]

 

A poor alternative is to inspect the set for the location of the one you want, and to pick off the value by using rhs of (here the second entry) entry indexed by position. It is a weak methodology because the set order depends on lexicographic ordering of the set (eg. alphapbetical ordering of the lhs names). If the names were changed the position of the one you wanted could change. So the code would be less robust.

 

rhs(sol1[1][2]);

(1/30)*10^(1/2)*(lambda11*(pi*lambda11+90))^(1/2)/lambda11

rhs(sol1[2][2]);

-(1/30)*10^(1/2)*(lambda11*(pi*lambda11+90))^(1/2)/lambda11

 

Another alternative is to specify the order of the solved variables, in the call to solve, by supplying the solved variables in a list rather than in a set. The effect is that the results will also contain lists, with the same ordering.

 

You can then pick off the values acording to the position of the one you want (in the list of variable names passed to solve).

 

But you still have to remember to re-use the same ordinal position, when picking off the values. Also, the purpose of the lines with rhs,etc may be less obvious to you or someone else at a later date.

 

Hence I prefer the methodology using eval, which involves less bookkeeping and works for both lists and sets of equations name=value.

 

# Notice the second argument passed to `solve` is now a list.
# Also, with this calling sequence `solve` returns a list of
# the solutions, so I didn't wrap the `solve` call in [] brackets.

sol1:=solve( sys, [lambda22,lambda11bar],explicit):

rhs(sol1[1][1]);

(1/30)*10^(1/2)*(lambda11*(pi*lambda11+90))^(1/2)/lambda11

rhs(sol1[2][1]);

-(1/30)*10^(1/2)*(lambda11*(pi*lambda11+90))^(1/2)/lambda11

 

Download solve_eval.mw

There are several ways to accomplish this, so I'll give one of several alternatives.

One of the reasons to split off the "first argument" passed to plots:-animate is to provide a (sometimes easier) mechanism that having to use unevaluation-quotes. The original authors of plots:-animate recognized that avoiding premature evaluation would very often be wanted -- because the animating parameter would not yet have a numeric value -- and hence allowed a syntax that split the command from its arguments.

For your example, that same mechanism can handle it directly, as opposed to also utilizing unevaluation-quotes.

restart;
with(DynamicSystems):
with(plots):

num:=10*(s+alpha);
den:=s*(s^2+4*s+8);

animate(RootLocusPlot@NewSystem,
        [subs(alpha=a,num/den)],
        a=1..10);

f:=unapply(num/den, alpha):
animate(RootLocusPlot@NewSystem,[f(a)],
        a=1..10);

Of course there will be other examples where other flavours are easier.

A subtlety is that you do want f(a) or the subs call to evaluated up-front. It'd be (here, slightly) less efficient to have that be done upon each new value taken by the parameter a. The solutions discussed so far in this thread all do that part well.

The first step is to ensure that your procedure(s) can run under hardware double-precision mode, the simplest mechanism for which is by using the evalhf interpreter.

The plot3d command will attempt to use evalhf mode, but we need to make the TF procedure "evalhf'able". The following takes 4 seconds on my machine, while the original took 412 seconds.

restart

kernelopts(version);

`Maple 18.02, X86 64 LINUX, Oct 20 2014, Build ID 991181`

eps := 0.1e-1;

0.1e-1

N := 20;

20

maxIter := 100;

100

bail := 0.1e11;

0.1e11

orb := Array(1 .. N, datatype = complex[8]):

TF := proc (x, y, orb, N, eps) local n, p, z, c, lgc; z := 0; c := x+I*y; lgc := ln(c); for n to N do orb[n] := 0 end do; for n to maxIter while abs(z) < bail do z := exp(lgc*z); if maxIter-N < n then orb[n-maxIter+N] := z end if end do; p := 0; for n to N-1 do if abs(orb[N]-orb[N-n]) <= eps then p := n; break end if end do; log(p) end proc:

evalhf(TF(2.1, .3, orb, N, eps));

1.60943791243410028

p := CodeTools:-Usage(plot3d(('TF')(a, b, orb, N, eps), a = -4 .. 4, b = -4 .. 4, grid = [400, 400], lightmodel = none, shading = ZHUE, style = PATCHNOGRID, orientation = [-90, 0], axes = BOX, labels = [x, y, iter])):

memory used=3.12MiB, alloc change=1.22MiB, cpu time=3.96s, real time=3.96s, gc time=0ns

 

NULL

Download tetrationFractal18_ac.mw

The next step is to ditch using plot3d as a mechanism for generating the data values, since we can avoid having to invoke evalhf being called for every x-y point. Instead we could write a single evalhf'able procedure that operates in-place on a predefined 2-dimensional Array (or Matrix).

While re-writing in that way, the coloring scheme could be revamped.

Further improvement would be to make that paralellizable (each thread operating on only a subrange of the x-values, say, but writing to the shared result Matrices), and compile'able under Compiler:-Compile.

note: In more recent Maple versions tetration can be attainted using the IterativeMaps:-Escape command, although that is with a fixed iteration limit (and not so easily with the fancy orbit check). I already had an example of that, but it won't run in Maple 18 since the package was introduced in Maple 2015.

After you assign a value to A the thing assigned to V still contains A.

If you then access individual elements of  V then evaluation occurs, and you'd get what you expect.

If you subsequently assign a different value to A the you could get a different result when accessing the relevant entries of V. (You might want to do this.)

But that evaluation isn't occurring when you print V.

You can map the eval command over V, or apply the rtable_eval command to V, to see the printing display all the elements after evaluation.  Eg,

     map(eval, V);

     rtable_eval(V);

Those both create a new instance of the thing, that contains the value you assigned to A rather than the actual name A. (If you subsequently assign a different value to A then these new instances aren't affected.)

The above describes "how" the behavior happens. Do you also want the (involved) explanation of the design decisions behind rtable evaluation rules?

Here are illustrations of the explanation above (of "how").

restart

with(VectorCalculus)

V := `<,>`(A, B, C)

Vector(3, {(1) = A, (2) = B, (3) = C})

A := 1

1

V

Vector(3, {(1) = A, (2) = B, (3) = C})


We can see what is actually stored in V.

lprint(V)

Vector[column](3,{1 = A, 2 = B, 3 = C},datatype = anything,storage =
rectangular,order = Fortran_order,attributes = [coords = cartesian],shape = [])

 

Accessing the individual elements (entries) gets the evaluated result.

 

V[1]

1


We can force the evaluation of elements (also making a new structure).

W := map(eval, V)

Vector(3, {(1) = 1, (2) = B, (3) = C})

Y := rtable_eval(V)

Vector(3, {(1) = 1, (2) = B, (3) = C})

 

If we now assign a different value to A then we can observe that
doing so affects the evaluated elements of V, but not those of Y (or W).

 

A := 17

17

rtable_eval(V)

Vector(3, {(1) = 17, (2) = B, (3) = C})

rtable_eval(Y)

Vector(3, {(1) = 1, (2) = B, (3) = C})

lprint(Y)

Vector[column](3,{1 = 1, 2 = B, 3 = C},datatype = anything,storage =
rectangular,order = Fortran_order,attributes = [coords = cartesian],shape = [])

 

Download VectorCalculusQuestion_ac.mw

Yes, the parameter names generated by LinearSolve are indexed, and can be found by their type.

You can also programatically generate the base name, as a "new" (previously unassigned) name.

But I don't like your methodology of finding them, as it relies on each such name being present as an entry value (alone) of sol. That might be true as consequence of the LinearSolve algorithm, but it's an unnecessary weakness. And there are also unnecessary mappings and conversions.

You can choose the base name, and pass that to LinearSolve using its free option. And you can even programmatically generate such base names -- unassigned. That's pretty decent, and safe, control.

This alleviates the difficulty you've mentioned with not knowing what name like _t, _t2, etc, that will be indexed in the result. If you don't generate the base name using `tools/genglobal` then LinearSolve will. That's the same mechanism as LinearSolve uses.

So, generate the base name yourself, and pass it with the free option, and then you'll be sure to know what it is.

restart;

randomize():

for iter from 1 to 5 do
  A := LinearAlgebra:-RandomMatrix(10,density=0.15):
  v := Vector(10):

  # This generates an unassigned global name.
  frname := `tools/genglobal`('K');

  sol := LinearAlgebra:-LinearSolve(A,v,free=frname);
  found := indets(sol,'specindex'(frname));
  print(found);
end do:

{K[4], K[5], K[8]}

{K0[4], K0[7]}

{K1[3], K1[5], K1[8], K1[9]}

{K2[3], K2[5]}

{K3[10]}

 

Download specindexLS.mw

You don't need to generate discrete data (and interpolate it), in order to obtain the smooth plot of the lambda operator.

The smooth plot of lambda is assigned to P2, below. You can use the so-called parametric calling sequence of the plot command to get the orientation for lambda(t) versus t.

   plot( [ lambda(t), t, t=0..1 ] );

You don't even need the discrete data here. I included it simply to overlay the point-plot with the smooth curve, to show you that they agree.

restart

lambda := proc (x) options operator, arrow; (-1.565845910+2.393779704*x^2+1.564996800*x^4+1.800900000*x^6)/(x^2+2)^4 end proc

proc (x) options operator, arrow; (-1.565845910+2.393779704*x^2+1.564996800*x^4+1.800900000*x^6)/(x^2+2)^4 end proc

data := [seq([lambda(x), x], x = 0 .. 1, .1)]

[[-0.9786536938e-1, 0], [-0.9445602702e-1, .1], [-0.8473253124e-1, .2], [-0.7004169612e-1, .3], [-0.5215959052e-1, .4], [-0.3283205351e-1, .5], [-0.1345044703e-1, .6], [0.5065808511e-2, .7], [0.2221891338e-1, .8], [0.3780341321e-1, .9], [0.5177568635e-1, 1.0]]

P1 := plot(data, style = point, symbol = solidcircle, symbolsize = 15, color = blue)

P2 := plot([lambda(t), t, t = 0 .. 1])

plots:-display(P1, P2)

 

Download help2_acc.mw

Yes, as you mention you could freeze the function calls. The type system lets you target them quite well, and so subsindets is useful.

I don't like using frontend for this, since the set of things not to freeze varies by example. It's ad hoc. Consider what you'd need to do to handle,
   {x(0) + y(0) = 5, x(0) - sin(y(0)) = 3}
And now realize that it's much more straightforward to freeze specifically the function calls (of type constant) in question than it is to programmatically determine example by example what kinds of symbolic calls and subexpressions (depending on those function calls) which must be excluded from being frozen.

sys:={x(0) + y(0) = 5, x(0) - y(0) = 3}:

thaw(fsolve(subsindets(sys,And(anyfunc(constant),
                               Not(complex(realcons))),
                       freeze)));

              {x(0) = 4., y(0) = 1.}

sys:={x(0) + y(0) = 5, x(0) - sin(y(0)) = 3}:

thaw(fsolve(subsindets(sys,And(anyfunc(constant),
                               Not(complex(realcons))),
                       freeze)));

            {x(0) = 3.893939842, y(0) = 1.106060158}

Im not a big fan of the idea that these function calls (of type constant) should be considered as "variables". It's unclear what problems might arise, were fsolve to support this directly. Next might come unknowns x(a) and y(b), but then fsolve would need be careful to reject examples with unknown x(a) alongside bare unknown a.

Could you draw it twice, with two different padding sizes and the layout the same (so they are overlaid)?

The colors might be chosen with different chroma but similar intensity, so that the text labels are legible.

restart;

with(GraphTheory):

G := CycleGraph([v__1,v__2,v__3,v__4]):
H:=Graph({{u__1,u__2}}):

DrawGraph(G,size=[250,250],stylesheet=[vertexborder=false,vertexpadding=10,
vertexcolor="Orange",edgethickness=3]);

DrawGraph(H,size=[250,250],
          stylesheet=[vertexborder=false,vertexpadding=10,
                      vertexcolor="Cyan",edgethickness=3]);

GH:=CartesianProduct(G,H):

 

plots:-display(
  DrawGraph(GH,style=spring,
            stylesheet=[vertexborder=false,vertexpadding=4,
                        vertexcolor="Orange",edgethickness=3]),
  DrawGraph(GH,style=spring,
            stylesheet=[vertexborder=false,vertexpadding=14,
                        vertexcolor="Cyan",edgethickness=3])
);

 

Download graph_twice.mw

The evalf (kernel builtin) extension mechanism invokes the extension procedure `evalf/XXX` even when the function call to XXX is a call to the local of that name.

restart;

local gamma:

trace(`evalf/gamma`):

evalf(gamma(1));

{--> enter \`evalf/gamma\`, args = 1

<-- exit \`evalf/gamma\` (now at top level) = -0.7281584548e-1}

-0.7281584548e-1

restart;

proc(x)
  local sin;
  evalf(sin(Pi/4));
end proc();

.7071067813

restart;

m := module()
  option package;
  export sin;
end module:
with(m);

[sin]

K:=sin(Pi/4);
op(0,%);
addressof(%), addressof(:-sin);

sin((1/4)*Pi)

sin

18446884039962474558, 18446884039962016894

evalf(K);

.7071067813

 

Download evalf_extension.mw

I have a very vague recollection of reporting/complaining about this -- but that would be almost two decades ago. IIRC I was concerned with some cases of what happened (unexpectly, for me) under evalhf.

You can enter your expression using the inert operator `%*` , as Kitonum shows.

You can also enter your expression as a string, and use InertForm:-Parse (though it's trickier to target just the powers in question, if part of a larger expression).

InertForm:-Parse("x*x*x*x");

`%*`(x, x, x, x)


And you can cobble together mechanisms to decompose. (Adjust or robustify as needed...)

And you can convert back to the active `*` form, in a few ways.

restart;

expr := x^2 + 1/y^5 + (z+1)^2 + Sum(i,i=1..N);

x^2+1/y^5+(z+1)^2+Sum(i, i = 1 .. N)

T := subsindets(expr, anything^integer,
           proc(u) if op(2,u)>0 then
                     `%*`(op(1,u)$op(2,u));
                   else
                     1/`%*`(op(1,u)$(-op(2,u)));
                   end if;
           end proc);

`%*`(x, x)+1/`%*`(y, y, y, y, y)+`%*`(z+1, z+1)+Sum(i, i = 1 .. N)

# To display without the gray operator.

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

`%*`(x, x)+1/`%*`(y, y, y, y, y)+`%*`(z+1, z+1)+Sum(i, i = 1 .. N)

InertForm:-Value(T);

x^2+1/y^5+(z+1)^2+(1/2)*(N+1)^2-(1/2)*N-1/2

# And, without affecting the `Sum`

subsindets(T, specfunc(anything,`%*`), u->`*`(op(u)));

x^2+1/y^5+(z+1)^2+Sum(i, i = 1 .. N)

 

Download inert_pow_prod.mw

First 108 109 110 111 112 113 114 Last Page 110 of 336