acer

17617 Reputation

29 Badges

14 years, 190 days

On google groups. (comp.soft-sys.math.maple and sci.math.symbolic)

On stackoverflow.

On math.stackexchange.com.

MaplePrimes Activity


These are answers submitted by acer

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

Here are some ideas:

restart;

T:=325:
a12:=2305.28444347652 - 9.14490843016421*T + 0.00680052257590234*T^2:
a21:=-6665.24838284836 + 46.0897018087247*T - 0.0694991633494123*T^2:
alfa:=0.3:
x2:=1-x1:
tau12:=a12/T:
tau21:=a21/T:
G12:=exp(-alfa*tau12):
G21:=exp(-alfa*tau21):
lng1:=x2^2*(tau21*(G21/(x1+x2*G21))^2+tau12*(G12/((x2+x1*G12)^2))):
lng2:=x1^2*(tau12*(G12/(x2+x1*G12))^2+tau21*(G21/((x1+x2*G21)^2))):
lnga1:=subs(x1=xa1,lng1):
lngb1:=subs(x1=xb1,lng1):
lnga2:=subs(x1=xa1,lng2):
lngb2:=subs(x1=xb1,lng2):
r1:=lnga1+ln(xa1)=lngb1+ln(xb1):

r2:=lnga2+ln(1-xa1)=lngb2+ln(1-xb1):

ans1 := fsolve(abs~((lhs-rhs)~([r1,r2])));
 

{xa1 = -7.290257776, xb1 = -7.290257776}

map(lhs-rhs, eval([r1,r2], ans1));

[0.+0.*I, 0.]

ans2 := fsolve(abs~((lhs-rhs)~([r1,r2])), {xa1,xb1}, avoid={ans1});

{xa1 = -7.290476897, xb1 = -7.290476897}

map(lhs-rhs, eval([r1,r2], ans2));

[0.+0.*I, 0.]

ans3 := fsolve(abs~((lhs-rhs)~([r1,r2])), {xa1,xb1}, avoid={ans1, ans2});

{xa1 = -1.173033030, xb1 = -1.173033030}

map(lhs-rhs, eval([r1,r2], ans3));

[0.+0.*I, 0.]


 

ansC := fsolve([r1,r2],{xa1=-10-I..10+I, xb1=-10-I..10+I});

{xa1 = 1.136495104-.9634251620*I, xb1 = 1.136495104-.9634251620*I}

map(lhs-rhs, eval([r1,r2], ansC));

[0.+0.*I, 0.+0.*I]

 

Download NRTL_ac.mw

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