acer

32622 Reputation

29 Badges

20 years, 44 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

There are a few things that can be said.

1) No, it is not necessary to use active int for this. In fact doing so can be unnecessarily very expensive, since it can cause the integral to be attempted symbolically for each new call to the objective procedure (where each might fail, after a while). See first example in the attached.

There can sometimes be a mixup for Mimimize when called with procedure (rather than expression) form for the objective. It happens more rarely then it did long ago. It's caused mostly by some muddle about generating a procedure for evaluating the derivative numerically, via automatic differentiation. My first example does a few things to avoid that.

2) Your second bullet point sounds like it might be an evaluation mix-up. See the last example in the attachment with proc4, which may be what you were looking for here. The point is to keep straight that the K,r,x in globally defined df are themselves global (whereas the K,r in proc4 are procedure parameters). It helps to always test the objective procedure at a single candidate point, as sanity-check. Sorry if I misunderstand you here.

3) There's a battle between the numerical optimizer's wish to establish an optimum that satisfies its optimalitytolerance control option's value versus the numeric integrator's specification to only provide a result that needs to be convergent up to its own accuracy tolerance option epsilon. I know that sounds like a mouthful, but analogous struggles can happen between evalf/Int, fsolve, and Minimize. The "outer" one needs the "inner" one to be accurate enough and stable enough that it can establish its own criteria for convergence. So, in your later example I have forced values for those options. Sometimes it requires fiddling, also because evalf/Int might also require high enough working precision to attain its target accuracy.

Unfortunately, I lost patience fiddling and threw in a scaling factor, to get the integrand to a nice range. It works, but muddies the question of what is really needed.

I have used unapply within the call to evalf(Int(...)) since I've found that all too often it is a useful way to prevent the numeric integrator's zealous approach to checking for discontinuities when presented an integrand in expression form. I didn't check whether it was completly necessary here.

4) I didn't accomodate the case that the numeric integration ever fails to produce a numeric value. The proc3 could (and probably should) test that the numeric integration has returned a numeric value. If it hasn't then it your luck depends on which minimizer you're using, and how it handles ever getting a nonnumeric value for the objective. Some won't handle the discontinuity that a default "large" fallback value can cause. Some won't handle a nonnumeric Float(undefined). DirectSearch handles occasional nonnumeric values from the objective better. 

5) Keep in mind that for the multivariate case Optimization is a local rather than a global optimization package. (The add-on package DirectSearch provides global optimizers.)

The attached worksheet runs pretty quickly. But I had to fiddle to get that.

Practice problem using Optimization:-Minimization

restart

with(Optimization)

proc1 := proc (K, r) local f1, f2, x; if not [K, r]::(list(numeric)) then return ('procname')(args) end if; f1 := exp(.3*x); f2 := K/((K-1)*exp(-r*x)+1.0); evalf(Int((f1-f2)^2, x = 0 .. 3.2)) end proc

Minimize(proc1, 10 .. 100, 0 .. 1)

[0.243305733083052001e-4, Vector[column](%id = 18446883878292498062)]


Actual problem:

f1 := exp(-2.52428056431082+.468950959387562*x-0.403934021717953e-2*x^2); f2 := K*P0/(exp(-r*x)*K-exp(-r*x)*P0+P0); P0 := eval(f1, x = 0)

0.8011593034e-1

An approximation by hand:

plot([f1, eval(f2, [K = 0.25e5, r = .34])], x = 0 .. 38, 0 .. 0.15e5, legend = ["f1", "f2"])

``

proc3 := proc (K, r) local f1, f2, x; global P0; Digits := 15; if not [K, r]::(list(numeric)) then return ('procname')(args) end if; P0 := 0.8011593034e-1; f1 := exp(-2.52428056431082+.468950959387562*x+(-1)*0.403934021717953e-2*x^2); f2 := K*P0/(exp(-r*x)*K-exp(-r*x)*P0+P0); 0.1e-6*evalf(Int(unapply((f1-f2)^2, x), 0 .. 38, epsilon = 0.1e-5, method = _Dexp)) end proc

