acer

32405 Reputation

29 Badges

19 years, 351 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

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

You should be able to use the results from the ifactors command, programmatically. (You can trivially massage its output into similar formats of your own invention.)

ifactors(725);

        [1, [[5, 2], [29, 1]]]

ifactors(1125);

         [1, [[3, 2], [5, 3]]]

ifactors(2048);

             [1, [[2, 11]]]

If I understand your sentences rightly then yes, the inner procedure that is defined in the body of the outer procedure will be recreated each time the outer is invoked/called, and this is unnecessary overhead (computational effort and memory management) that could be avoided.

If instead the procedure A were simply defined at a higher level than procedure B, then you should be able to call A from within B, under the lexical scoping rules.

You could also implement both A and B within the same module, if you would prefer some compartmentalization.

You could also make the outer procedure be the local ModuleApply of a module B, where A is another local procedure within that module. That would make B into a so-called appliable module, where the name B could be called with arguments.

Would you like examples of all this? Have you read the Programming Guide (esp. the section on modules?)

If you call the ZeroMatrix command without specifying its option

    compact=false

then (by default) the Matrix will be constructed with a special indexing function.

The effect of that is that the Matrix can only contain zeroes. (This is great for efficiency in some programming situations, since very large zero-matrices can take up little/constant memory.)

Much simpler is just to use the Matrix command instead of the ZeroMatrix command. The default fill-values of a call to Matrix are zero, and those don't have to be specified explicitly.

This is a guess, since you did not provide you code. It's almost always better to provide full code to reproduce. (The green up arrow in the Mapleprimes editor allows you to upload and attach documents. )

restart;

MZ1 := LinearAlgebra:-ZeroMatrix(3, 4);

Matrix(3, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0})

MZ1[2,3] := 17;

Error, attempt to assign non-zero to zero Matrix entry

MZ2 := LinearAlgebra:-ZeroMatrix(3, 4, compact=false);

Matrix(3, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0})

MZ2[2,3] := 17;

17

MZ2;

Matrix(3, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 17, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0})

MZ3 := Matrix(3, 4);

Matrix(3, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0})

MZ3[2,3] := 17;

17

MZ3;

Matrix(3, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 17, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0})

 

If I set the "rtable" (ie. Matrix) associated with the above DataTable to be MZ1 then I'll get a popup error "Unable to change value" whenever I attempt to edit an entry and hit the Return/Enter key,

 

But if I associate MZ2 or MZ3 with the DataTable instead then I can edit the entries.

 

Creating MZ3 as above is much simpler than creating MZ2 as above.

 

Download zeromatrix_edit.mw

For the 2D case here are two short versions. (They could be shorter still, but I wanted to show a few of the intermediate assignments.)

First, one using operators.

with(plots):
F:=(x,y)->4*x^2+9*y^2;
G:=[D[1],D[2]](F)(2,1);
display(arrow([2,1],G),
        implicitplot(F(x,y)=F(2,1),x=-4..4,y=-2..2),
        scaling=constrained);

And one using expression form.

with(plots):
f:=4*x^2+9*y^2;
fval:=eval(f,[x=2,y=1]);
gval:=eval(zip(diff,f,[x,y]),[x=2,y=1]);
display(arrow([2,1],gval),
        implicitplot(f=fval,x=-4..4,y=-2..2),
        scaling=constrained);

Both of those can be adjusted using various options for the arrow dimensions, etc.

And 3D versions can also be made. The 2D implicit plot can be raised to a level height on the 3D surface, or a single contour line can be specified at the specific height. See attached.

gradient_arrow_operator.mw

gradient_arrow.mw

I think that your example is interesting, since the "size" of the desired target expression is smaller than that of the original expression, using any of these as norm:
   length`simplify/size/size`MmaTranslator:-Mma:-LeafCount
So it seems reasonable to ask why the "smaller" form is not found directly by some flavour of call to simplify.

My compressor oracle found these two ways:

ex1:=(2*sqrt(N__s) + N__s + 1)/(2*N__s + 2);

(2*N__s^(1/2)+N__s+1)/(2*N__s+2)

eval(collect(collect(ex1,[N__s],__P),[__P]),
     __P=(u->u));

1/2+N__s^(1/2)/(N__s+1)

thaw(collect(subsindets(ex1,radical,freeze),
     [freeze~(indets(ex1,radical))[]],simplify));

1/2+N__s^(1/2)/(N__s+1)

Download simp_radical_ex.mw

You probably already realize this, but you should not use the symbolic option of the simplify command if you care about results branch cuts and always being correct. However, your current example and desired target are equivalent. I just wanted to mention it, since you made it sound as if simplify(...,symbolic) might be generally applicable/valid.

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