acer

32363 Reputation

29 Badges

19 years, 332 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Are these the kinds of things you're trying to accomplish?

restart;

f := -4*sin(x) + 2*exp(y^2) + 5 - 5*cos(x^3)*sin(y^2) + 5*sinh(x^2);

-4*sin(x)+2*exp(y^2)+5-5*cos(x^3)*sin(y^2)+5*sinh(x^2)

p := proc(e)
  indets(f,And(`*`,satisfies(u->member(e,[op(u)]))))[];
end proc:

 

p(sin(x));

-4*sin(x)

p(sin(y^2));

-5*cos(x^3)*sin(y^2)

p(exp(y^2));

2*exp(y^2)

p(sinh(x^3)); # NULL

 

q := proc(e)
  indets(f,And(`*`,satisfies(u->membertype(specfunc(e),[op(u)]))));
end proc:

q(sin);

{-5*cos(x^3)*sin(y^2), -4*sin(x)}

q({sin,cos,sinh,exp});

{-5*cos(x^3)*sin(y^2), 2*exp(y^2), -4*sin(x), 5*sinh(x^2)}

map(q,[sin,cos,sinh,exp]);

[{-5*cos(x^3)*sin(y^2), -4*sin(x)}, {-5*cos(x^3)*sin(y^2)}, {5*sinh(x^2)}, {2*exp(y^2)}]


Download indets_ex.mw

If your Matrix M is datatype=float[8] then you could accelerate using either evalhf or Compile.

I use one of dharr's quicker procedures for the evalhf variant, but now faster.

The Compiled version is about 30 times faster than the "original"; it varies by run, and version.

I keep the original non-hardware-datatype constructions only for comparison.

restart;

kernelopts(version);

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

N := 10^5:
L := 5:
M := LinearAlgebra:-RandomMatrix(N, L, generator=0..1):

P := [1,0,1]:
nC := numelems(P):

C:=[2,4,5];

[2, 4, 5]

Original

pattern := proc(n) local k;
  `if`(evalb(`and`(seq(M[n, C[k]]=P[k], k=1..nC))), 1, 0);
end proc:
t0 := time[real]():
add(pattern(n), n=1..N);
time[real]()-t0

12585

.256

Suppose we had instead constructed the data this way.

M8:=Matrix(M,datatype=float[8]):
C8:=Vector(C,datatype=float[8]):
P8:=Vector(P,datatype=float[8]):

pattern:=proc(n,nC,M,C,P)
  local k;
  for k to nC do
    if M[n,C[k]]<>P[k] then return 0 end if;
  end do;
  1
end proc:
t0 := time[real]():

trunc(evalhf(add(pattern(n,nC,M8,C8,P8), n=1..N)));
time[real]()-t0

12585

0.49e-1

Suppose we had instead construced the data in this way

M8:=Matrix(M,datatype=float[8]):
C4:=Vector(C,datatype=integer[4]):
P4:=Vector(P,datatype=integer[4]):

pattern:=proc(N::integer,nC::integer,M::Matrix(datatype=float[8]),
                                     C::Vector(datatype=integer[4]),
                                     P::Vector(datatype=integer[4]))
  local n::integer, k::integer, S::integer, ok::integer;
  S := 0;
  for n from 1 to N do
    ok := 1;
    for k to nC do
      if M[n,C[k]]<>P[k] then k := nC; ok := 0; end if;
    end do;
    if ok=1 then S := S + 1; end if;
  end do;
  S;
end proc:

cpattern:=Compiler:-Compile(pattern, 'inmem'=false):
t0 := time[real]():
cpattern(N,nC,M8,C4,P4);
time[real]()-t0

12585

0.9e-2


Download tests_ac1.mw

These are just toy examples, deliberately intended for illustration and comparison rather than some "best" way to do the computation. (Please don't respond with a better way to do it. That's not my point.)

One possibility is to split out your purely floating-point and evalhf'able computations into their own dedicated procedure(s), and then do any non-evalhf'able or symbolic operations separately. I don't bother showing that for the following example, because it's trivial.

But sometimes such splitting involves complicated rewriting of one's existing procedures (which already mix those two aspects in more involved ways). Sometimes one of the following simpler code revisions is possible instead of such inconvenience.


No special acceleration.

restart;
f0 := proc()
  local x,res,y;
  res := add( sin(1.0*x), x=1..10^5);
  res + fsolve(y^3-1);
end proc:
CodeTools:-Usage( f0() );

memory used=0.62GiB, alloc change=33.00MiB, cpu time=4.38s, real time=4.39s, gc time=337.42ms

2.847777198