sol := Minimize(proc3(K, r), K = 0.10e5 .. 0.60e5, r = .2 .. .4, optimalitytolerance = 0.1e-8)

[.184857934786779005, [K = HFloat(17139.92551737672), r = HFloat(0.34703798628407145)]]

plots:-display(
  plot3d(proc3(K, r), K = 1000. .. 60000., r = 0.2 .. 0.4),
  plots:-pointplot3d([eval([K,r,sol[1]],sol[2])],
                     symbolsize=20, symbol=solidsphere, color=red),
  view=0.0 .. 1e1,
  size=[400,400]);

df := (f1-f2)^2

(exp(-2.52428056431082+.468950959387562*x-0.403934021717953e-2*x^2)-0.8011593034e-1*K/(exp(-r*x)*K-0.8011593034e-1*exp(-r*x)+0.8011593034e-1))^2

proc4 := proc (K, r) Digits := 15; if not [K, r]::(list(numeric)) then return ('procname')(args) end if; 0.1e-6*evalf(Int(unapply(eval(df, [:-K = K, :-r = r]), x), 0 .. 38, epsilon = 0.1e-5, method = _Dexp)) end proc

 

proc3(1010, .35), proc4(1010, .35)

35.0925066706428, 35.0925066706428

Minimize(proc4(K, r), K = 0.10e5 .. 0.60e5, r = .2 .. .4, optimalitytolerance = 0.1e-8)

[.184857934786779005, [K = HFloat(17139.92551737672), r = HFloat(0.34703798628407145)]]

 

Download MaplePrimes_Minimization_ac.mw

restart;

expr := y[1](t)*y[2](t)*(y[1](t)+y[2](t))^3:

frontend(diff, [expr, y[1](t)]);

                                  3                                        2
       y[2](t) (y[1](t) + y[2](t))  + 3 y[1](t) y[2](t) (y[1](t) + y[2](t))

Physics:-diff(expr, y[1](t));

                                  3                                        2
       y[2](t) (y[1](t) + y[2](t))  + 3 y[1](t) y[2](t) (y[1](t) + y[2](t))

subs(y1t=y[1](t),diff(subs(y[1](t)=y1t,expr),y1t));

                                  3                                        2
       y[2](t) (y[1](t) + y[2](t))  + 3 y[1](t) y[2](t) (y[1](t) + y[2](t))

thaw(diff(subs(y[1](t)=freeze(y[1](t)),expr),freeze(y[1](t))));

                                  3                                        2
       y[2](t) (y[1](t) + y[2](t))  + 3 y[1](t) y[2](t) (y[1](t) + y[2](t))

By "group" do you mean sort, or separate?

If you mean sort then you could use the sort command. If you mean separate then you could use the select (or selectremove command).

In any case, it would help if you explained how you wish to compare (for sorting) or discriminate (for selecting) between the eigenvalues.

Are they all purely real?

Here is an example of sorting a list of lists by the real value in each inner list's second entry.

L := [ ["SA",3], ["HF",2], ["JU",7], ["AK",1] ]:

sort( L, (a,b) -> a[2] < b[2] );

    [["AK", 1], ["HF", 2], ["SA", 3], ["JU", 7]]

sort( L, (a,b) -> a[2] > b[2] );

    [["JU", 7], ["SA", 3], ["HF", 2], ["AK", 1]]

If the second entries were complex-valued then you could also sort them by modulus or argument (or some combination, using and, etc).

The sort command has some powerful functionality, and I've only given a simple example. There are also matters like efficiency (eg. see the key option), ease-of-use, etc, which may or may not be important to your case.

sort( L, key = (a -> a[2]) );

    [["AK", 1], ["HF", 2], ["SA", 3], ["JU", 7]]

