acer

32368 Reputation

29 Badges

19 years, 333 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

I've tried to adjust the methodology to get better performance from the numerical integration computations of eq7 & eq8 after the substitutions.

Getting that fast (for individual values of alpha and beta) is, naturally, key to figuring out efficiently from purely numerical computations where the inequality might hold.

I'm pretty sure that the following adjustments to calling evalf/Int, etc, gets a faster computation than merely removing the unapply, as CW suggested in his Answer.

More importantly, your original pointplot approach would compute something like 26*900 times, for alpha and beta sampled in a evenly spaced grid. A fine grid sampling is often an unnecessarily expensive way, and gets only a coarse visual of an inequality.

Let me know if this below is helpful. (You were originally sampling 10^(-9) to 1=10^0, etc. I've changed the ranges a bit.)

restart;

eq1:=diff(a(t),t)/a(t)=alpha*t-beta*t^3:

eq2:=dsolve({eq1,a(0)=1}):

eq3:=diff(lambda(t),t)=rhs(eq2):

eq4:=dsolve({eq3,lambda(0)=0}):

sm:=exp(-lambda^2/(2*tau^2))*1/(sqrt(2*Pi)*tau):

eq5:=sm*diff(rhs(eq1),t)/rhs(eq2):

eq6:=combine(eval(eq5,[lambda=rhs(eq4),tau=eval(rhs(eq4),t=1)])):

eq7:=subs(Int(exp(-_z1^2*(_z1^2*beta - 2*alpha)/4), _z1 = 0 .. 1)=L,eq6):
eq7 := subsindets(eq6, specfunc(Int),
                  u->op(0,u)(op(u), digits=15, epsilon=1e-11, method=_d01ajc));

(1/2)*2^(1/2)*(-3*beta*t^2+alpha)*exp(-(1/2)*(Int(exp(-(1/4)*_z1^2*(_z1^2*beta-2*alpha)), _z1 = 0 .. t, digits = 15, epsilon = 0.1e-10, method = _d01ajc))^2/(Int(exp(-(1/4)*_z1^2*(_z1^2*beta-2*alpha)), _z1 = 0 .. 1, digits = 15, epsilon = 0.1e-10, method = _d01ajc))^2+(1/4)*t^2*(beta*t^2-2*alpha))/(Pi^(1/2)*(Int(exp(-(1/4)*_z1^2*(_z1^2*beta-2*alpha)), _z1 = 0 .. 1, digits = 15, epsilon = 0.1e-10, method = _d01ajc)))

GG:=Int(exp(-_z1^2*(_z1^2*beta - 2*alpha)/4), _z1 = 0 .. 1, digits=15, epsilon=1e-11, method=_d01ajc);

Int(exp(-(1/4)*_z1^2*(_z1^2*beta-2*alpha)), _z1 = 0 .. 1, digits = 15, epsilon = 0.1e-10, method = _d01ajc)

eq8:=1/(8*L^2);

(1/8)/L^2

HH := proc(a,b) option remember; local res;
  if not [a,b]::list(numeric) then return 'procname'(args); end if;
  res := evalf(eval(Int(eval(eval(eq7,L=GG),[:-alpha=a,:-beta=b]), t=-25..25,
                        epsilon=1e-5, method=_Dexp))
               - eval(GG,[:-alpha=a,:-beta=b]));
  if res::numeric then evalf[5](res); else undefined; end if;
end proc:

 

CodeTools:-Usage(
  plots:-inequal(HH(alpha,beta)>=0
         , alpha=-1.5..0, beta=-2..0
         , optionsimplicit=[gridrefine=1]
         )
);

memory used=14.96GiB, alloc change=148.53MiB, cpu time=3.44m, real time=3.09m, gc time=29.88s

 

Download NestedError_accc.mw

ps. Your original approach was also flawed in how it generated the pointplot data pairs. The table M was converted to a list, but the i values (indices) were being dropped, that way. If the inequality didn't hold (or the integration failed) then there'd be no value stored in M. So not all i values are even implicit in the M data. The generation of the pointplot pairs was thus all muddled up, using the wrong i values for the stored j values in M.

Purely numeric integration can be quicker, here. (The later plot done in 22ms.)

Is this result as expected?

Gaussian_integration_ac.mw

Is this the kind of thing that you're trying to accomplish?

restart;

P := expand((cosh(k*eta) + alpha*m*sinh(k*eta)/eta)^3);

cosh(k*eta)^3+3*cosh(k*eta)^2*alpha*m*sinh(k*eta)/eta+3*cosh(k*eta)*alpha^2*m^2*sinh(k*eta)^2/eta^2+alpha^3*m^3*sinh(k*eta)^3/eta^3

P := convert(P, exp);

((1/2)*exp(k*eta)+(1/2)*exp(-k*eta))^3+3*((1/2)*exp(k*eta)+(1/2)*exp(-k*eta))^2*alpha*m*((1/2)*exp(k*eta)-(1/2)*exp(-k*eta))/eta+3*((1/2)*exp(k*eta)+(1/2)*exp(-k*eta))*alpha^2*m^2*((1/2)*exp(k*eta)-(1/2)*exp(-k*eta))^2/eta^2+alpha^3*m^3*((1/2)*exp(k*eta)-(1/2)*exp(-k*eta))^3/eta^3

frontend(coeff,[P, eta, 0]);

((1/2)*exp(k*eta)+(1/2)*exp(-k*eta))^3

frontend(coeff,[P, eta, -1]);

3*((1/2)*exp(k*eta)+(1/2)*exp(-k*eta))^2*alpha*m*((1/2)*exp(k*eta)-(1/2)*exp(-k*eta))

frontend(coeff,[P, eta, -2]);

3*((1/2)*exp(k*eta)+(1/2)*exp(-k*eta))*alpha^2*m^2*((1/2)*exp(k*eta)-(1/2)*exp(-k*eta))^2

frontend(coeff,[P, eta, -3]);

alpha^3*m^3*((1/2)*exp(k*eta)-(1/2)*exp(-k*eta))^3

combine(expand(frontend(coeff,[P, eta, -3])));

(1/8)*alpha^3*m^3*exp(3*k*eta)-(3/8)*alpha^3*m^3*exp(k*eta)+(3/8)*alpha^3*m^3*exp(-k*eta)-(1/8)*alpha^3*m^3*exp(-3*k*eta)

Q := combine(collect(expand(P),eta));

(1/8)*exp(3*k*eta)+(3/8)*exp(k*eta)+(3/8)*exp(-k*eta)+(1/8)*exp(-3*k*eta)+((3/8)*alpha*m*exp(3*k*eta)+(3/8)*alpha*m*exp(k*eta)-(3/8)*alpha*m*exp(-k*eta)-(3/8)*alpha*m*exp(-3*k*eta))/eta+((3/8)*alpha^2*m^2*exp(3*k*eta)-(3/8)*alpha^2*m^2*exp(k*eta)-(3/8)*alpha^2*m^2*exp(-k*eta)+(3/8)*alpha^2*m^2*exp(-3*k*eta))/eta^2+((1/8)*alpha^3*m^3*exp(3*k*eta)-(3/8)*alpha^3*m^3*exp(k*eta)+(3/8)*alpha^3*m^3*exp(-k*eta)-(1/8)*alpha^3*m^3*exp(-3*k*eta))/eta^3

frontend(coeff,[Q, eta, -3]);

(1/8)*alpha^3*m^3*exp(3*k*eta)-(3/8)*alpha^3*m^3*exp(k*eta)+(3/8)*alpha^3*m^3*exp(-k*eta)-(1/8)*alpha^3*m^3*exp(-3*k*eta)

Download coeff_ex.mw

I'm not sure what collect has to do with this.

ps. The name eta also appears in the trigh/exp calls, eg. sinh(k*eta). Is it possible that you might want/need some kind of (truncated) series in eta?

Using active product instead of inert Product, the following can be done.

I think that it needs Maple 2018 or later, though, and not the Maple 2015 you used.

restart;

f := Product(x[i]^a*(1-x[i])^b, i);

Product(x[i]^a*(1-x[i])^b, i)

Lf := ln(f);

ln(Product(x[i]^a*(1-x[i])^b, i))

expand(simplify(expand(value(Lf))))
  assuming x[i] > 0, x[i] < 1, a > 0, b > 0;

a*(sum(ln(x[i]), i))+b*(sum(ln(1-x[i]), i))

Download From2to7_ac.mw

Is this the kind of thing that you want to accomplish?

restart

V := Vector(2, [25.0*Pi*r/(r^2+z^2), 25.0*Pi*z/(r^2+z^2)])

Vector[column](%id = 36893628704703683148)

LinearAlgebra:-Norm(V, 2)

(6168.502752*abs(r/(r^2+z^2))^2+6168.502752*abs(z/(r^2+z^2))^2)^(1/2)

 

I got the next 2D Input by using the double-bar entry
in the Common Symbols palette, and then Ctl-Shift-underscore 2
to get the subscript for the 2-norm.

 

LinearAlgebra[Norm](V, 2)

(6168.502752*abs(r/(r^2+z^2))^2+6168.502752*abs(z/(r^2+z^2))^2)^(1/2)


I got the next 2D Input by typing Norm in 2D Input mode, and

then the Esc key to get the command-completion templates
for that. I chose the second item in the ensuing popup, for the
2-norm.

LinearAlgebra[Norm](V, 2)

(6168.502752*abs(r/(r^2+z^2))^2+6168.502752*abs(z/(r^2+z^2))^2)^(1/2)

Download 2nrm.mw

ps. If you've loaded the LinearAlgebra package then you could simply call Norm(V,2) .

The plots:-display command can also handle a sequence of plots.

And the seq command can turn your Vector of plots into a sequence of plots.

restart

with(Physics); with(LinearAlgebra); with(VectorCalculus); with(Optimization); with(Statistics); with(ArrayTools); with(plots); with(plottools); with(Threads); with(MmaTranslator[Mma]), with(StringTools); with(CodeGeneration); with(ImageTools); with(ImageTools:-Draw); VectorCalculus:-`*`(Setup(mathematicalnotation = true), Setup(coordinatesystems = cartesian))

`Default differentiation variables for d_, D_ and dAlembertian are:`*{X = (x, y, z, t)}

`Systems of spacetime Coordinates are:`*{X = (x, y, z, t)}

radius := 1.0

alpha := sin((1/4)*Pi)

step := .2

radius := 100.0

l1 := line([0, 0], [100, 100])

loAng := 0.

hiAng := (1/2)*Pi

c1 := arc([0, 0.], radius, loAng .. Pi, color = "red")

c2 := arc([step, 0.], -alpha*step+radius, loAng .. hiAng, color = "blue")

c3 := arc([2*step, 0.], -2*alpha*step+radius, loAng .. Pi, color = "purple")

c4 := arc([3*step, 0.], -3*alpha*step+radius, loAng .. Pi, color = "black")

`cir&Assign;arc`([x, 0.], radius-step, loAng .. hiAng, i, color = "red")

xLimitWall := 500.0; nScatter := 20; step := xLimitWall/(nScatter+1); x := 0.; for i from 0 to nScatter do x := x+step; cir := arc([x, 0.], radius-step, 0 .. Pi, color = "red"); if i = 0 then lc := Vector([cir]) else i; cir; Append(lc, cir) end if end do

plots:-display(seq(lc))

Download JFKWEdgeDifractionDirection_ac.mw

You could also turn your Vector into a list and pass that instead, eg. you could do one of these:
    plots:-display([seq(lc)]);
or,
    plots:-display(convert(lc,list));

(While repeatedly appending to a Vector is not the most efficient way to construct your Vector of plots, I didn't change that aspect.)

ps. Your attached worksheet was last saved in Maple 2018, and not Maple 18 (released: year 2014).

Is this the kind of thing you're trying to accomplish?

fcn1 := sum(N[i]*z[i]^2, i);

sum(N[i]*z[i]^2, i)


The terms N__i and z__i (typed with double-underscore) are just whole
names, and don't depend on i.

They're not what you want, since they are unrelated to the sum index i.

fcn2 := sum(N__i*z__i, i);

N__i*z__i*i


Is this the kind of thing you want to accomplish?

diff(fcn1, N[i]);

sum(z[i]^2, i)


In 2D Input mode, you can enter N[i] by using the keystrokes
    N Ctl-Shift-underscore i

diff(fcn1, N[i])

sum(z[i]^2, i)

Download subsc_diff_ex.mw

The unapply command is what you need here.

It has been available in every Maple version for decades.

Below, I use the unapply command to construct a re-usable operator (procedure) from the output expression than was assigned to fin.

restart

with(Statistics)

 

X := [0, 8, 10, 14, 25, 28]

Y := [50, 250, 300, 400, 700, 775]

 

fin := Fit(a*x+b, X, Y, x)

HFloat(26.09176748303333)*x+HFloat(42.86662732369449)

f := unapply(fin, x)

f(25)

HFloat(695.1608143995278)


Note that you don't have to assign to fin, since
you could also act directly on Fit call.

f2 := unapply(Fit(a*x+b, X, Y, x), x)

f2(25)

HFloat(695.1608143995278)


If you do assign the expression to fin, you could
also use the eval command to evaluate the
expression at values of x.

 

eval(fin, x = 25)

HFloat(695.1608143995278)

map(f, X)

[HFloat(42.86662732369449), HFloat(251.60076718796114), HFloat(303.7843021540278), HFloat(408.15137208616113), HFloat(695.1608143995278), HFloat(773.4361168486278)]

[seq(fin, x = X)]

[HFloat(42.86662732369449), HFloat(251.60076718796114), HFloat(303.7843021540278), HFloat(408.15137208616113), HFloat(695.1608143995278), HFloat(773.4361168486278)]

Download unapply_ex.mw

Instead of trying to reset that numbering you should improve your program so that it is robust with respect to the naming of such introduced variables.

Let's suppose that you have an expression assigned to expr, which might be a solution from solve (or a list of set of such, etc). Then you can issue,

    indets(expr, suffixed(_Z,posint))

and get a set of such names present in the expression.

You can then programatically substitute for the introduced name(s) into the expression, without having to hard-code the name.

For example,

expr := solve(cos(x*Pi),x,allsolutions);

        expr := 1/2 + _Z5~

nms := indets(expr, suffixed(_Z,posint));

          nms := {_Z5~}

eval(expr, nms[1]=6);;

              13/2

seq(eval(expr, nms[1]=v), v=1..6);;

      3/2, 5/2, 7/2, 9/2, 11/2, 13/2

Notes:

1.) The need for robust handling of indeterminate names is standard stuff.

2.) Resetting that counter isn't a great idea in general. Different solutions with potentially differently assumed introduced-names should have such names be distinct, regardless of whether the assumptions happen to be. (I'm deliberately going to not mention how to do it.)

