acer

32343 Reputation

29 Badges

19 years, 327 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

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

You could try renaming your .mw worksheet with a file-name that doesn't include special characters (`+`, "(), etc). Try it with only alpha-numeric characters, underscores, and a period.

Some special characters can confuse the backend of this forum, preventing downloads.

@Carl Love Does your Windows Task Manager show a similar CPU load if the option engine=triade is used instead?

@MaPal93 It is counter-productive and unnecessarily difficult for the readership if a derivative followup query is separated from its parent, since that obscures connections amongst the information and details given by the original poster as well as the responders.

@Carl Love Despite the wording of the error message I don't see that the sign command is necessarily the problem here.

I'd guess that signum is being used at some internal juncture, and the calling routine is not properly prepared to deal with its result...

@MaPal93 Please do not put that followup in a (duplicate) separate Question thread.

How would you feel about taking limits instead of using eval? Or do you want to know precisely how dsolve's internals are handling it differently?

note: in Maple 2017.3 (and 2022.1) a single call to limit can do the same.

restart

kernelopts(version);

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

sys_1 := {diff(x(t), t$2)=sin(t)-x(t), x(0)=0, D(x)(0)=0}:

sol_1 := dsolve(sys_1);

x(t) = (1/2)*sin(t)-(1/2)*cos(t)*t

sys_2 := {(M__p+M__a)*diff(x(t), t$2)=M__p*sin(t)-x(t), x(0)=0, D(x)(0)=0}:
sol_2 := dsolve(sys_2)

x(t) = sin(t/(M__p+M__a)^(1/2))*M__p*(M__p+M__a)^(1/2)/(M__p+M__a-1)-M__p*sin(t)/(M__p+M__a-1)

limit(eval(sol_2, M__p=1), M__a=0);

x(t) = (1/2)*sin(t)-(1/2)*cos(t)*t

limit(eval(sol_2, [M__a=0]), M__p=1);

x(t) = (1/2)*sin(t)-(1/2)*cos(t)*t

Download SomethingWrong_ac.mw

Shout if you have reservations.

@jalal Here is my suggestion for 2D math in the "headers" together with Carl's suggestion for additional substitution of complex values.

DocumentTools:-Tabulate(subsindets[flat](
                          Matrix([<["",F(x)[]]>,
                                  ((Vector@[u->u,F[]])~(C))[]]),
                          {undefined,And(complexcons,Not(realcons))},
                          ()->dne),
                        width=50,
                        fillcolor=((T,i,jj)->`if`(i=1 or jj=1,
                                                  [210,230,255],
                                                  [255,255,235]))):

You could put that into a reusable procedure, if you'd like, which could be hidden in a document's Startup code or a collapsed Code-Edit-Region. Such a procedure could simply take F and C as its arguments.

@jalal You could instead use a Matrix for the whole thing, with data and headers together, to get the header expressions to render in pretty-printed 2-D math.

DocumentTools:-Tabulate(subs(undefined=dne,
                             Matrix([<["",F(x)[]]>,
                                     ((Vector@[u->u,F[]])~(C))[]])),
                        width=50,
                        fillcolor=((T,i,jj)->`if`(i=1 or jj=1,
                                                  [210,230,255],
                                                  [255,255,235]))):

That summation index variable n1 is not entirely like the variable n. That name n1 represents a bound variable, since your sum call denotes that it would only take integer values between 1 and n-1.

There's not much significant difference between using the name n1 instead of the name i, or what have you. For example, these next two are conceptually similar:

    sum(f(i), i = 1 .. n-1)
    sum(f(n1), n1 = 1 .. n-1)

The end result of both of those would usually be the same; the results depend on (any value of) n in the upper end-point of the index, but the results don't depend on which choice of name is used for the summation index.

Maple knows a bit of that, which can be useful:

depends(sum(f(n1),n1=1..n-1), n);

         true

depends(sum(f(n1),n1=1..n-1), n1);

         false

depends(sum(f(i),i=1..n-1), i);

         false

expr := sum(f(i),i=1..n-1):

indets(expr, And(name,
                 satisfies(nm->depends(expr,nm))));

          {n}

@David Mazziotti Maple's GUI renders a certain kind of markup for its 2D Math pretty-printing, a.k.a. typesetting. There is a strong resemblance to MathML. The Maple specifics are not user-documented.

The "mo" stands for math-operator, and "mi" for math-identifier. (See here, for similar details.)

It can come in two forms. One construction consists of specially formed names, and the other of special function calls using exports of the Typesetting module. For example,

This first is all one long name.

`#mrow(mi("a",mathcolor="FF00BB"),msub(mi("x"),mi("i")))`;

`#mrow(mi("a",mathcolor="FF00BB"),msub(mi("x"),mi("i")))`

This next consists of nested function calls.

Typesetting:-mrow(mi("a",mathcolor="FF00BB"),
                  Typesetting:-msub(Typesetting:-mi("x"),
                                    Typesetting:-mi("i")));

a*x[i]

Download typemkex.mw

It can be very handy to be able to programmatically construct things which print in special ways (including the case where there is no usual Maple structure that would normally print in that manner).

Quite often one can infer how these forms work by creating some 2D Input, then using the right-click context-menu item,
   2-D Math -> Convert To -> Atomic Variable
and then executing, and then calling lprint on the result to show its structure.

See also the results of,
   exports(Typesetting);
which are mostly undocumented.

@MalakMMK I don't think that you coded the formula from the Hafiz paper correctly; at least some brackets are missing.

For example you had, as the first step,
   u - (2*F(u)/3*dF(u))
but the paper shows the denominator of that as all of 3*dF(u), hence,
   u - 2*F(u)/(3*dF(u))
Also, you had,
   v - F(u)*W/dF(u)+dF(v)
but the paper show the denominator of that as all of dF(u)+dF(v), hence,
   v - F(u)*W/(dF(u)+dF(v))

The paper defines "eta = f'(x)/f'(y) without the index n", which seems a little confusing. Below I used eta := dF(u)/dF(v) and get apparent success. Please edit and explain, if you think that's wrong.

[edit: Below I used a value of 20 for the maximal iteration cutoff, ie the third argument passed to H. So for a sensible comparison with the earlier Newton's and King's methods you might want to re-run those with that same value.]

Haf_Newton_ac1.mw

 

First 98 99 100 101 102 103 104 Last Page 100 of 592