A procedure can "escape" evalhf mode temporarily, even for
non-evalhf'able computation (which returns a float result or NULL,
yet which might internally use some symbolic names), by using eval.

restart;
f1 := proc()
  local x,res,y;
  res := add( sin(1.0*x), x=1..10^5);
  res + eval(fsolve(y^3-1));
end proc:
CodeTools:-Usage( evalhf(f1()) );

memory used=2.86MiB, alloc change=0 bytes, cpu time=35.00ms, real time=36.00ms, gc time=0ns

2.84777710363033965


By having option hfloat on the procedure, its float computations
get done using HFloat's, which can be fast because quite a lot of
arithmetic, and elementary functions can deal with those quickly.

It's not as fast as operating in evalhf mode, but often can be faster
for some basic floating-point computations than when done with no
special acceleration.

restart;
f2 := proc()
  option hfloat;
  local x,res,y;
  res := add( sin(1.0*x), x=1..10^5);
  res + fsolve(y^3-1);
end proc:
CodeTools:-Usage( f2() );

memory used=12.28MiB, alloc change=2.00MiB, cpu time=114.00ms, real time=114.00ms, gc time=8.30ms

HFloat(2.847777103630354)

 

Download hfloat_evalhf_none.mw

restart;

TSelided := proc(ee, L:=infinity) local foo, oldom, u;
   oldom := kernelopts(':-opaquemodules'=false);
   try foo := SubexpressionMenu:-truncateTypeMK(ee, L);
   catch: finally kernelopts(':-opaquemodules'=oldom); end try;
   eval(parse(nprintf("%s",(String(foo)[3..-2]))),
     [seq(u=Typesetting[u],
          u=[':-mfenced',':-mfrac',':-mi',':-mio',':-mn',':-mo',':-mover',':-mroot',':-mrow',
             ':-ms',':-mscripts',':-mspace',':-msqrt',':-mstyle',':-msub',':-msubsup',
             ':-msup',':-mtable',':-mtd',':-mtext',':-mtr',':-munder',':-munderover'])]);
end proc:

QQFProj := proc(q12::algebraic, q23::algebraic,
                q34::algebraic, q14::algebraic,
                prnt::boolean:=true, #{columns:=[QQFproj,Q13proj,Q24proj]}
                L::posint:=32)
  description "Projective quadruple quad formula and intermediate 13 and 24 quads.
               Useful for cyclic quadrilaterals";
  local qqf,q13,q24, sub1,sub2,sub3, R,values,DF,lens;
  uses DocumentTools;
  sub1:= (q12 + q23 + q34 + q14)^2 - 2*(q12^2 + q23^2 + q34^2 + q14^2) ;
  sub2:=-4*(q12*q23*q34+q12*q23*q14+q12*q34*q14+q23*q34*q14)+8*q12*q23*q34*q14;
  sub3:=64*q12*q23*q34*q14*(1-q12)*(1-q23)*(1-q34)*(1-q14);
  qqf:=((sub1+sub2)^2=sub3);
  q13:=((q12-q23)^2-(q34-q14)^2)/(2*(q12+q23-q34-q14-2*q12*q23+2*q34*q14));#check this
  q24:=((q23-q34)^2-(q12-q14)^2)/(2*(q23+q34-q12-q14-2*q23*q34+2*q12*q14));#check this
  if prnt then
   values:=<qqf,q13,q24>;
   DF:=DataFrame(TSelided~(<values>,L), columns=[`"Values Equations"`],
                 rows=[`#1  QQF`,`#2  Q13`,`#3  Q24`]);
   lens := [4 +8* max(op(length~(RowLabels(DF)))),
            min(ceil(500*sqrt(L)/max(op(length~(RowLabels(DF))))),1000)];
   Tabulate(DF,width=add(lens),widthmode = pixels,weights = lens);
  end if;
  return qqf,q13,q24
end proc:

 

q12:=1/2:q23:=9/10:q34:=25/26:q41:=9/130: #Cyclic quadrilateral
q12:=sqrt(17+a)/2:q23:=expand(r^2+t^2)^2/10:q34:=expand((a+b+c)^4/26):q41:=sqrt(17+b)/130:

 

Q := [QQFProj(q12,q23,q34,q41,true)];

 

QQFProj(q12,q23,q34,q41,true,80);

QQFProj(q12,q23,q34,q41,true,12);

Q[1];

