acer

32338 Reputation

29 Badges

19 years, 326 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

I found it tough to get Maple to look at y=7, x=1..3 as one of the canditate ranges for which the objective is purely real. I was trying to code it in such a way as to give Maple as little help and prodding as I could.

One thing that made this tougher than it ought to be is that I wasn't sure how hard I should try and get Maple to go from something like, {x>=1,x<=3} to x=1..3.  (See the member of newS in the attached worksheet, which correspondes to this.) This aspect makes a difference because, Optimization:-Maximize seems to go awry when this restriction is passed as a pair of constraints, but succeeds when it is passed as simple bounds. See end of attached worksheet.

 

restart:

interface(warnlevel=0):

obj:=sqrt((x-1)*(y-x))+sqrt((7-y)*(1-x))+sqrt((x-y)*(y-7));

((x-1)*(y-x))^(1/2)+((7-y)*(1-x))^(1/2)+((x-y)*(y-7))^(1/2)

plot3d([Re,Im](obj), x=-2..3, y=0..11, color=[red,blue], axes=box);

S:=[solve(evalc(Im(obj)))]:  # takes a little while

newS:={seq(solve(s union {x>=-2,x<=3,y>=0,y<=11}),s in S)};

{{x = 1, y = 1}, {x = 1, y = 4}, {x = 1, y = 7}, {x = 1, 1 <= y, y <= 7}, {x = y, 0 <= y, y <= 1}, {y = 7, x <= 3, 1 < x}}

for s in newS do
   try
      seql,sineq:=selectremove(type,s,`=`);
      ans:=Optimization:-Maximize(eval(obj,seql), sineq,
                                  initialpoint=[x=1/2,y=1/2]);
      print(ans[1],[ans[2][],seql[]]);
   catch:
   end try:
end do;

2.64575131106459072, [y = HFloat(5.551115123125783e-17), x = y]

for s in newS do
   try
      seql,sineq:=selectremove(type,s,`=`);
      ans:=Optimization:-Maximize(eval(obj,seql), sineq,
                                       initialpoint=[x=1,y=5]);
      print(ans[1],[ans[2][],seql[]]);
   catch:
   end try:
end do;

3., [y = HFloat(4.000000000000039), x = 1]

-0., [y = HFloat(1.0), x = y]

#plots:-display(Array([
#   plot([Re,Im](eval(obj,x=1)),y=0..11,view=0..5),
#   plot([Re,Im](eval(obj,y=7)),x=-2..3,view=0..5),  # method above missed this one
#   plot([Re,Im](eval(obj,x=y)),y=0..1,view=0..5) ]));

 

Download toughopt1.mw

For some reason maplenet is not inlining the resst of the sheet. See the link for a bit more.

acer