You might consider looking at the Help page for the sort command.

Here is a very simple example of separation according to a test on the second entry of each inner list. Again, you could make the selection procedure more complicated.

(L1,L2) := selectremove(a->type(a[2],odd), L):

L1;

      [["SA", 3], ["JU", 7], ["AK", 1]]

L2;

                [["HF", 2]]

If by "group" you meant something more complicated then you should tell us what you mean in detail. It might be that you're after functionality attainable using a command from the ListTools package (eg, its Classify, or Group commands, etc).

Check your syntax. It looks like some of your bracketed terms might accidentally be entered as part of function calls rather than as parts of any multiplication.

For example, it might be that you have the function call,
    UA_C(T_c-T_C)
while what you intended was the multiplication,
    UA_C*(T_c-T_C)

Otherwise, you won't be able to eliminate T_c, because it's part of an abstract function call.

This syntax issue might be a problem in more than one place (eg, also on UA_H(...) at least...).

We can't tell for sure because you have only provided us with an image of your code rather than actual code. You can upload and attach a Document using the green up-arrow in the Mapleprimes editor's menubar.

Your condition that all the variables are positive integers violates the constraint. It is violated even if the values are all as small as possible, ie. they have a value of 1.

restart

with(Optimization):

Maximize(5*x1+2*x2+7*x3+6*x4+x5+6*x6+8*x7+6*x8, {x1 >= 1, x2 >= 1, x3 >= 1, x4 >= 1, x5 >= 1, x6 >= 1, x7 >= 1, x8 >= 1, 10*x1+9*x2+15*x3+3*x4+11*x5+6*x6+3*x7+4*x8 <= 22}, assume = nonnegint)

Error, (in Optimization:-LPSolve) no feasible point found for LP subproblem

eval(10*x1+9*x2+15*x3+3*x4+11*x5+6*x6+3*x7+4*x8, [x1 = 1, x2 = 1, x3 = 1, x4 = 1, x5 = 1, x6 = 1, x7 = 1, x8 = 1]);

61

``

Positive_integer_ac.mw

Perhaps you meant this:

Maximize(5*x1+2*x2+7*x3+6*x4+x5+6*x6+8*x7+6*x8,
         {10*x1+9*x2+15*x3+3*x4+11*x5+6*x6+3*x7+4*x8<=22},
         assume = nonnegint);

    [56, [x1 = 0, x2 = 0, x3 = 0, x4 = 0, x5 = 0,
          x6 = 0, x7 = 7, x8 = 0] ]

If you state which solution you want (and what ranges for x and m) then you could adapt the rootfinding call appropriately.

The `plot` below could also be made more clear (filling in the gap). But you could also figure out how much accuracy you want.

