acer

32333 Reputation

29 Badges

19 years, 320 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Inside main's catch clause the callstack is no longer what it was when the original exception was thrown.

You can access the procedure name, using lastexception[1]. For example,

restart;

A := module()  
  local B:=module()
     export foo:=proc(n)
         1/n;
     end proc;
  end module;

  export main:=proc()
    try
       B:-foo(0);
    catch:
      print(sprintf("error happened in %a: %s",
            lastexception[1],
            StringTools:-FormatMessage(lastexception[2..-1])));
    end try;
  end proc;

end module:

A:-main()

"error happened in foo: numeric exception: division by zero"

Download try_loc.mw

If your procedures/modules are being accessed from a Library archive into which they were saved (and if the source exists so that line-number information can be used) then you may be able to access that source via lasterrorlocus and showsource. (See the Description of the Help page for topic error.)

For example, I have a file "foo.mpl" with the following plaintext Maple code in it:

A := module()  
  local B:=module()
     export foo:=proc(n)
         1/n;
     end proc;
  end module;
  export main:=proc()
    try
       B:-foo(0);
    catch:
      if lasterrorlocus<>'lasterrorlocus' then
         print(lasterrorlocus);
         showsource();
      end if; 
      print(sprintf("error happened in %a: %s",
            lastexception[1],
            StringTools:-FormatMessage(lastexception[2..-1])));
    end try;
  end proc;
end module:

I can then load that into a .mla Library archive file, and access the relevent souce upon catching an error. For example:

restart;

LB:=cat(currentdir(),"/nm.mla"):
read "foo.mpl";
if FileTools:-Exists(LB) then
  FileTools:-Remove(LB);
end if;
LibraryTools:-Create(LB);
LibraryTools:-Save(A,LB);

restart;

libname := libname,cat(currentdir(),"/nm.mla"):

A:-main()

["foo.mpl", 4, false]

"error happened in A:-main: numeric exception: division by zero"


--- foo.mpl -------------------------------------------------------------------
1  A := module()  
2    local B:=module()
3       export foo:=proc(n)
4!          1/n;
5       end proc;
6    end module;
7  
8    export main:=proc()
9      try

 

Download try_locb.mw

I'll note that the lastexception[1] shows as A:-main in that last example, while more usefully it shows foo in the first example. But the point of the second example is mostly about printing the source snippet with error location indicated.

If you want the result from calling komb to include indexed names then you could form then as such.

Your approach of concatenating produces names that pretty-print with subscripts but are not indexed names.

There are two kinds of name that pretty-print with subscripts:
1) indexed names like a[2,2] for which the subscripts pretty-print in upright Roman characters.
2) special names with double-underscore like say `a__2,2` for which the subscript pretty-prints in italic characters.
These two kinds of name might appear similar, except for the italic aspect and that is easy to overlook. But they are quite distinct from each other (ie. will not cancel upon subtraction).

That approach of concatenation will produce global names. If you want the result from komb to not have escaped local (indexed) names then you can pass in the base-symbol as an additional argument (instead of declaring that base-symbol as local). Below, I do this by having var be an additional parameter of the procedure komb, so that the base-symbol (eg. a) can be passed in. The effect is that the results from komb have the same kind of global name a as are in the Matrix A.

NULL

restart

komb := proc (n::integer, m::integer, var::symbol) local L, j; L := [seq(var[j, j], j = 1 .. n)]; add(mul(k), `in`(k, combinat:-choose(L, m))); return % end proc

komb(5, 5, a)

a[1, 1]*a[2, 2]*a[3, 3]*a[4, 4]*a[5, 5]

lprint(komb(5, 5, a))

a[1,1]*a[2,2]*a[3,3]*a[4,4]*a[5,5]

A := Matrix(5, 5, shape = diagonal, symbol = a)

Matrix(%id = 36893628247904552828)

abs(A)

a[1, 1]*a[2, 2]*a[3, 3]*a[4, 4]*a[5, 5]

lprint(abs(A))

a[1,1]*a[2,2]*a[3,3]*a[4,4]*a[5,5]

simplify(komb(5, 5, a)-abs(A))

0

NULL

Download indexed_ac.mw

The colored portions shown by sparsematrixplot are filled polygons, not symbols, ie. it's not a point-plot, which is why passing the symbol option doesn't work for you.

Here's another approach:

restart;

plots:-setoptions(size=[500,500]):

kernelopts(version);

`Maple 2021.2, X86 64 LINUX, Nov 23 2021, Build ID 1576349`

