acer

17602 Reputation

29 Badges

14 years, 188 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

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

There is no need to use the implicitplot command for this, and indeed it makes little sense to do so. The result using implicitplot computes less efficiently, is less smooth without resorting to extra options (eg. gridrefine), and is generally problematic for export as regular numeric data.

The above is also true if you made a simple sign mistake in the left-hand side of your equation.

It's easy enough to flip the axes (without requiring any additional plottools command) by simply using the parametric calling sequence of the plot command.

It's easy enough to apply either Re or Im to the parts of the explicit solution, if solving B explicitly in terms of a.

Or you could solve a explicitly in terms of B and plot directly, and omit any use of the Re command. (As suggested, if there is a sign mistake then a plot similar to described arises.)

Ensure that the expression is as you want it, ie. no mistakes. Then figure out which portion of these you want to plot.

Look at the explicit formula which you choose to plot (whether it be B in terms of a, or a in terms of B) and check whether it makes sense to you.

You can exclude any of the options like view, colorlabels, gridlines, axis[2], tickmarks, etc, which are there mostly to modify the aesthetic. If in doubt, call the plot command without them, and add as desired.
 

restart;

Digits:=15:

eqn := (7.72-7.72*B)*(-7.717267500*a) = 662204.4444*B^2;

-7.717267500*(7.72-7.72*B)*a = 662204.4444*B^2

A := solve(eqn,a);

11115.0452892842*B^2/(-1.+B)

plot([A,B,B=0..1], labels=[a,B],
     view=[default,default], gridlines,
     axis[2]=[tickmarks=[seq(0..1.0,0.1)],
              location=low]
    );

S:=[solve(eqn,B)];

[0.449840722180451e-4*a+0.226516148099715e-12*(0.394383920330943e17*a^2-0.175343805433755e22*a)^(1/2), 0.449840722180451e-4*a-0.226516148099715e-12*(0.394383920330943e17*a^2-0.175343805433755e22*a)^(1/2)]

plot(map(Re,S),a=-80000..80000, labels=[a,Re(B)], color="Burgundy");

plot(S[1],a=-400000..0, labels=[a,B],
     view=[-400000..0,0..1], gridlines,
     axis[2]=[tickmarks=[seq(0..1.0,0.1)],
              location=low]
    );

plot(Re(S[1]),a=0..20000, labels=[a,B],
     view=[default,0..1], gridlines,
     axis[2]=[tickmarks=[seq(0..1.0,0.1)],
              location=low]
    );

 

Download flipped.mw

Do you mean the Help page with Topic  rtable_indexing ?

Or are you looking for help with function calls of an indexed name, such as, say,

F[bar]( some_args )

Or are you looking for the (uncommon) prefix-form indexing operator,  `?[]` , eg,

`?[]`(K,[1,2,3]);

               K[1, 2, 3]

You can use solve(E) to obtain all the roots of this particular example exactly.

You could also use fsolve(E,x,complex) if floating-point approximations were acceptable to you.

The select and selectremove commands can be used to separate the list of roots for various properties.

restart;

E:=2*x^(7)-x^(6)-4*x^(5)-15*x^(4)-14*x^(3)+19*x^(2)+28*x+30;

2*x^7-x^6-4*x^5-15*x^4-14*x^3+19*x^2+28*x+30

All:=[solve(E)];

[5/2, -1/2-((1/2)*I)*11^(1/2), -1/2+((1/2)*I)*11^(1/2), 2^(1/2), -2^(1/2), -1/2-((1/2)*I)*3^(1/2), -1/2+((1/2)*I)*3^(1/2)]

R,C := selectremove(is,All,real);

[5/2, 2^(1/2), -2^(1/2)], [-1/2-((1/2)*I)*11^(1/2), -1/2+((1/2)*I)*11^(1/2), -1/2-((1/2)*I)*3^(1/2), -1/2+((1/2)*I)*3^(1/2)]

CP,CN := selectremove(c->is(Im(c)>0), C);

[-1/2+((1/2)*I)*11^(1/2), -1/2+((1/2)*I)*3^(1/2)], [-1/2-((1/2)*I)*11^(1/2), -1/2-((1/2)*I)*3^(1/2)]

# real roots
R;