The code has two ways to do the rootfinding (Calculus1:-Roots which utilizes fsolve's avoid option, and fsolve with its maxsols option which utilizes RootFinding:-NextZero).

restart;

N := q -> (exp(q)/sqrt(q))*sqrt(3/2)*Dw(sqrt(3*q/2));

proc (q) options operator, arrow; exp(q)*sqrt(3/2)*Dw(sqrt((3/2)*q))/sqrt(q) end proc

Dw := x -> exp(-x^2)*Int(exp(t^2), t=0..x);

proc (x) options operator, arrow; exp(-x^2)*(Int(exp(t^2), t = 0 .. x)) end proc

this := Dw(sqrt(3*q/2));

exp(-(3/2)*q)*(Int(exp(t^2), t = 0 .. (1/2)*6^(1/2)*q^(1/2)))

eq3a := q=(24/(n-3))*(x*k*m)/t;
eq3b := eval(eq3a, [t=6/10,n=6,k=13/10]);

q = 24*x*k*m/((n-3)*t)

q = (52/3)*x*m

eq4a := m = 1/2*((exp(q)/(q*N(q)))-1-1/q);

m = (1/6)*6^(1/2)/(q^(1/2)*exp(-(3/2)*q)*(Int(exp(t^2), t = 0 .. (1/2)*6^(1/2)*q^(1/2))))-1/2-(1/2)/q

eq4b := combine(value(eval(eq4a, eq3b))) assuming real;

m = (1/26)*26^(1/2)*exp(26*x*m)/(erfi((x*m)^(1/2)*26^(1/2))*(x*m*Pi)^(1/2))-(3/104)/(x*m)-1/2

M := proc(X)
  local res;
  option remember;
  res:=Student:-Calculus1:-Roots(eval(eq4b,x=X),m=-0.6 .. 0.6,numeric);
  #res:=[fsolve(eval(eq4b,x=X),m=-0.6 .. 0.6,maxsols=4)];
  res := [res[], undefined$(4-nops(res))];
end proc:
M1 := proc(X)
  if not X::numeric then return 'procname'(X); end if;
  M(X)[1];
end proc:
M2 := proc(X)
  if not X::numeric then return 'procname'(X); end if;
  M(X)[2];
end proc:
M3 := proc(X)
  if not X::numeric then return 'procname'(X); end if;
  M(X)[3];
end proc:
M4 := proc(X)
  if not X::numeric then return 'procname'(X); end if;
  M(X)[4];
end proc:

M(1.1);
M1(1.1), M2(1.1), M3(1.1), M4(1.1);

[-.4404638037, -0.4037552807e-1, .1382804516, .3786155189]

-.4404638037, -0.4037552807e-1, .1382804516, .3786155189

plot([M1(x), M2(x), M3(x), M4(x)], x = -3 .. 3, adaptive=false, numpoints=60);

plots:-implicitplot(rhs(eq4b)-lhs(eq4b), x=-3 .. 3, m=-0.5 .. 0.5, gridrefine=3);

 

Download dwson_m_x.mw

You can also do a similar thing for finding x as function of m. (It's actually easier and slightly quicker, avoids the gap, and the plot's axes could even be transposed).  dwson.mw

You could use an Empirical distribution or choose for this from RandomTools (or ProbabilityTable from the Statistics package).

I don't know whether you want individual (scalar) results, or lists, or Vectors. I don't know whether you want integers or floats. I don't know whether you'd prefer using Statistics or RandomTools, and I don't know whether you want easy ways to get just a few values or efficient ways to get many values.

So here are various ways.

restart;

with(RandomTools):

G := distribution(Empirical(<1,2,3,4,5>,
                            probabilities=[0.2,0.5,0.1,0.1,0.1])):

Generate(G);

HFloat(4.0)

seq(trunc(Generate(G)), i=1..10);

5, 1, 5, 2, 1, 2, 2, 5, 5, 1

L1 := Generate(list(choose([1$2,2$5,3,4,5]),10));

[2, 3, 2, 4, 5, 1, 3, 5, 4, 1]

# This is not fast
N := 10^5:
L2 := Generate(list(choose([1$2,2$5,3,4,5]),N)):

map(s->lhs(s)=evalf[2](rhs(s)/N), Statistics:-Tally(L2));

[1 = .20, 2 = .50, 3 = 0.99e-1, 4 = .10, 5 = .10]

with(Statistics):

H := RandomVariable(Empirical(<1,2,3,4,5>,
                              probabilities=[0.2,0.5,0.1,0.1,0.1])):

Sample(H, 1);

Vector(1, {(1) = 2.0})

trunc(Sample(H, 1)[1]);

3

V1 := trunc~(Vector(Sample(H, 10),datatype=anything));

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

L3 := trunc~(convert(Sample(H, 10),list));

[4, 4, 1, 2, 2, 2, 5, 3, 2, 1]

# This is not terribly slow
N := 10^5:
V2 := trunc~(Vector(Sample(H, N),datatype=anything)):

map(s->lhs(s)=evalf[2](rhs(s)/N), Statistics:-Tally(V2));

[1 = .20, 2 = .50, 3 = .10, 4 = .10, 5 = .10]

CodeTools:-Usage( Generate(list(choose([1$2,2$5,3,4,5]), 2*10^5)) ):

memory used=331.13MiB, alloc change=1.53MiB, cpu time=2.88s, real time=2.67s, gc time=439.41ms

CodeTools:-Usage( trunc~(Vector(Sample(H, 2*10^5),datatype=anything)) ):

memory used=19.85MiB, alloc change=6.11MiB, cpu time=68.00ms, real time=68.00ms, gc time=0ns

 

Download empirical_ac1.mw

You might also consider throwing in additional arrows normal to the line whose slope is the linear objective c1*x1+2*x2 .

One nice thing to illustrate is that the maximum is attained by moving [x1,x2] as far as possible (within the feasible region) in the direction of the normal to the objective line.

Below, I translate the line (so that it passes through the optimial point, give the current c1 value) since that gives a visual cue.

There are all kinds of variations on doing this.

If you really want to get fancy you could animate by the rotation angle (and compute the effective c1 value on the fly from that), which would produce a more even visual experience.

Another finesse might be to have a chunk of frames for the special case that the line is parallel to a side of the feasible region. In that case the whole edge of the region is optimal, and could be briefly all colored in green. Having a chunk of constant frames would provide a visual pause.

restart;

cnsts := [ 4*x1 + x2 <= 12,
             x1 - x2 >= 2,
                  x1 >= 0,
                  x2 >= 0 ]:

feasibleRegion := plots:-inequal(cnsts, x1 = 1 .. 4, x2 = -0.5 .. 1.5, nolines):

L := c1*x1+2*x2;

c1*x1+2*x2

m := solve(L,x2)/x1;

-(1/2)*c1

b := x2 - m*x1;

(1/2)*c1*x1+x2

Lform := Typesetting:-Typeset(Typesetting:-EV(L)):

 

All := proc(C1)
  local Nor,opt,P1,P2;
  uses plots, Optimization;
  opt := Maximize(eval(L,c1=C1), cnsts, x1=1..4, x2=-0.5..1.5);
  P1 := plot(eval(m*x1 + eval(b,opt[2]), c1=C1),
             x1=1..4, x2=-0.5..1.5, color=red,
             title=typeset('max'(L)=evalf[5](opt[1]),"\n",
                           "occurs at ",evalf[5](fnormal(opt[2])),"\n",
                           "when ", 'c1'=evalf[5](C1)));
  Nor := arrow(eval([x1,x2],opt[2]),
               <[C1,2]>, length=0.3, head_width=0.075, color=blue);
  P2 := pointplot([rhs~(opt[2])], color=green,
                  symbolsize=15, symbol=solidcircle);
  display(P1, P2, Nor, feasibleRegion,
          scaling=constrained);
end proc:

All(-0.4);

plots:-animate(All, [c1], c1=-20..20, paraminfo=false);

 

 

Download LP_anim2.mw

The basics of showing the feasible region, objective line, and optimal point are not complicated. There are several straightforward approaches to that, depending on how clear you want it to be for a reader who does not already know the underlying math.

Naturally, the code gets longer as you add in normal lines, a complicated title, etc.

Eventually you might discover that it's easier to learn and use a little Maple language syntax than to try and use the context-menus (Clickable Math, right-panel) for everything.

But here are some ways, including a little more use of the context-menus. (There are, of course, many other ways, including doing it fully by language commands.)

restart

with(LinearAlgebra)

tau := Matrix(3, 3, {(1, 1) = 200, (1, 2) = 0, (1, 3) = 0, (2, 1) = -50, (2, 2) = -800, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 200})

 

tau-lambda*(Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = 0, (3, 2) = 1, (3, 3) = 1})) = Matrix(%id = 18446883947032617798)"(->)"-(-200+lambda)^2*(800+lambda)"(->)"-(-200+lambda)^2*(800+lambda) = 0"(->)"[[lambda = -800], [lambda = 200], [lambda = 200]]"(->)"results``NULL