H := proc(M::Matrix) local i,m,n,L,V;
  uses LA=LinearAlgebra, IT=ImageTools, P=plots;
  (m,n):=LA:-Dimensions(M):
  (L,V):=((u->[lhs(u)[2],m-lhs(u)[1]+1])~,rhs~)([rtable_elems(M)[]]):
  P:-pointplot(L, style=point, symbol=cross, symbolsize=50,
               color=COLOR(HUE,0.75*IT:-FitIntensity(Array(1..nops(V),1..1,V,
                                                           datatype=float[8]))),
               axis[1]=[tickmarks=[seq(floor(i)=floor(i), i=1..n, (n-1)/4)]],
               axis[2]=[tickmarks=[seq(floor(i)=floor(m-i+1), i=1..m, (m-1)/4)]],
               axes=box, view=[0.5..n+0.5,0.5..m+0.5],
               labels=["column","row"], labeldirections=[horizontal,vertical],
               scaling=constrained, _rest);
end proc:

 

A := Matrix([[2, 1, 0, 0, 3],[0, 2, 1, 0, 0],
             [0, 0, 2, 1, 0],[0, 0, 0, 2, 1],[0, 0, 0, 0, 2]]);

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

plots:-sparsematrixplot(A, matrixview, color = "Red", scaling=constrained);

H(A);

B := LinearAlgebra:-RandomMatrix(40,80,generator=0..10,density=0.2):

plots:-sparsematrixplot(B, matrixview, color = "Red", scaling=constrained);

H(B, symbolsize=8);

H(B, symbolsize=8, symbol=solidbox);

 

Download sparse_mat_plot3.mw

Another example:

A := Matrix(100,(i,j)->sin(Pi*(i+j)/20)) *~
     LinearAlgebra:-RandomMatrix(100, generator=1..1,
                                 density=0.1):

H(A, symbolsize=6, symbol=soliddiamond,
  background="#777777", size=[500,500]);

If you do decide that you want just filled polygons rather than symbols -- but colored according to value -- then another approach would be to use the Statistics:-HeatMap command, optionally with a colorscheme option.

It wasn't clear whether you were going to deal with other, larger, Matrices. Once the size gets moderately large (more than 200, say) then switching to a density-style plot using surfdata might be better in terms of GUI performance. And if the size got large (more than 400, say) then switching from a plot to an image might be better.

Here's a way, that has reasonably decent speed to get 500x500 resolution.

It also uses the iteration count to adjust the intensity of the color, retaining a little more information. I like the extra pop it provides. Some alternates approaches are commented out.

The timing is not bad: in total about 4.2 sec for the whole computation on my machine.