Is it all T's except at the "square" positions? Let's consider H=0 and T=1, (This could be quite wrong, I've only had one coffee.)

> restart:   
                               
> seq((numtheory[tau](i)-1 mod 2),i=1..50); 

0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1,

    1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1

> evalf[30](sum((1/2^i)*'(numtheory[tau](i)-1 mod 2)',i=1..1000));

                        0.435531586394061420665270725728

> evalf[30](1-sum(1/2^(i^2),i=1..infinity));                      

                        0.435531586394061420665270725728

acer

So, your `eq` is just an equality constraint.

restart;

eq := x = 11/500+6*u*(1/25)-u^2
          +surd(-174200000000*x^5+13045131*E^4*u^4,5)/(65225655*E^4*u^4):

E:=0.1;
ans:=Optimization:-Maximize(x,{eq},u=0..1,x=0..1,initialpoint=[x=1,u=0.5]);
eval((rhs-lhs)(eq),ans[2]);
E:='E':

sols:=seq([E,Optimization:-Maximize(x,{eq},u=0..1,x=0..1,
                                    initialpoint=[x=1,u=0.5])],
          E=0.1..1,0.1);

seq( eval((rhs-lhs)(eq),[s[2,2][],E=s[1]]), s in sols );

acer

Try it with evalf(M) as the first argument to SingularValues.

(This is just an easily fixable bug.)

acer

This appears to be another instance of `sum` versus `add`.

acer

The solutions you've found are all "pretty good" fits to the data. It's possible to obtain a close fit with Optimization:-Minimzie rather than NonlinearFit, but the optimality tolerance has to be tightened.

Also, the solution space has more than just one local minimum, which makes it tough for a local optimizer such as Optimization:-Minimize which will return whichever local minimum it converges to (and not necessarily the globally optimal of all such local minima).

I've created a univariate objective, which fortunately happens to be pretty smooth for your example. But in general there's no reason for it to be so well behaved. (I feel that I got lucky, with the tricks.)

This is just the kind of example where I would expect a curvefitting routine based on global optimization techniques (such as DirectSearch offers) to do much better by default.

 

restart; with(Statistics)

Edata := [-.462124699683824428309502010103, -.527751016544377322226200091423, -.598698303952974890384265763998, -.674809251210729440119324699197, -.75601431564264356773357700158, -.84227676885295360042125209345, -.93357527229747985234963542418, -1.02989666231439422969552894410, -1.13123248073042925355643641564]:

Ez := proc (Z, Zc) options operator, arrow; -(1/2)*Zc^2-a*(Z-Zc)-b*(Z-Zc)^(3/2)-c*(Z-Zc)^2-d*(Z-Zc)^(5/2)+e*(Z-Zc)^3+f*(Z-Zc)^(7/2) end proc:

Ezpw := proc (Z, Zcc) options operator, arrow; piecewise(Zcc < Z, Ez(Z, Zcc)) end proc:

obj := add((eval(Ezpw(Z, Zc), Z = Zdata[i])-Edata[i])^2, i = 1 .. nops(Zdata)):

``

best := infinity:

Optimization:-Minimize(G, .9 .. .95);

[3.59737560152029*10^(-16), Vector(1, {(1) = .910849324507901})]

sol;

[0.359738773129374952e-15, [Zc = HFloat(0.9108493350527227), a = HFloat(1.1425316489552044), b = HFloat(0.1742451497775256), c = HFloat(0.7695947603372874), d = HFloat(0.14059245764891992), e = HFloat(0.023045433394426117), f = HFloat(0.008513627717854673)]]

plots:-display(Array([plot(G, .900 .. .95, view = 0 .. 0.1e-5, style = line), plot(G, .900 .. .93, view = 0 .. 0.1e-7, style = line), plot(G, .900 .. .92, view = 0 .. 0.1e-8, style = line)]))

 

Download maphelp_modif.mw

acer

I would disagree that the natural "factor" in your example is 1/3.

The current posting relates to a recent Question here. Alejandro mentioned that using ``() to block automatic simplification had issues because the reversing operation `expand` could hit too much. But this is also true of prefix `%` and `value`, since applying `value` can meddle with other inert subexpressions such as calls to `Int`, etc. The essential problem stems from trying to make a mechanism serve two purposes. It so often turns out in the long run to be a bad idea whenever anyone tries to hijack a pre-existing facility and try and make it do double-duty. So below I'll avoid both those mentioned ways, and try and use some made-up wrapping procedure (the concept of which has also been done a few times before on this site). As name for the wrapping and unwrapping functions one can try and choose names that aren't already used for other purposes.

restart:

`print/_`:=proc(t) ``(t); end proc:

block:=proc(expr) local c,r;
           c:=content(expr);
           r:=expr/c;
           _(c)*r;
       end proc:

unblock:=proc(expr)
             subsindets(expr,'specfunc(anything,_)',t->op(t));
         end proc:

a:= (1/3)*x^2 + (1/6)*x - (1/12);

                        1  2   1     1 
                        - x  + - x - --
                        3      6     12

A:=block(a);

                      /1 \ /   2          \
                      |--| \4 x  + 2 x - 1/
                      \12/   
              
unblock(A);

                        1  2   1     1 
                        - x  + - x - --
                        3      6     12

acer

Procedures are of type last_name_eval. (Visit that link and see near the end for a tip on returning a procedure as the result from another procedure.)

> Phi := unapply(A1[1], x1, x2);

                           Phi := (x1, x2) -> A1[1]

> Phi;

                                      Phi

> eval(Phi,1);

                               (x1, x2) -> A1[1]

> type(Phi,last_name_eval);

                                     true

acer

Try it with the exponent ``(1/2) instead of 0.5 where `` is two single left-quotes (name-quotes).

If it gets the behaviour you are after, then try applying the expand command to the entire expression and see if that gets you back to normal. (If it works then the sqrt(2) should pop out from the rest.)

acer

If your real, floating-point data is symmetric then try,

B := Matrix(..., shape=symmetric);

Eigenvalues(B);

so as to get an appropriate solver which returns a real result. The argument to `Matrix` above (where I had ... ellipsis) could be an Array, or listlist, or another Matrix the previous version of B, etc.

It's difficult to answer about that error message (indicating a missing argument) without seeing your code. Could you upload a self-contained worksheet that demonstated the problem?

acer

There certainly would be good uses for an "Interpolate" command which, given data, could return an interpolating procedure. Such a returned interpolating procedure might even be made clever in that it could know how to apply the integrated version of its own simple interpolatory formulae (splines, say). But such a thing does not yet exist, as described.

As you say, one need is for very low overhead in on-the-fly interpolation at a given (single) point, which might be done repeatedly but cannot all be done in bulk at once.

Having said that, it might still be possible to speed up your application a bit. If you don't have a great deal of data points then you might be able to use BSplineCurve and instead get an interpolating piecewise polynomial (which could be integrated, once). Or, if sticking with ArrayInterpolation then a little relief might be had with reusable 1-element float[8] Arrays for input/output. Or, possibly, for a fixed degree spline use in ArrayInterpolation then a particular differencing scheme might be able to quickly integrate "exactly".

Any chance you could upload a self-contained worksheet that implemented your current methodology?

acer

I think that your second example is misprogrammed. Since `result` is a table the procedure `make` should return `eval(result,1)` or `eval(result)`.

[edited] A Maple table has type `last_name_eval`. See the help-page last_name_eval which states, "A procedure that returns a table, module, or another procedure normally requires an evaluation at the return site, using the two-argument form of eval in which the second argument is 1."

acer

ee := cos(x^2+x)-sin(x+1/6*Pi)-cos(x+1/3*Pi);

                          2                 Pi              Pi
               ee := cos(x  + x) - sin(x + ----) - cos(x + ----)
                                            6               3

> ff : combine(expand(ee));                    

                                  2
                             cos(x  + x) - cos(x)

I'm not sure that we really need to rely on Maple to get from there to,

            x^2 + x = x + (2*n)*Pi

where n is any integer. Anyway, we can get some solutions, x = +-sqrt(2*n*Pi) for n=0,1,2...And there are other solutions which may be found numerically. There's no reason to stop this search at 2*Pi, as below. Note this this also produces the solutions described exactly above.

ee := cos(x^2+x)-sin(x+1/6*Pi)-cos(x+1/3*Pi);

Digits:=20: # your choice
FF:=unapply(combine(expand(ee)),x):
NZ:=RootFinding:-NextZero:         
S:={0.0}: last:=0.0:               
while is( last < 2*Pi ) do         
  last := NZ(FF,last);
  S := S union {NZ(FF,last)};
end do:
S;

  {0., 2.5066282746310005024, 2.6832554370229568635, 3.4552840449895851091,

    3.5449077018110320546, 4.1120192907224387652, 4.3416075273496059562,

    4.6934986199961384012, 5.0132565492620010048, 5.2208610210386085156,

    5.6049912163979286993, 5.7068843101888305450, 6.1399602476789309309,

    6.1599917917157343495, 6.5860838226726890216, 6.6319150439565422209}

So, is there are nice way to represent these,

evalf[10](S) minus evalf[10]({seq(evalf[20](sqrt(2*n*Pi)), n=0..7)});

   {2.683255437, 3.455284045, 4.112019291, 4.693498620, 5.220861021, 5.706884310,
    6.159991792, 6.586083823}

acer

All:=map(`*`,combinat[partition](15,5),10):

F:=(n::posint)->select(u->nops(u)=n,All):

F(3);

                         [[50, 50, 50]]

F(4);

    [[30, 40, 40, 40], [30, 30, 40, 50], [20, 40, 40, 50], 

      [20, 30, 50, 50], [10, 40, 50, 50]]

F(5);

         [[30, 30, 30, 30, 30], [20, 30, 30, 30, 40], 

           [20, 20, 30, 40, 40], [10, 30, 30, 40, 40], 

           [10, 20, 40, 40, 40], [20, 20, 30, 30, 50], 

           [10, 30, 30, 30, 50], [20, 20, 20, 40, 50], 

           [10, 20, 30, 40, 50], [10, 10, 40, 40, 50], 

           [10, 20, 20, 50, 50], [10, 10, 30, 50, 50]]

and so on up to F(15).

acer

I don't understand why the maplenet server used by Mapleprimes is causing the gridlines to show up below. The gridlines don't appear when run in 64bit Maple 16.01 on Windows 7.

 

restart:

with(Statistics):

(a,b):=-1,1:
X:=RandomVariable(Uniform(a,b)):

plot( PDF(X,x), x=5*a/2..5*b/2, view=0..3/2*1/(b-a), labels=[``,``],
      axes=box, discont=[symbolsize=1], thickness=2,
      tickmarks=[[a=a,b=b],[1/(b-a)=1/(b-a)]]);

DensityPlot( X, range=5*a/2..5*b/2, view=0..3/2*1/(b-a), axes=box,
             tickmarks=[[a=a,b=b],[1/(b-a)=1/(b-a)]]);

PDF(X,x); # why not piecewise(-1>x,0,1>=x,1/2,0) ??

piecewise(x < -1, 0, x < 1, 1/2, 0)

G:=piecewise(-1>x,0,1>=x,1/2,0);

G := piecewise(x < -1, 0, x <= 1, 1/2, 0)

plot( G, x=5*a/2..5*b/2, view=0..3/2*1/(b-a), labels=[``,``],
      axes=box, discont=[symbolsize=15,symbol=circle], thickness=2,
      tickmarks=[[a=a,b=b],[1/(b-a)=1/(b-a)]]);

a,b:=-1.0,1.0: # floats

plot( G, x=5*a/2..5*b/2, view=0..3/2*1/(b-a), labels=[``,``],
      axes=box, discont=[symbolsize=15,symbol=circle], thickness=2,
      tickmarks=[[a=a,b=b],[1/(b-a)=1/(b-a)]]);

 

 

Download uniform.mw

acer

First 255 256 257 258 259 260 261 Last Page 257 of 336