sigma := map[2](eval, lambda, results)

[-800, 200, 200]

sigma[1]

-800

sigma[2]

200

sigma[3]

200

eval(lambda, results[1])

-800

tau-lambda*(Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = 0, (3, 2) = 1, (3, 3) = 1})) = Matrix(%id = 18446883947011398582)"(->)"-(-200+lambda)^2*(800+lambda)"(->)"-(-200+lambda)^2*(800+lambda) = 0"(->)"eqn

R := [solve(eqn, lambda)]

[-800, 200, 200]

R[1]

-800

``

Download maple_ac.mw

You could try changing the syntax to, say,

   evalf(Int((r,a)->P(r)*a, [0..1, 0..1])):

or,

   evalf(Int('P'(r)*a, [r=0..1, a=0..1]));

Do you see that when you changed from single-argument to two-argument you also changed the calling sequence? You changed it to an invalid mashup of operator and expression form.

Alternatively, you could use P(r)*a as the expression, but you'd first have to modify the defn of P so that it returned unevaluated if passed symbol r rather than a number. One of your attempts this way encountered what's known as premature evaluation, where P received nonnumeric r and was not set up to deal with that. For example,

P := proc(r)
   if not r::numeric then return 'procname'(r); end if;
   ...proceed as usual...   
end proc