But the classifier that matches points to closest root of the polynomial is not optimized, and that step alone takes a 2.2 sec portion of that total! (I couldn't find my old code that did that step better, sorry...)

restart:

kernelopts(version);

`Maple 2021.2, X86 64 LINUX, Nov 23 2021, Build ID 1576349`

with(IterativeMaps): with(ImageTools):

#ee := (z^6-1)/(5*z^2+1)^(1/5);
ee := z^3-1:

st := time[real]():

rts:=Vector[row]( sort( [ fsolve(numer(ee),z,complex) ] ),
                  datatype=complex[8]):
evalf[8](%);

Vector[row]([1.+0.*I, -.500000000000000-.866025400000000*I, -.500000000000000+.866025400000000*I])

(minr,maxr) := -2,2:
(mini,maxi) := -2,2:

dee := diff(ee, z):
Z := eval( z - ee/dee, z=zr+zi*I):

fzi := evalc( Im( Z ) ):
fzi := simplify(combine(fzi)) assuming real:

fzr := evalc( Re( subs(zi=zit, Z) ) ):
fzr := simplify(combine(fzr)) assuming real:

N := 500;
maxiter := 20; # increasing this will reduce nonconvergence
h,w := N, floor(N*(maxr-minr)/(maxi-mini));

500

20

500, 500

newt := Escape( height=h, width=w,
                       [zi, zr, zit],
                       [fzi, fzr, zi],
                       [y, x, y], abs(eval(ee,[z=zr+zit*I]))^2<1e-10,
                       minr, maxr, mini, maxi,
                       iterations = maxiter,
                       redexpression = zr, greenexpression = zit, blueexpression=i ):

 

## Use full saturation and full intensity
img:=Array(1..h,1..w,1..3,(i,j,k)->1.0,datatype=float[8]):

temp:=copy(newt[..,..,3]):

 

## Uncomment either scheme 1 or scheme 2.
## Optional, set intensity to zero for non-convergence.
## scheme 1
#temp:=Threshold(temp,1,low=maxiter+1):
#img[..,..,3]:=Threshold(temp,maxiter-1,high=0,low=1):

 

## Optional, set intensity by count (to converge).
## scheme 2
temp:=ln~(map[evalhf](max,temp,1.0)):
temp2:=Threshold(temp,0.01,low=1):
temp:=zip(`*`,temp,temp2):
temp:=FitIntensity(temp):
img[..,..,3]:=temp:

 

## Optional, scale the saturation by count (to converge).
## This works as a possible addition to scheme 2.
#img[..,..,2]:=1-~img[..,..,3]:

 

## Done inefficiently, set the hue by position in the list of
## all roots.
G:=Array(zip((a,b)->min[index](map[evalhf](abs,rts-~(a+b*I))),
             newt[..,..,1],newt[..,..,2]),
         datatype=float[8],order=C_order):
img[..,..,1]:=240*FitIntensity(G):

 

final := HSVtoRGB(Flip(img,'vertical')):

(time[real]() - st) * 'seconds';

4.208*seconds

Embed(final);

NULL

newtbasin2.mw

The above shows an image, in the ImageTools sense. You could turn that into a kind of density-style plot with axes as follows. And you could optionally display this along with a pointplot of the roots, etc. (But note that rendering such high-resolution density-style plots uses a great deal of GUI resources. Embedding images is much gentler on the GUI.)

PLOT(GRID(-2..2,-2..2,Matrix(N,N,datatype=float[8]),
     COLOR(RGB,Rotate(final,-90,degrees))),
     STYLE(PATCHNOGRID),AXESSTYLE(BOX));

Of course one could optionally scale the intensity (brightness) in reverse, so that points which took longer to converge were darker. (That might jibe better with the non-converged points showing as black.)

And here is the same thing, with the option of using the escape-count to adjust the color saturation as well as its intensity.

And here is the other option, using a simpler "scheme 1" for the coloring.

It looks as if you might be trying to compute the following kind of thing.

(But you might be unclear on the differences between definite and indefinite integration).

restart;

kernelopts(version);

`Maple 2015.2, X86 64 LINUX, Dec 20 2015, Build ID 1097895`

Digits := 15;

15

yyExact := 1/(m*c)*sqrt((m*c)^2+p^2+q^2);

(c^2*m^2+p^2+q^2)^(1/2)/(m*c)

ff11 := cos(k*(p+q+x+t));

cos(k*(p+q+x+t))

JJx11 := simplify(combine( QQ/m * Int(p/yyExact * ff11, [p=0..1, q=0..1]) ));

Int(cos(k*p+k*q+k*t+k*x)*QQ*p*c/(c^2*m^2+p^2+q^2)^(1/2), [p = 0 .. 1, q = 0 .. 1])

eval(JJx11, [t=0, m=9.11e-31, c=3e8, QQ=1.6e-19, k=0.5, x=0]);

Int(0.48e-10*cos(.5*p+.5*q)*p/(p^2+q^2+0.7469289e-43)^(1/2), [p = 0 .. 1, q = 0 .. 1])

evalf(%);

0.265894013394288e-10

# I'm not sure why you were raising Digits as high as 30.
Digits := 30;
eval(Int(op(JJx11), epsilon=1e-7, method=_CubaCuhre),
     [t=0, m=9.11e-31, c=3e8, QQ=1.6e-19, k=0.5, x=0]);
evalf[7](evalf[30](%));
Digits := 15;

30

Int(0.48e-10*cos(.5*p+.5*q)*p/(p^2+q^2+0.7469289e-43)^(1/2), [p = 0 .. 1, q = 0 .. 1], epsilon = 0.1e-6, method = _CubaCuhre)

0.2658940e-10

15

simplify(combine( Int(p/yyExact * ff11, [p=0..1, q=0..1]) ));

Int(cos(k*p+k*q+k*t+k*x)*p*m*c/(c^2*m^2+p^2+q^2)^(1/2), [p = 0 .. 1, q = 0 .. 1])

eval(%, [t=0, m=9.11e-31, c=3e8, QQ=1.6e-19, k=0.5, x=0]);

Int(0.2733e-21*cos(.5*p+.5*q)*p/(p^2+q^2+0.7469289e-43)^(1/2), [p = 0 .. 1, q = 0 .. 1])

evalf(%);

0.151473858761892e-21

 

Download question11_ac.mw

If you can figure out what kind of integration you really intend then perhaps it's worth investigating why you were raising Digits so high.

Perhaps you're after this kind of functionality,

solve( Or(And(x + y = 2,y = 1,x > 0),And(x^2 = 4, x < 0, y = 0)) );
 
         {x = 1, y = 1}, {x = -2, y = 0}

restart;

eq := A=P*(1+r/n)^(n*t);

A = P*(1+r/n)^(n*t)

ans:=combine(simplify(evalindets(solve(eq,n,allsolutions),suffixed(_Z),u->0)));

-r*ln(A/P)/(r*t*LambertW(-ln(A/P)*(A/P)^(-1/(r*t))/(r*t))+ln(A/P))

evalf(eval(eval(eq, n=ans), [P=5,A=10,r=1/10,t=2]));

10. = 10.00000001

Download solve_ex2.mw

This is the same kind of mistake as your other example, and is not a bug.

A successful way is to use map[2], to map applyrule over B using the entire list half_angle_rule as the first argument (ie. passed to applyrule, for each entry in B).

What your orignal code instructs the elementwise operator applyrule~ to do is more like an interleaved zip action.

restart;

half_angle_rule := [
        sin(x::name) = 2*sin(x/2)*cos(x/2),
        cos(x::name) = 1 - 2*sin(x/2)^2
]:

A := < sin(u), sin(u) >;

A := Vector(2, {(1) = sin(u), (2) = sin(u)})

map[2](applyrule, half_angle_rule, A);

Vector[column]([[2*sin((1/2)*u)*cos((1/2)*u)], [2*sin((1/2)*u)*cos((1/2)*u)]])

applyrule~(half_angle_rule, A);

Vector[column]([[2*sin((1/2)*u)*cos((1/2)*u)], [sin(u)]])

zip(applyrule, <half_angle_rule>, A);

Vector[column]([[2*sin((1/2)*u)*cos((1/2)*u)], [sin(u)]])

Download map2_ex2.mw

The error message states the problem. It is not a bug.

You goal appears to be to map applyrule over each of the entries of C, with all of double_angle_rule as the first argument.

The elementwise operator applyrule~ is not correct for doing that. As you wrote it Maple is trying something more like an interleaved zip operation, but the dimensions of C and double_angle_rule do not match.

One way to accomplish the presumed goal is to utilize map[2].

restart;

double_angle_rule := [
        sin(x::name/2)*cos(x::name/2) = 1/2*sin(x),
        sin(x::name/2)^2 = 1/2*(1-cos(x)),
        cos(x::name/2)^2 = 1/2*(1+cos(x))
]:

C := < cos(1/2*u)*sin(1/2*u), cos(1/2*u)^2 >;

C := Vector(2, {(1) = cos((1/2)*u)*sin((1/2)*u), (2) = cos((1/2)*u)^2})

map[2](applyrule, double_angle_rule, C);

Vector[column]([[(1/2)*sin(u)], [1/2+(1/2)*cos(u)]])

applyrule(double_angle_rule, C[1]);

(1/2)*sin(u)

applyrule(double_angle_rule, C[2]);

1/2+(1/2)*cos(u)

Download map2_ex.mw

That message (about int being protected) can occur with Grid:-Map in at least in some old versions, eg. Maple 16.

You might pass a user-defined procedure instead, eg.

restart;

kernelopts(version);

`Maple 16.01, X86 64 LINUX, May 6 2012, Build ID 744592`

A := Matrix([[x,x^2],[2*x,sin(x)]]);

A := Matrix(2, 2, {(1, 1) = x, (1, 2) = x^2, (2, 1) = 2*x, (2, 2) = sin(x)})

Grid:-Map(e->int(e,x=0..2,numeric=true),A);

Matrix(2, 2, {(1, 1) = 2., (1, 2) = 2.666666667, (2, 1) = 4., (2, 2) = 1.416146837})

Download Gr_M16.mw

You could find both tangents to the hyperbola at x=0.504244923, and then select which touches close to y=0.3781836925.

restart;

(x0,y0) := 0.504244923, 0.3781836925:

expr := 7*x^2-7*y^2-12.0*x+9.0*y+2.25-2*x*y:

both := Student:-Calculus1:-Tangent~([solve(expr,y)],x=x0):

ans := y=~select(t->abs(eval(t,[x=x0])-y0)<1e-5,both)[];

y = 2.112372426*x-.6869693794

plots:-display(
  plots:-pointplot([[x0,y0]],symbol=solidcircle,symbolsize=15),
  plot(eval~(y,[ans]),x=-1..2,color=blue),
  plots:-implicitplot(expr,x=0..2,y=-0.5..1.5),
  view=[0..2,-0.5..1.5],scaling=constrained,size=[300,300]);

Download hyp_tang.mw

The Compiler does not know how to handle the syntax ++n and val += FN, and that's where that error message arises.

But it hardly matters, since the Compiler is not aware of combinat:-fibonnacci. You could have that problematic statement "escape" temporarily back to Maple proper, by wrapping in an eval call.

But adjusting the syntax and escaping the combinat:-fibonacci call will not -- in itself -- result in faster performance.

This looks like a contrived example, perhaps constructed to investigate compilation, etc. If so, then it might be better to tell us any other goals in detail.

It can be shown to be true (via is and combine) in a 1-line statement for your conditions r>0, r<1. (Also for the wider r>-1, r<1.)

C1 := sqrt(r+1)*sqrt(1/(r-1)):
C2 := -sqrt(r^2-1)/(r-1):

is(combine(C1-C2)=0) assuming r>0, r<1;

             true

is(combine(C1-C2)=0) assuming r>-1, r<1;

             true

note: Such situations in which an extra call to combine (or simplify, say) is needed (in order for is to get the result) are quite rare in my experience. Like obscure corners of lepidopterology. I have reported only a few such, over the years. I'll submit a report for this one too.

I prefer solutions in which the "4" does not get rendered in italics (as if it were a name...).

plots:-loglogplot(x, x = 0.1*10^(-5) .. 10,
                  title = 10^Typesetting:-Typeset(-4));

plots:-loglogplot(x, x = 0.1*10^(-5) .. 10,
                  title = 10^(-`#mn(4)`));

plots:-loglogplot(x, x = 0.1*10^(-5) .. 10,
                  title = `#msup(mn("10"),mrow(mo("&uminus0;"),mn("4")));`);

plots:-loglogplot(x, x = 0.1*10^(-5) .. 10,
                  title = InertForm:-Display(`%^`(10,-4),inert=false));

You won't be able to get 23333.33 after rounding to 6 digits, since that requires 7 digits.

Here are a few ways:

restart;

K := evalf[6](Matrix(3,[1/3,-20/3,-20/3,200/3,70000/3,1.44]));

K := Matrix(3, 3, {(1, 1) = .333333, (1, 2) = -6.66667, (1, 3) = -6.66667, (2, 1) = 66.6667, (2, 2) = 23333.3, (2, 3) = 1.44, (3, 1) = 0., (3, 2) = 0., (3, 3) = 0.})

map(x->parse(sprintf("%.2f",x)),K);

Matrix([[.33, -6.67, -6.67], [66.67, 23333.30, 1.44], [0., 0., 0.]])

K := evalf[10](Matrix(3,[1/3,-20/3,-20/3,200/3,70000/3,1.44]));

K := Matrix(3, 3, {(1, 1) = .3333333333, (1, 2) = -6.666666667, (1, 3) = -6.666666667, (2, 1) = 66.66666667, (2, 2) = 23333.33333, (2, 3) = 1.44, (3, 1) = 0., (3, 2) = 0., (3, 3) = 0.})

new := map(x->parse(sprintf("%.2f",x)),K):
new;

Matrix([[.33, -6.67, -6.67], [66.67, 23333.33, 1.44], [0., 0., 0.]])

The result here is actually as its displayed.

new[1,2]

-6.67

You could also get K printed with the effect. This does
not affect the actual values.

interface(displayprecision=2):
K;
interface(displayprecision=-1):

Matrix([[.3333333333, -6.666666667, -6.666666667], [66.66666667, 23333.33333, 1.44], [0., 0., 0.]])

A similar effect can be obtained by right-clicking on the
following outout and changing its numeric-formatting
(in-place) to "Fixed" with 2 decimal places.

K;

Matrix(3, 3, {(1, 1) = .33, (1, 2) = -6.67, (1, 3) = -6.67, (2, 1) = 66.67, (2, 2) = 23333.33, (2, 3) = 1.44, (3, 1) = 0., (3, 2) = 0., (3, 3) = 0.})

K[1,2];

-6.666666667

Download Matrix_dec.mw

First 73 74 75 76 77 78 79 Last Page 75 of 336