Based on your description so far, I'd suggest using a PlotComponent within your document, to get that aspect of interactivity.

If you read the Help page for that embedded component you'll see that the Component Properties (clickx, clicky) can provide you with access to the position of the mouse cursor when you Click (or Click&Drag) on a 2D plot frame. Similarly for the properties endx and endy, when a Drag ends, or hoverx and hovery when the pointer hovers over.

That access is also conveniently programmatic, and fits in well with triggering/working with other interactivity. The x- and y-values can be easily used for other computations, in a manner straightforwardly consistent with other interactive purposes (and documented for such).

The plot `annotation` mechanism is more convenient (IMO) for displaying the x- and y-values in a formatted "hover-over" frame over the plot (though it's possible it might be used programmatically to trigger an embedded component action elsewhere). It has more functionality related to closeness of plotted features in the plotting frame.

I'm not able to tell from your Question and details so far which is closer to your situation and wider goals.

Using space " " is not reliable here, in my experience.

In the past, I've resorted to using one of:
    &#8201;      thin space
    &#8194;      en space
    &#8195;      em space


IIRC, I was not able to get those as named HTML entities, eg.   &emsp;    

Here are three ways. (There are others.)

restart

y := x+1

Explore(InertForm:-Display(%y, inert = false) = y, x = 0 .. 5, echoexpression = false)

Explore(`#mi("y");` = y, x = 0 .. 5, echoexpression = false)

F := proc (X) options operator, arrow; 'y' = eval(y, x = X) end proc; Explore(F(x), x = 0 .. 5, echoexpression = false)

Explore_with_Inert_Form_of_display_ac.mw

ps. I suspect that the uneval-quoting didn't work because Explore makes an effort (in order to overcome its special evaluation rules, sigh...) to resolve assigned names.

The m in the denominator of 200*Unit(kg)/m^3 was not entered as units.

It's just the name m. The visual cue of that is that it was rendered in italics in your output.

Below, I entered in as the unit m, instead.

restartNULL

Pi

Pi

"(->)"

3.1415926535897932385

 ``

Dyson[Inside] := 355016117130*Unit('ft')NULL

NULLDyson[Outside] := 355016117280*Unit('ft')

355016117280*Units:-Unit(ft)

NULLNULL

Vol[Dyson] := (4/3)*Pi*Dyson[Outside]^3-(4/3)*Pi*Dyson[Inside]^3 = 75621866085188578208340000*Pi*Units:-Unit(ft)^3"(->)"0.23757309894397957472e27*Units:-Unit(ft)^3

   

Total_Weight := 200*Unit('kg')*Vol[Dyson]/Unit('m')^3

15124373217037715641668000000*Units:-Unit(kg)*Pi*Units:-Unit(ft)^3/Units:-Unit(m)^3

"(=)"

(6691789940422039289112993734304/15625)*Pi*Units:-Unit(kg)

"(->)"

0.13454641994365413638e28*Units:-Unit(kg)NULL

Download Q_for_Dyson_Shell_Weight_ac.mw

Are you trying for one of these?

nprintf("#mrow(%a,mo(\"&equals;\"),%a,mo(\"&equals;\"),%a)",
        Typesetting:-Typeset(P(X <= 5)),
        Typesetting:-Typeset(P(-5 <= -X)),
        Typesetting:-Typeset(P(E(X)-5 <= E(X)-X)));
 

`#mrow(Typesetting:-mrow(Typesetting:-mi("P"),Typesetting:-mo("&ApplyFunction;"),Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mi("X"),Typesetting:-mo("&le;"),Typesetting:-mn("5")))),mo("&equals;"),Typesetting:-mrow(Typesetting:-mi("P"),Typesetting:-mo("&ApplyFunction;"),Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mn("&uminus0;5"),Typesetting:-mo("&le;"),Typesetting:-mrow(Typesetting:-mo("&uminus0;"),Typesetting:-mi("X"))))),mo("&equals;"),Typesetting:-mrow(Typesetting:-mi("P"),Typesetting:-mo("&ApplyFunction;"),Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi("E"),Typesetting:-mo("&ApplyFunction;"),Typesetting:-mfenced(Typesetting:-mi("X"))),Typesetting:-mo("&minus;"),Typesetting:-mn("5")),Typesetting:-mo("&le;"),Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi("E"),Typesetting:-mo("&ApplyFunction;"),Typesetting:-mfenced(Typesetting:-mi("X"))),Typesetting:-mo("&minus;"),Typesetting:-mi("X"))))))`

`And`(P(X <= 5) = P(-5 <= -X),
      P(-5 <= -X) = P(E(X)-5 <= E(X)-X));
 

And(P(X <= 5) = P(-5 <= -X), P(-5 <= -X) = P(E(X)-5 <= E(X)-X))

`%and`(P(X <= 5) = P(-5 <= -X),
       P(-5 <= -X) = P(E(X)-5 <= E(X)-X));
  

%and(P(X <= 5) = P(-5 <= -X), P(-5 <= -X) = P(E(X)-5 <= E(X)-X))

Download ts_inert_ex.mw

Note that the latter two are inert, ie. not merely unevaluated by quotes (which might be subsequently ruined by inadvertant evaluation if passed around as arguments, etc).

Is this the kind of thing you're after?

It's not entirely clear what precisely you mean by, "...if R is not zero". Do you mean just when R is not a 4x4 Matrix with all zero as entries, or do you mean whenever R is some other 4x4 Matrix, or when R is any other size Matrix, or when R<>0 (scalar), or something else...?

restart;

 

p := proc(R)
  if R::'Matrix(4,4)' and LinearAlgebra:-Norm(R)=0 then
      A;
  else
      B;
  end if;
end proc:

 

p( Matrix(4,4) );

A

p( Matrix(4,4,12.3) );

B

# adjust code if you'd rather this threw error, or other
p( Matrix(2,2) );

B

Download Matrix_type_chk.mw

First 32 33 34 35 36 37 38 Last Page 34 of 336