[5/2, 2^(1/2), -2^(1/2)]

# complex roots with strictly positive imaginary component
CP;

[-1/2+((1/2)*I)*11^(1/2), -1/2+((1/2)*I)*3^(1/2)]

 

Download select_example.mw

It is just as easy to separate out the roots which are real and positive, etc.

restart;

h1 := -I*(c1*exp(A)-c2*exp(-A)):
h2 := c2*exp(A)+c1*exp(-A):
A := (I*a/2)*(x+b*t):

extra := [c1=sqrt(a+2*lambda)/a,
          c2=sqrt(a-2*lambda)/a]:

h1*h2;

-I*(c1*exp(((1/2)*I)*a*(b*t+x))-c2*exp(-((1/2)*I)*a*(b*t+x)))*(c2*exp(((1/2)*I)*a*(b*t+x))+c1*exp(-((1/2)*I)*a*(b*t+x)))

ans := simplify(eval(combine(convert(h1*h2,trig)),extra));

(2*(a+2*lambda)^(1/2)*(a-2*lambda)^(1/2)*sin(a*(b*t+x))-(4*I)*lambda)/a^2

simplify( eval(h1*h2,extra) - ans );

0

 

Download simptrig_example.mw

You could get that even smaller under certain assumptions on a and lambda. For example,

simplify(combine(ans)) assuming a>2*lambda;

(2*(a^2-4*lambda^2)^(1/2)*sin(a*(b*t+x))-(4*I)*lambda)/a^2

 

There were not many edits or syntax mistakes, to my eye.

I had to change two instances of square brackets to round brackets, for delimiting a term of an arithmetic expression.

I had to change the procedure parameters from r[D] and t[D] to r__D and t__D (as typed in 2D Input mode, for literal subscripts which are symbols in the proper sense).

Your use of exp(...) in 2D Input mode was fine. (Evidence for that, in advance, was the fact that the "e" rendered in Roman upright rather than italic.)

I turned the a(i) computations into a loop, to shorten it.

Please check it over.

restart

 

Find 1st 40 roots of this equation.

 

BesselY(1, alpha[n])*BesselJ(1, alpha[n]*r[eD])-BesselJ(1, alpha[n])*BesselY(1, alpha[n]*r[eD]) = 0

 

r[eD] := 4.5

4.5

"rootfunction( beta ) :=(Y)[1](beta)* (J)[1](beta *r[eD]) - (J)[1](beta) * (Y)[1](beta * r[eD]) "

proc (beta) options operator, arrow, function_assign; BesselY(1, beta)*BesselJ(1, beta*r[eD])-BesselJ(1, beta)*BesselY(1, beta*r[eD]) end proc

Define vector to store roots.

 

a := Vector(40, fill = 0)``

NULL

Initialize the root finding package.

==========================

``

with(RootFinding)

``

Find roots for the given r[eD].

====================

 

a(1) := NextZero(rootfunction, 0.1e-3)

"for i from 2 to 40 do    a(i):=NextZero(rootfunction,a(i-1)); end do:"
   

a

Vector[column](%id = 18446884500815317822)

NULL

a(1)

.9609795905

a(40)

35.90623590

 

Define this function:

p__D := proc (r__D, t__D) options operator, arrow; ((1/2)*r__D^2+2*t__D)/(r[eD]^2-1)-r[eD]^2*ln(r__D)/(r[eD]^2-1)-(1/4)*(3*r[eD]^4-4*r[eD]^4*ln(r[eD])-2*r[eD]^2-1)/(r[eD]^2-1)^2+Pi*add(BesselJ(1, a(i)*r[eD])^2*(BesselJ(1, a(i))*BesselY(0, a(i)*r__D)-BesselY(1, a(i))*BesselJ(0, a(i)*r__D))*exp(-a(i)^2*t__D)/(a(i)*(BesselJ(1, a(i)*r[eD])^2-BesselJ(1, a(i))^2)), i = 1 .. 40) end proc

p__D(1, 1)

.8021455142

p__D(.5, .5)

1.384740257

``

Download DifficultFunctionac.mw

You could try to install the last point-release update to Maple 2019.2.1 .

(I am not sure that it will help your case, but it's possible that point-release fixed a rare issue with licensing misidentification.)

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