I prefer this last approach.

But you cannot sensibly use P*a as the integrand, where P is a procedure and a is an expression in one of more of the variables of integration.

If you upload your actual data file as an attachment in a Comment on the Question then I can easily show you all three methods at work.

As a side note, procedure P could be far more efficient if you wrote it to not have to re-form the data every time it gets called. I'd prefer to have the data, to check, before re-writing it all. You could substitute the actual data in lieu of a dummy name in the proc, or you could do it like, say,
   unapply('CurveFitting:-ArrayInterpolation'(...),numeric)
where use is again made of uneval quotes.

It's difficult to be precise when you haven't shown us the details of the actual data.

But I'm curious, did you want to escape a backslash in "VAGRANTASP\Adam" ?

Your final output seems like it may contain "VAGRANTASP\\Adam" with an escaped backslash. But your 2D Input contains a piecewise that compares against "VAGRANTASP\Adam" instead. Is that mismatch an oversight?

The expression is zero, for y=0 and 3 <= x <= 4 . So its maximum is not negative.

The RealDomain:-solve and solve commands are not super strong when tackling trig equations over specified ranges.

You could try this,

Student:-Calculus1:-Roots(sec(x) = sqrt(2), x=3*Pi*(1/2) .. 2*Pi);

            [7   ]
            [- Pi]
            [4   ]

I am not sure how it works on your Windows (command shell), but have you tried this kind of redirection?

ssystem("ffmpeg -i c:/test/test.mp3 2>&1");

What you were missing -- I suspect -- is that ffmpeg usually sends its media stream to stdout, and its diagnostics to stderr. I just checked ffmpeg's manual-page, and it states, "By default the program logs to stderr."

So, I had to redirect stderr in order to get ffmpeg's messages from Maple's ssystem. That was on Linux, and I expect it behaves similarly on Windows.

I'm also assuming that ssystem allows forward slash as directory separator on your Windows, and that part of your call was ok.

You could consider why this is ok here. And you might also apply sqrt afterwards.

minimize(1-(4*(x-1))*(1-y)/(x^2*(1+y)), x=3..4, y=0..1, location);

            1   /[                1]\ 
            -, { [{x = 3, y = 0}, -] }
            9   \[                9]/ 

I'll try and find time to look at what causes the original failure; it may be part of a pattern. (I was making a study of problematic cases a while back, but got busy.)

First 133 134 135 136 137 138 139 Last Page 135 of 339