(((1/2)*(17+a)^(1/2)+(1/10)*(r^2+t^2)^2+(1/26)*a^4+(2/13)*a^3*b+(2/13)*a^3*c+(3/13)*a^2*b^2+(6/13)*a^2*b*c+(3/13)*a^2*c^2+(2/13)*a*b^3+(6/13)*a*b^2*c+(6/13)*a*b*c^2+(2/13)*a*c^3+(1/26)*b^4+(2/13)*b^3*c+(3/13)*b^2*c^2+(2/13)*b*c^3+(1/26)*c^4+(1/130)*(17+b)^(1/2))^2-35921/4225-(1/2)*a-(1/50)*(r^2+t^2)^4-2*((1/26)*a^4+(2/13)*a^3*b+(2/13)*a^3*c+(3/13)*a^2*b^2+(6/13)*a^2*b*c+(3/13)*a^2*c^2+(2/13)*a*b^3+(6/13)*a*b^2*c+(6/13)*a*b*c^2+(2/13)*a*c^3+(1/26)*b^4+(2/13)*b^3*c+(3/13)*b^2*c^2+(2/13)*b*c^3+(1/26)*c^4)^2-(1/8450)*b-(1/5)*(17+a)^(1/2)*(r^2+t^2)^2*((1/26)*a^4+(2/13)*a^3*b+(2/13)*a^3*c+(3/13)*a^2*b^2+(6/13)*a^2*b*c+(3/13)*a^2*c^2+(2/13)*a*b^3+(6/13)*a*b^2*c+(6/13)*a*b*c^2+(2/13)*a*c^3+(1/26)*b^4+(2/13)*b^3*c+(3/13)*b^2*c^2+(2/13)*b*c^3+(1/26)*c^4)-(1/650)*(17+a)^(1/2)*(r^2+t^2)^2*(17+b)^(1/2)-(1/65)*(17+a)^(1/2)*((1/26)*a^4+(2/13)*a^3*b+(2/13)*a^3*c+(3/13)*a^2*b^2+(6/13)*a^2*b*c+(3/13)*a^2*c^2+(2/13)*a*b^3+(6/13)*a*b^2*c+(6/13)*a*b*c^2+(2/13)*a*c^3+(1/26)*b^4+(2/13)*b^3*c+(3/13)*b^2*c^2+(2/13)*b*c^3+(1/26)*c^4)*(17+b)^(1/2)-(1/325)*(r^2+t^2)^2*((1/26)*a^4+(2/13)*a^3*b+(2/13)*a^3*c+(3/13)*a^2*b^2+(6/13)*a^2*b*c+(3/13)*a^2*c^2+(2/13)*a*b^3+(6/13)*a*b^2*c+(6/13)*a*b*c^2+(2/13)*a*c^3+(1/26)*b^4+(2/13)*b^3*c+(3/13)*b^2*c^2+(2/13)*b*c^3+(1/26)*c^4)*(17+b)^(1/2)+(1/325)*(17+a)^(1/2)*(r^2+t^2)^2*((1/26)*a^4+(2/13)*a^3*b+(2/13)*a^3*c+(3/13)*a^2*b^2+(6/13)*a^2*b*c+(3/13)*a^2*c^2+(2/13)*a*b^3+(6/13)*a*b^2*c+(6/13)*a*b*c^2+(2/13)*a*c^3+(1/26)*b^4+(2/13)*b^3*c+(3/13)*b^2*c^2+(2/13)*b*c^3+(1/26)*c^4)*(17+b)^(1/2))^2 = (8/325)*(17+a)^(1/2)*(r^2+t^2)^2*((1/26)*a^4+(2/13)*a^3*b+(2/13)*a^3*c+(3/13)*a^2*b^2+(6/13)*a^2*b*c+(3/13)*a^2*c^2+(2/13)*a*b^3+(6/13)*a*b^2*c+(6/13)*a*b*c^2+(2/13)*a*c^3+(1/26)*b^4+(2/13)*b^3*c+(3/13)*b^2*c^2+(2/13)*b*c^3+(1/26)*c^4)*(17+b)^(1/2)*(1-(1/2)*(17+a)^(1/2))*(1-(1/10)*(r^2+t^2)^2)*(1-(1/26)*a^4-(2/13)*a^3*b-(2/13)*a^3*c-(3/13)*a^2*b^2-(6/13)*a^2*b*c-(3/13)*a^2*c^2-(2/13)*a*b^3-(6/13)*a*b^2*c-(6/13)*a*b*c^2-(2/13)*a*c^3-(1/26)*b^4-(2/13)*b^3*c-(3/13)*b^2*c^2-(2/13)*b*c^3-(1/26)*c^4)*(1-(1/130)*(17+b)^(1/2))

TSelided(Q[1], 32);

(((1/2)*sqrt(17+a)+(1/10)*(r^2+t^2)^2+(1/26)*a^4+2*a^3*`&hellip;`*(1/13)+`&hellip;`)^2-`&hellip;`)^2 = `&hellip;`

 

 

