acer

26592 Reputation

29 Badges

17 years, 0 days

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

With the H I supplied then --  given the form of f -- the simplest check is likely the following (Maple 2022.1),

    simplify(diff(H,x) - f);
                      0

For the F you gave, one slightly terser check is,

    (simplify@@2)(expand(diff(F,x) - f))

My earlier concoction manipulates diff(H,x) in order to obtain f, which (as is quite often the case) is a more involved task that simplifying their difference to zero.

The integration should be fun to see.

@MaPal93 That old posting you cited is (IMO) not the best way to install the DirectSearch v.2 package, in modern Maple versions. Here are some alternatives, in a descending priority that I would prefer:

1) Open the Maple cloud directly from your Maple GUI, and install that package directly from there.
2) Visit the Maple cloud webpage for that package and download from there.
3) Use the installer that comes in the .zip file bundled at the Application Center page for that package.

All of those approaches should put the relevent files into some special location that Maple looks for (so that you don't have to manually adjust libname in an edited initialization file).

You may have to completely close the Maple GUI and relaunch in order to see the Help pages that come with it.

If none of those approaches work then it would be interesting to know what is in the folder/directory that'd be named as the output of this Maple command:
    cat(kernelopts(homedir),"/maple/toolbox");

I have not yet found a satisfactory way.

I got something backwards (a simplification challenge):

H := 1/2*(-arctan((cos(x) - sin(x))/sqrt(2 + sin(2*x)))
     - log(cos(x) + sin(x) + sqrt(2 + sin(2*x))))

-(1/2)*arctan((cos(x)-sin(x))/(2+sin(2*x))^(1/2))-(1/2)*ln(cos(x)+sin(x)+(2+sin(2*x))^(1/2))

combine(simplify(expand(rationalize(diff(H,x)))));

sin(x)/(2+sin(2*x))^(1/2)

Download backwds.mw

@C_R The second argument to the indets command specifies a type, and that is documented. That is why Tom's answer doesn't use the word "type" explicitly-- it's implied by that particular choice of command.

And type suffixed is documented.

note: The top-level (global) name suffixed is not protected (Maple 2022.1) and so could be assigned some value. Some modest defensive programming could avoid accidental failure in such a case, eg. by using unevaluation quotes as Carl did.

   indets(A, 'suffixed'(t, :-nonnegint))


By the way, if you wish to use select here then you can get by without the extra layers of a user-defined anonymous operator (arrow) and an `if` call, both shown in an earlier suggestion. You had,

   select(X -> `if`(type(X, 'suffixed(t, integer)'), true, false), lst)

The `if`(...,true,false) layer is not needed because type itself returns true or false. And the anonymous operator is not needed because any additional arguments (3rd onward) passed to select will get passed along to the 1st argument (as predicate). So you can get by with the terser and more efficient,

   select(type, lst, 'suffixed'(t,:-nonnegint))

That kind of use of (only) kernel builtins can bring a small absolute time savings on average. But the relative time savings can be large, and so it can often scale better to large examples. So it can be worthwhile to program in such a way when possible.

@nm Your odetest example induces the following call to solve, whose result is different on the second attempt (with Carl's redefinition in play...).

restart;

kernelopts(version);

`Maple 2022.1, X86 64 LINUX, May 26 2022, Build ID 1619613`

#from https://www.mapleprimes.com/questions/230035-Error-in-SolveToolsCancelInverses
CI__orig:= eval(SolveTools:-CancelInverses):
unprotect(SolveTools:-CancelInverses):
SolveTools:-CancelInverses:= e->
local _A;
    eval(CI__orig(eval(e, _Z= _A), _rest), _A= _Z)
:
protect(SolveTools:-CancelInverses, CI__orig):

P := (-42*((-2*x^2+3*x*y(x)+y(x)^2)/x^2)^(1/2*17^(1/2))*2^(-3/2*17^(1/2))*((-11/21*17^(1/2)-17/7)*x^(2*17^(1/2)+1)+y(x)*x^(2*17^(1/2))*(17^(1/2)+85/21))+4*(1/8*(-2*x^2+3*x*y(x)+y(x)^2)/x^2)^(-1/2*17^(1/2))*((2*y(x)+3*x)*17^(1/2)+17*x))/(-5/4*x^(2*17^(1/2))*2^(-3/2*17^(1/2))*((-2*x^2+3*x*y(x)+y(x)^2)/x^2)^(1/2*17^(1/2))*(17^(1/2)+21/5)+(1/8*(-2*x^2+3*x*y(x)+y(x)^2)/x^2)^(-1/2*17^(1/2))),
y(x):

a1 := solve(P): # first call

a2 := solve(P): # second call

nops([a1]);

1

nops([a2]);
a2[2];

2

-(1/2)*(17^(1/2)+3)*x

# The above call to `solve` happens here, for the given odetest example.
#
showstat(`ODEtools/Solve/EnvDropMultiplicity`,29);


`ODEtools/Solve/EnvDropMultiplicity` := proc(ee, X := NULL, {keepalreadysolveduntouched::truefalse := false, removelabel::truefalse := false})
local EE, sol, N, sys, zn, vars, sysaux, tmp, to_jet, j, k, _sum, Q, Q_eqs, P,
  copy_P, PL, P2, PL2, PL21, PQ, ee_signum, sol_signum;
       ...
  29                   sol := :-solve(ee,X,_rest)
       ...
end proc
 

Download cancelinverses_rum_ex.mw

@MalakMMK Here is something like the previous simple (manually written) looped code. Of course it now has to store the converged value since that's the kind of image you requested.

So I adjusted the looped code to store the converged values as a side-effect in a new Array. That Array is then used for the hue values in the image, according to classification according to nearest actual root. The iteration counts are also (still) returned, and those are extracted (in bulk) from the plot3d result and used for the intensity values in the image.

The classification by root takes slightly longer than the iterative process (which is slower in this looped code than it is using Escape, contrary to your claim).

Of course one could improve the performance of the code that classifies by root. I don't have the time for that. Perhaps you could do it; this is your project work, after all.

King_iter_Ma2_umm.mw

[edit: there are also examples for which symmetry can provide additional savings.]

@MalakMMK Do you understand than the two images compute two different things, which is one reason why they differ in speed (as well as the main reason why they different in coloring)?

The first image has four colors because it shades according to which root, and your current example has four roots. It takes some extra time to classify the points -- by root.

You wrote, "...in the second image I'm Always getting 3 colors" but that is not true. It has a variety of colors, whose hue corresponds to the number of iterations counted.

In my earlier code I tried to simplify the real and imaginary parts of the iterative formula. That simplification can be too expensive, and even unnecesary as my later code showed for the Hafiz method.

You should try to understand the earlier code to the degree that you could verify that it's working properly. (I might have made mistakes!) And if you cannot edit and adjust the code then how can you adequately understand it? If you cannot adequately understand the code then how can you use it in your coursework?

@MalakMMK In your most recent followup query (ie. the Repley immediately above) you have not made it completely clear what you want.

You have not provided the code you used to generate those images, which doesn't help with any guessing on out part.

You do not have to inline all of any attachments, here. You can upload it and then insert a link (instead of the whole contents, inlined).

I think that you might be trying to ask for something like this:
  - utilize the looped code (as used by my densityplot examples)
  - instead of producing the HUE-based densityplot image (colored by number of iterations to attain tolerance), produce the other kind of image I previously showed where the color is according to which root is closest to the converged value.

If that guess by me is correct then it could be done by extracting the hue data from the densityplot result and then treating it similarly to this line in my earlier procedure:

   # ...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):

I may have some time to look at that tomorrow evening, or possibly the day after.

ps. I don't quite understand why you might prefer the manually written looped code. That Escape command provides a much faster (multithreaded, then compiled or evalhf'd) implementation of the loop's same iterate-while-counting-until-tolerance-met idea. The faster way worked fine for the Hafiz 4th order method.

@Stretto You wrote, "printing has nothing to do with mouse hovering".

I can think of one mechanism where a displayed list of point data could cause arbitrary extra information to be revealed upon mouse hover, another for mouse single left-click, another for mouse right-click, and another for mouse left-button-selection. None of those need show positional "graphs" of the data.

@PhearunSeng 

Tom's worksheet has eqs as a list.

Your image shows that you rewrote that so that eqs is a Matrix.

You claimed that your image shows you running Tom's worksheet. But your image shows that the command's that you are using are different from his.

Download Tom's attachment. Open that .mw file, and run it. Don't paste in something else.

Also, it makes little sense to have the equations be some scalars expressions added to Vectors. Be more careful about your syntax.

@PhearunSeng You are still calling restart after you define eqs.

You wrote, "Is there a way to print points that have..." (italics mine). But then you talk of "graphs".

So, are you trying to plot the point data, or print it like output, or...?

And you want hover-over information to popup, for different points?

(You description of how you're reading in the data, and the slowness of that, seems to be a quite different query. But it's difficult to offer best advice when you don't show code-to-reproduce.)

@Carl Love The OP sent me a link to another paper.

https://www.sciencedirect.com/science/article/pii/S0377042708003555

He said, "I used equation 14, Method 1 , and took values for alpha and gamma = 1".

And he sent me a worksheet that had an implementation of the formulas. As a rough test of that I plotted its "escape-time" using densityplot.

restart:

ee := z^3-1:

__FF := proc(z) option inline; z^3-1 end proc:
Wei1Count := subs(g=3, F=__FF, dF=D(__FF),
  proc(x0,y0) local u,i,v,A,p,fpu,fpv,fpuu,ut,W;
  u := evalf(x0+y0*I);
  for i from 1 to 20 do
    v := u - F(u)/dF(u);
    A := F(u)+(g-2)*F(v);
    p := v - 1/A*(F(u)+g*F(v))*F(v)/dF(u);
    fpu := (F(u)-F(p))/(u-p);
    fpv := (F(v)-F(p))/(v-p);
    fpuu := (dF(u)-fpu)/(u-p);
    ut := F(p)/F(u);
    W:= 1+2*ut+ut^2+ut^3;
    u := p - F(p)*W/(fpv +fpuu*(p-v));
    if abs(F(u))<1e-3 then return i; end if;
  end do:
  return i-1;
end proc):

 

plots:-densityplot(Wei1Count, -2..2, -2..2, grid=[201,201],
                   title="Wei's method (1)", size=[300,300],
                   colorstyle=HUE, style=surface);

NULL

Download Wei_method_ac2.mw

I haven't time to figure out why my earlier EscapeTime template runs away on it...

[note: in my original answer here I preprocessed the Newton's method (compounded) formula and hit its Im and Re components with evalc and simplify. I happened to do the same for my King's method implementation in that same Answer. But I did not do that in my followup implementation of the "Haf1" (Hafiz paper) 4th order method. There I left the Im and Re calls alone, since the simplify/evalc time is large, it's of little merit, and these days evalhf and the Compiler can handle those better.]

@Dkunb Here are a couple of ways to get that first image -- the contour N(beta,alpha)=0.

I used Axel's idea of integrating only out to 6 instead of infinity, and using oscillating method _d01akc instead of semi-infinite method _d01amc.

I used a list of values by rootfinding (in two ways), instead of contourplot or implicitplot.

I also show a simple (simplistic) idea of using Interpolation and using a coarser list of points. But the final curve appears quite smooth.

I also tried Grid:-Seq, and you might try tweaking that.

I may not have any more spare time for looking at this interesting problem. But this might illustrate a few alternate approaches.

restart;

F:=kappa->kappa:
f:=(alpha,delta)->exp(-abs(F(kappa))^2*(1+delta^2)/2-abs(F(kappa))*alpha)/abs(F(kappa)):

L:=(alpha,delta,Lambda) ->
  (lambda^2*exp(-alpha^2/2)/4)*(Int(f(alpha,delta),kappa= -infinity..-Lambda)+Int(f(alpha,delta),kappa= Lambda..infinity)):

CodeTools:-Usage(evalf(L(4,1,0.001)));

memory used=1.88MiB, alloc change=0 bytes, cpu time=65.00ms, real time=65.00ms, gc time=0ns

0.8209373770e-3*lambda^2

X:=Int(exp(-kappa^2*(delta^2 + 1)/2)*cos(kappa*beta)*cos(kappa*gg)/(kappa),
       kappa = Lambda .. infinity):

Y:='Int(-exp(-kappa^2*(delta^2 + 1)/2)*cos(kappa*beta)*Im(exp(kappa*gg*I)*erf( (gg+kappa*I)/sqrt(2)))/(kappa), kappa = Lambda .. infinity)':

theJ:='sqrt(X^2+Y^2)/2*exp(-alpha^2/2)*lambda^2':

J:=unapply(subs(infinity=6, theJ), [alpha,delta,Lambda,beta,gg]): # cut off

thisJ:=unapply(subsindets(J(alpha,1,0.001,beta,3),
                          'specfunc(anything,Int)',
                          'q -> Int(op(q), digits=15, epsilon = 1e-4, method = _d01akc)'),[alpha,beta]):

evalf(thisJ(4,8));

0.7304273940e-3*lambda^2

N:=proc(beta,alpha) option numeric; local res;
     Digits:=15;
     if not [beta,alpha]::[numeric,numeric] then
     return 'procname'(args); end if;
     res:=thisJ(alpha,beta) - L(alpha,1,0.001);
     res:=evalf[15](res);
     res:=evalf[10](res);
     if type(res/lambda^2,numeric) then
     res/lambda^2; else undefined; end if;
end proc:

forget(evalf):
CodeTools:-Usage( evalf(N(4,8)) );

memory used=13.14MiB, alloc change=4.00MiB, cpu time=89.00ms, real time=90.00ms, gc time=0ns

0.7385644022e-14

#forget(evalf);
#CodeTools:-Usage( plots:-contourplot(N(beta,alpha), beta=0..10, alpha=0..2, grid=[7,7]) );

#forget(evalf);
#CodeTools:-Usage( plots:-contourplot(N(beta,alpha), beta=0..10, alpha=0..10,
#                                     contours=[-1e-8, 1e-8], grid=[17,17]) );

N:=subsop(4=NULL,eval(N)): # clear N's remember table
forget(evalf);
time[real](assign('PP',[seq([beta,fsolve(alpha->N(beta,alpha),1..10)],
                            beta=0..10,1.0)]))*`seconds real time`;

108.788*`seconds real time`

#eval(PP,1);

PPSPL:=Interpolation:-Interpolate(PP[..,1],PP[..,2]):

plot(PPSPL, 0..10, view=[0..10,0..10],
     filled=true, color=gray, labels=[beta,alpha], title=(gamma=3),
     transparency=0, background="aa0000");

N:=subsop(4=NULL,eval(N)): # clear N's remember table
forget(evalf);
time[real](assign('PNZ',
                  [seq([bval,
                        RootFinding:-NextZero(subs(beta=bval,
                                                   alpha->N(beta,alpha)),
                                              1.0)],bval=0..10,1.0)]) )*`seconds real time`;

106.904*`seconds real time`

#eval(PNZ,1);

PNZSPL:=Interpolation:-Interpolate(PNZ[..,1],PNZ[..,2]):

plot(PNZSPL, 0..10, view=[0..10,0..10],
     filled=true, color=gray, labels=[beta,alpha], title=(gamma=3),
     transparency=0, background="aa0000");

N:=subsop(4=NULL,eval(N)):
forget(evalf);
Grid:-Set(N);Grid:-Set(L);Grid:-Set(thisJ);Grid:-Set(J);Grid:-Set(f);Grid:-Set(F);
time[real](assign('PGR',
                  [Grid:-Seq['tasksize'=max(1,floor(10/kernelopts('numcpus')))]([beta,
                                                                                 fsolve(N(beta,alpha),
                                                                                        alpha=1..10)],
                                                       beta=0..10,1.0)]) )*`seconds real time`;

57.881*`seconds real time`

#eval(PGR,1);

PGRSPL:=Interpolation:-Interpolate(PGR[..,1],PGR[..,2]):

plot(PGRSPL,0..10, view=[0..10,0..10],
     filled=true, color=gray, labels=[beta,alpha], title=(gamma=3),
     transparency=0, background="aa0000");

 

Download Negativity_v1_acU0.mw

Please don't create duplicate Question threads of your followup (below) about an 8th order method.

(I hope to find time to look at your formulation. But the fact that nobody's yet done so here is not justification to submit duplicates.)

1 2 3 4 5 6 7 Last Page 2 of 494