2024-05-14_Q_Data_Table_Limit_Equation_lentgh_ac.mw

 

You asked, "What is the problem?"

The problem is that the answer depends on the variables. The answer changes according to the signs of both a and b, even if all of a,b,h are taken as real.

restart;

r := a + (b - a)*z/h:

x1 := sqrt(r^2 - y^2):


The result depends on whether a and b are positive or negative.

simplify( int(int(int(1, x = -x1 .. x1), y = -r .. r), z = 0 .. h) )
   assuming h::real, a>0, b>0;

(1/3)*(a^2+a*b+b^2)*Pi*h

simplify( int(int(int(1, x = -x1 .. x1), y = -r .. r), z = 0 .. h) )
   assuming h::real, a<0, b>0;

-Pi*(a^3+b^3)*h/(-3*b+3*a)

simplify( int(int(int(1, x = -x1 .. x1), y = -r .. r), z = 0 .. h) )
   assuming h::real, a<0, b<0;

-(1/3)*(a^2+a*b+b^2)*Pi*h


Download trip_int_ex.mw

Your call InitState(5,5) calls Coherent(1,Pi/2,5) and that returns a 1-Vector.

Your call Coherent(5) returns a 5-Vector.

You've set up both InitState and Coherent with keyword parameters. You haven't set it up with optional ordered parameters.

When your InitState calls Coherent(1,Pi/2,5) then that "1" becomes the value of n_max (the first required positional parameter). The other two arguments, Pi/2 and 5, are thrown away. In particular, that's not the way to pass them along for the zeta and phi keyword parameters of your Coherent.

Perhaps you instead intended something like this?

InitState := proc({zeta:=1,phi:=Pi/2,L:=1,sigma:=0.01,beta:=0.02,k0:=-5*Pi},d1,d2) 
  local z1, Z0;
  z1 := Coherent(d1,':-zeta'=zeta,':-phi'=phi);
end proc:

That would utilize the value of d1 for the required parameter of Coherent. It would also pass along the values of zeta and phi to Coherent.

See next attachment for more detail.

test3_2_ac.mw

If you wanted to use optional ordered parameters instead of keyword options then you'd also have to shift the required positional parameters to the front. See next attachment. Notice the absence of {...} braces around the declarations.

test3_2_acc.mw

See also here.

Your previous Question was marked as Maple 2023, so I'm guessing that's the version you're using now.

In Maple 2024.0 the plot3d command supports the colorbar option (introduced for 2D plots in M2023).

But in fact one can attach a COLBAR plotting substructure into a PLOT3D plotting structure, even in M2023.

In the following example the image of the plot is from a single file that I manually exported using right-click. That single image file contains both surface and colorbar.

Here I used the 2D contourplot command as a quick way to construct the colorbar. But it sound as if you've already constructed your own 2D color-bar (for use in your Tabulate). Just wrap that 2D color-bar plot in a call to COLBAR, display that to the desired smaller size, and add that as a last argument to the PLOT3D structure.

restart;

kernelopts(version);

`Maple 2023.2, X86 64 LINUX, Nov 24 2023, Build ID 1762575`

ars := sin(x*y), x=-3..3, y=-3..3, colorscheme=["Red","Green","Blue"]:
PLOT3D(op(plot3d(ars,size=[400,400])),
       subsindets(indets(plots:-contourplot(ars),specfunc(COLBAR))[1],
                  specfunc(PLOT),P->plots:-display(P,size=[200,200])));

 

Download colrbar3d.mw

You cannot do it directly as some option to Tabulate.

But you can programmatically get a zoom by using some of the internal commands which Tabulate uses, and making use of an additional option of the InlinePlot constructor.

restart;
with(plots,display): with(plottools,cylinder):

c1 := display(cylinder([1, 1, 1], 1, 3), orientation = [45, 70],
              scaling = constrained, color = red, size = [300, 300]):
c2 := display(cylinder([1, 1, 1], 1, 3), orientation = [45, 70],
              scaling = constrained, color = blue, size = [300, 300]):

with(DocumentTools): with(DocumentTools:-Layout):

InsertContent(Worksheet(
  Table(exterior = none, interior = none,
        Row(Cell(InlinePlot(c1, scale = 1.7)),
            Cell(InlinePlot(c2, scale = 0.4)))))):

However, I'm not sure that the size=[300,300] will be preserved properly upon worksheet Close/Reopen. It's possible the DocumentTools:-Components could also be used, with the Layout:-InlinePlot call further wrapped in a Components:-Plot call that forced the 300x300 dimensions as width and height. Then suppress the component boundaries, etc. (Sorry, I was once more familiar with these subtleties, but it's been a while...)

I would not be surprised if Edgardo puts a fix into a Physics update.

But, here is a hot patch (for Maple 2024.0) that lets ODESteps handle that example.

restart;

kernelopts(opaquemodules=false):
unprotect(Student:-ODEs:-IVPOrder1):
Student:-ODEs:-IVPOrder1:=FromInert(subs(_Inert_FUNCTION(_Inert_ASSIGNEDNAME("indets", "PROC", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))))), _Inert_EXPSEQ(_Inert_LOCAL(6), _Inert_NAME("name", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected")))))))=_Inert_FUNCTION(_Inert_ASSIGNEDNAME("indets", "PROC", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))))), _Inert_EXPSEQ(_Inert_LOCAL(6), _Inert_FUNCTION(_Inert_NAME("And", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))))), _Inert_EXPSEQ(_Inert_NAME("name", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))))), _Inert_FUNCTION(_Inert_NAME("Not", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))))), _Inert_EXPSEQ(_Inert_NAME("constant", _Inert_ATTRIBUTE(_Inert_NAME("protected", _Inert_ATTRIBUTE(_Inert_NAME("protected"))))))))))),ToInert(eval(Student:-ODEs:-IVPOrder1)))):
protect(Student:-ODEs:-IVPOrder1):
kernelopts(opaquemodules=true):

 

kernelopts(version);

`Maple 2024.0, X86 64 LINUX, Mar 01 2024, Build ID 1794891`

 

ode:=diff(y(x),x)+x*y(x)=1;
ic:=y(0)=0;

diff(y(x), x)+x*y(x) = 1

y(0) = 0

dsolve([ode,ic]);

y(x) = -((1/2)*I)*exp(-(1/2)*x^2)*Pi^(1/2)*2^(1/2)*erf(((1/2)*I)*2^(1/2)*x)

Student:-ODEs:-ODESteps([ode,ic]);

"[[,,"Let's solve"],[,,[(&DifferentialD;)/(&DifferentialD;x) y(x)+x y(x)=1,y(0)=0]],["&bullet;",,"Highest derivative means the order of the ODE is" 1],[,,(&DifferentialD;)/(&DifferentialD;x) y(x)],["&bullet;",,"Isolate the derivative"],[,,(&DifferentialD;)/(&DifferentialD;x) y(x)=-x y(x)+1],["&bullet;",,"Group terms with" y(x) "on the lhs of the ODE and the rest on the rhs of the ODE"],[,,(&DifferentialD;)/(&DifferentialD;x) y(x)+x y(x)=1],["&bullet;",,"The ODE is linear; multiply by an integrating factor" mu(x)],[,,mu(x) ((&DifferentialD;)/(&DifferentialD;x) y(x)+x y(x))=mu(x)],["&bullet;",,"Assume the lhs of the ODE is the total derivative" (&DifferentialD;)/(&DifferentialD;x) (mu(x) y(x))],[,,mu(x) ((&DifferentialD;)/(&DifferentialD;x) y(x)+x y(x))=((&DifferentialD;)/(&DifferentialD;x) mu(x)) y(x)+mu(x) ((&DifferentialD;)/(&DifferentialD;x) y(x))],["&bullet;",,"Isolate" (&DifferentialD;)/(&DifferentialD;x) mu(x)],[,,(&DifferentialD;)/(&DifferentialD;x) mu(x)=mu(x) x],["&bullet;",,"Solve to find the integrating factor"],[,,mu(x)=(e)^((x^2)/2)],["&bullet;",,"Integrate both sides with respect to" x],[,,[]=&int;mu(x) &DifferentialD;x+C1],["&bullet;",,"Evaluate the integral on the lhs"],[,,mu(x) y(x)=&int;mu(x) &DifferentialD;x+C1],["&bullet;",,"Solve for" y(x)],[,,y(x)=(&int;mu(x) &DifferentialD;x+C1)/(mu(x))],["&bullet;",,"Substitute" mu(x)=(e)^((x^2)/2)],[,,y(x)=(&int;(e)^((x^2)/2) &DifferentialD;x+C1)/((e)^((x^2)/2))],["&bullet;",,"Evaluate the integrals on the rhs"],[,,y(x)=(-(&ImaginaryI; sqrt(Pi) sqrt(2) erf(&ImaginaryI;/2 sqrt(2) x))/2+C1)/((e)^((x^2)/2))],["&bullet;",,"Simplify"],[,,y(x)=-((&ImaginaryI; sqrt(Pi) sqrt(2) erf(&ImaginaryI;/2 sqrt(2) x)-2 C1) (e)^(-(x^2)/2))/2],["&bullet;",,"Use initial condition" y(0)=0],[,,0=C1],["&bullet;",,"Solve for" `c__1`],[,,C1=0],["&bullet;",,"Substitute" `c__1`=0 "into general solution and simplify"],[,,y(x)=-&ImaginaryI;/2 (e)^(-(x^2)/2) sqrt(Pi) sqrt(2) erf(&ImaginaryI;/2 sqrt(2) x)],["&bullet;",,"Solution to the IVP"],[,,y(x)=-&ImaginaryI;/2 (e)^(-(x^2)/2) sqrt(Pi) sqrt(2) erf(&ImaginaryI;/2 sqrt(2) x)]]"

 

 

Download odesteps_fail_nm_May_10_2024_ac.mw

One thing that it shows is that the two results of using alone either simplify or combine coincide, since the set of two identical results has collapsed to a set with only one member.

{simplify,combine}(piecewise(x <= 0, 2*ln(1 - x), 0 < x, ln((x - 1)^2))) assuming (x, real);

{ln((x-1)^2)}

[simplify,combine](piecewise(x <= 0, 2*ln(1 - x), 0 < x, ln((x - 1)^2))) assuming (x, real);

 

[ln((x-1)^2), ln((x-1)^2)]

 

The sentence pertaining to this example in the Maple 2024 Advanced Math What's New page is [emphasis mine],

   "In Maple 2024, simplify and combine are both now better at
    recognizing when two piecewise branch values can be merged into one:"

I think that the main gist of that sentence is that simplify and combine are each separately ("alone") better at dealing with this example.

This is not a bug. Rather, it's a frequently asked question.

The sum command relies on Maple's standard evaluation rules for procedure calls, in which all the arguments are evaluated up front. But before p takes on actual integer values the first argument being passed to sum is actually zero.

In contrast, the add command has special evaluation rules, in which case evaluation of its first argument is delayed until the summing-index name p takes on its actual values.

restart;

e[1] := (t,x,y,z)-> vector(4,[1,0,0,0]):
e[2] := (t,x,y,z)-> vector(4,[0,1,0,-6*z^4*(1-z)^3*(1-x^2)^2*(1-y^2)^3*x]):
e[3] := (t,x,y,z)-> vector(4,[0,0,1,-6*z^4*(1-z)^3*(1-x^2)^3*(1-y^2)^2*y]):
e[4] := (t,x,y,z)-> vector(4,[0,0,0,1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3]):

lie := (j,u,v)-> u[1]*diff(v[j],t)-v[1]*diff(u[j],t)+u[2]*diff(v[j],x)-v[2]*diff(u[j],x)
+u[3]*diff(v[j],y)-v[3]*diff(u[j],y)+u[4]*diff(v[j],z)-v[4]*diff(u[j],z):

G:=(t,x,y,z)->matrix(4,4,[
[1,0,0,0],
[0,1+36*z^8*(1-z)^6*(1-x^2)^4*(1-y^2)^6*x^2,36*z^8*(1-z)^6*(1-x^2)^5*(1-y^2)^5*x*y,-6*z^4*(1-z)^3*(1-x^2)^2*(1-y^2)^3*x*(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3)],
[0,36*z^8*(1-z)^6*(1-x^2)^5*(1-y^2)^5*x*y,1+36*z^8*(1-z)^6*(1-x^2)^6*(1-y^2)^4*y^2,-6*z^4*(1-z)^3*(1-x^2)^3*(1-y^2)^2*y*(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3)],
[0,-6*z^4*(1-z)^3*(1-x^2)^2*(1-y^2)^3*x*(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3),-6*z^4*(1-z)^3*(1-x^2)^3*(1-y^2)^2*y*(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3),(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3)^2]
]):
'G(t,x,y,z)'=G(t,x,y,z):

cf:=(k,l,m,t,x,y,z) -> seq(lie(p,e[k](t,x,y,z),e[l](t,x,y,z))*G(t,x,y,z)[p,m],p=1..4):

Warning, (in cf) `p` is implicitly declared local

cf1:=(k,l,m,t,x,y,z) -> sum(lie(p,e[k](t,x,y,z),e[l](t,x,y,z))*G(t,x,y,z)[p,m],p=1..4):

cf1(2,4,2,t,x,y,z);

0

cf2:=(k,l,m,t,x,y,z) -> add(lie(p,e[k](t,x,y,z),e[l](t,x,y,z))*G(t,x,y,z)[p,m],p=1..4):

Warning, (in cf2) `p` is implicitly declared local

cf2(2,4,2,t,x,y,z);

-6*(-6*z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^2*x-6*z^4*(1-z)^3*(-x^2+1)^2*(-y^2+1)^3*x*(-3*z^2*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3-7*z^3*(z-1)^2*(y^2-1)^3*(x^2-1)^3-2*z^3*(7*z-4)*(z-1)*(y^2-1)^3*(x^2-1)^3)-(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3)*(-24*z^3*(1-z)^3*(-x^2+1)^2*(-y^2+1)^3*x+18*z^4*(1-z)^2*(-x^2+1)^2*(-y^2+1)^3*x))*z^4*(1-z)^3*(-x^2+1)^2*(-y^2+1)^3*x*(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3)


This is what gets passed to sum, which does not have special evaluation rules.

lie(p,e[2](t,x,y,z),e[4](t,x,y,z))*G(t,x,y,z)[p,2]

0

cf3:=(k,l,m,t,x,y,z) -> sum('lie'(p,e[k](t,x,y,z),e[l](t,x,y,z))*G(t,x,y,z)[p,m],p=1..4):

cf3(2,4,2,t,x,y,z);

-6*(-6*z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^2*x-6*z^4*(1-z)^3*(-x^2+1)^2*(-y^2+1)^3*x*(-3*z^2*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3-7*z^3*(z-1)^2*(y^2-1)^3*(x^2-1)^3-2*z^3*(7*z-4)*(z-1)*(y^2-1)^3*(x^2-1)^3)-(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3)*(-24*z^3*(1-z)^3*(-x^2+1)^2*(-y^2+1)^3*x+18*z^4*(1-z)^2*(-x^2+1)^2*(-y^2+1)^3*x))*z^4*(1-z)^3*(-x^2+1)^2*(-y^2+1)^3*x*(1-z^3*(7*z-4)*(z-1)^2*(y^2-1)^3*(x^2-1)^3)

 

Download bug_sum_ac.mw

When you use the key option the comparison `<` is done using evalb.

But evalb cannot resolve an inequality involving two items (the norms) which are not both of Maple type numeric.

The abs of the entries of your Vectors are all of Maple type numeric, so the infinity-norm works. But your other norm involves sqrt, and the ensuing radicals are not of type numeric. So comparison via evalb is inadequate.

But those radicals are of type realcons, hence applying evalf (in the key) produces corresponding floats, which are of type numeric.

So it's the very same reason as why you used is as your own check to compare 4th and 5th Vectors.

nb. When the key comparison fails an error is not thrown. But you can see the failure if using the equivalent two-element comparator instead of the key. option -- an error gets thrown,

restart:

V := [seq(LinearAlgebra:-RandomVector(2, generator=1..10), k=1..4)];

[Vector(2, {(1) = 3, (2) = 5}), Vector(2, {(1) = 8, (2) = 7}), Vector(2, {(1) = 4, (2) = 5}), Vector(2, {(1) = 6, (2) = 2})]

Sorting wrt L(+oo) norm

sort(V, key=(t -> norm(t, +infinity)));  # correct

[Vector(2, {(1) = 3, (2) = 5}), Vector(2, {(1) = 4, (2) = 5}), Vector(2, {(1) = 6, (2) = 2}), Vector(2, {(1) = 8, (2) = 7})]

map(norm,V);
andmap(type, %, numeric);

[5, 8, 5, 6]

true

Sorting wrt L(2) norm

sort(V, key=(t -> norm(t, 2))); # not correct

[Vector(2, {(1) = 3, (2) = 5}), Vector(2, {(1) = 4, (2) = 5}), Vector(2, {(1) = 8, (2) = 7}), Vector(2, {(1) = 6, (2) = 2})]

map(norm,V,2);
andmap(type, %, numeric);

[34^(1/2), 113^(1/2), 41^(1/2), 2*10^(1/2)]

false

evalb(norm(V[4], 2) < norm(V[3], 2)); # neither true nor false!

2*10^(1/2) < 41^(1/2)

is(norm(V[4], 2) < norm(V[3], 2));

true

sort(V, key=(t -> evalf(norm(t, 2)))); # correct

[Vector(2, {(1) = 3, (2) = 5}), Vector(2, {(1) = 6, (2) = 2}), Vector(2, {(1) = 4, (2) = 5}), Vector(2, {(1) = 8, (2) = 7})]

map(evalf@norm,V,2);
andmap(type, %, numeric);

[5.830951895, 10.63014581, 6.403124237, 6.324555320]

true

sort(V, (v1,v2) -> norm(v1,2) < norm(v2,2));

Error, sort: comparison function argument, (v1, v2) -> norm(v1,2) < norm(v2,2), must be a function that always returns true or false

sort(V, (v1,v2) -> evalf(norm(v1,2)) < evalf(norm(v2,2)));

[Vector(2, {(1) = 3, (2) = 5}), Vector(2, {(1) = 6, (2) = 2}), Vector(2, {(1) = 4, (2) = 5}), Vector(2, {(1) = 8, (2) = 7})]

sort(V, (v1,v2) -> is( norm(v1,2) < norm(v2,2) ));

[Vector(2, {(1) = 3, (2) = 5}), Vector(2, {(1) = 6, (2) = 2}), Vector(2, {(1) = 4, (2) = 5}), Vector(2, {(1) = 8, (2) = 7})]

PP := map(norm,V,2);
sort(PP); # done by address
sort(PP, key=(u->u));
sort(PP, key=evalf);
sort(PP, (a,b)->is(a<b));
sort(PP, (a,b)->a<b);

[34^(1/2), 113^(1/2), 41^(1/2), 2*10^(1/2)]

[34^(1/2), 41^(1/2), 113^(1/2), 2*10^(1/2)]

[34^(1/2), 41^(1/2), 113^(1/2), 2*10^(1/2)]

[34^(1/2), 2*10^(1/2), 41^(1/2), 113^(1/2)]

[34^(1/2), 2*10^(1/2), 41^(1/2), 113^(1/2)]

Error, sort: comparison function argument, (a, b) -> a < b, must be a function that always returns true or false

 

Download SortingVectors_ac.mw

You can do that as follows, using assumptions.

Doing it in two separate steps avoids difficulties in the interplay of assuming and the names within the objects.

restart;

with(Student:-MultivariateCalculus):

A := [0, 0, 0];
B := [a, 0, 0];
C := [a, a, 0];
DD := [0, a, 0];
S := [0, 0, h];
d1 := Line(A, C);
d2 := Line(S, DD);

A := [0, 0, 0]

B := [a, 0, 0]

C := [a, a, 0]

DD := [0, a, 0]

S := [0, 0, h]

Student:-MultivariateCalculus:-Line([0, 0, 0], Vector(3, {(1) = a, (2) = a, (3) = 0}), variables = [x, y, z], parameter = t, id = 1)

_m139664740801120

temp := Distance(d1, d2);

abs(a)^2*abs(h)/(2*abs(a*h)^2+abs(a)^4)^(1/2)

simplify(temp) assuming a>0 and h>0;

a*h/(a^2+2*h^2)^(1/2)

Download SMC_simplify_ex.mw

You do not need to use the symbolic option of simplify (which is dodgy in general).

You could also simplify the call directly, if using assume instead of assuming.

restart;

with(Student:-MultivariateCalculus):

assume(a>0);
assume(h>0);

interface(showassumed=0):

A := [0, 0, 0]: B := [a, 0, 0]: C := [a, a, 0]: DD := [0, a, 0]:

S := [0, 0, h]:

d1 := Line(A, C): d2 := Line(S, DD):

simplify(Distance(d1, d2));

a*h/(a^2+2*h^2)^(1/2)

Download SMC_simplify_ex2.mw

You could construct them separately, with the disk having no edge.

restart

with(plots)

with(plottools)

display(disk([1, 1], 5, color = red, style = polygon), circle([1, 1], 5, linestyle = dash, color = blue))

NULL

Download how_do_i_change_the_color_of_the_line_around_a_disk_while_using_plottools_ac.mw

You can also compare the effect of having the inner color extend exactly as far as that of the edge (circle).

restart

with(plots)

with(plottools)

setoptions(size = [400, 400])

display(disk([1, 1], 5, color = pink, style = polygon), circle([1, 1], 5, linestyle = dash, color = blue, thickness = 3))

display(disk([1, 1], 5, color = pink, style = polygon), circle([1, 1], 5, color = pink, thickness = 3), circle([1, 1], 5, linestyle = dash, color = blue, thickness = 3))

 

Download disk_edge_2.mw

If one simply shows the plots generated by your code then it doesn't appear that the numeric solutions for u(y) correspond to the purported "Analytic solution".  So, what is that "Analytic solution" formula supposed to represent? How did you obtain it? Be specific and clear.

Here's an adjustment to get results from your printf attempts. I adjusted the code to keep the dsolve,numeric solutions generated in its loop, so that these could be re-used when constructing the tabular data.

You haven't told us what you mean by "error". Do you mean relative absolute error, or absolute error, or something else? Be specific and clear.

HelpAhsan_ac.mw

First 22 23 24 25 26 27 28 Last Page 24 of 336