Robert Israel

6567 Reputation

21 Badges

18 years, 39 days
University of British Columbia
Associate Professor Emeritus
North York, Ontario, Canada

MaplePrimes Activity


These are replies submitted by Robert Israel

@Preben Alsholm :  Oops, yes.  That "output=list" was from an earlier version that somehow didn't get deleted.  Thanks for catching this.

@Markiyan Hirnyk : There's no need to switch to simplex: LPSolve gets it right.  But there's confusion about the second constraint.  You stated it as x2 + x5 = 2, but your matrices A and AA have corresponding row [1,0,0,0,1] instead of [0,1,0,0,1].  In any case: yes, the optimal solution of the modified problem is an alternative optimal solution of the original.

@Markiyan Hirnyk : There's no need to switch to simplex: LPSolve gets it right.  But there's confusion about the second constraint.  You stated it as x2 + x5 = 2, but your matrices A and AA have corresponding row [1,0,0,0,1] instead of [0,1,0,0,1].  In any case: yes, the optimal solution of the modified problem is an alternative optimal solution of the original.

@Markiyan Hirnyk: You are using inequality constraints rather than equality constraints.  Try:

LPSolve(cc,[NoUserValue,NoUserValue,AA,bb],maximize,assume=nonnegative);

@Markiyan Hirnyk: You are using inequality constraints rather than equality constraints.  Try:

LPSolve(cc,[NoUserValue,NoUserValue,AA,bb],maximize,assume=nonnegative);

If the dual problem has a nondegenerate optimal solution, the primal problem has a unique optimal solution.  If the dual problem has a degenerate optimal solution, the optimal solution of the primal may or may not be unique.  To see if the primal has a nonunique solution, you may try something like this.  Suppose the primal is

maximize cTx
subject to

A x = b

x >= 0

and you get an optimal solution with objective value v in which the variables xi for i in set N are 0.  Then consider the problem

maximize Sum(x[i], i = N)
subject to

A x = b

cTx = v


x >= 0

If the optimal solution of the original problem is unique, it will also be the optimal solution of this problem.  If the optimal solution of the original problem is nonunique, this one will either be unbounded or have another optimal solution, which will be a different solution of the original problem.

Caution: this is assuming exact arithmetic.  Roundoff errors in floating-point arithmetic complicate matters.  In particular, it can be difficult to know whether something is really 0 or just close to 0.

If the dual problem has a nondegenerate optimal solution, the primal problem has a unique optimal solution.  If the dual problem has a degenerate optimal solution, the optimal solution of the primal may or may not be unique.  To see if the primal has a nonunique solution, you may try something like this.  Suppose the primal is

maximize cTx
subject to

A x = b

x >= 0

and you get an optimal solution with objective value v in which the variables xi for i in set N are 0.  Then consider the problem

maximize Sum(x[i], i = N)
subject to

A x = b

cTx = v


x >= 0

If the optimal solution of the original problem is unique, it will also be the optimal solution of this problem.  If the optimal solution of the original problem is nonunique, this one will either be unbounded or have another optimal solution, which will be a different solution of the original problem.

Caution: this is assuming exact arithmetic.  Roundoff errors in floating-point arithmetic complicate matters.  In particular, it can be difficult to know whether something is really 0 or just close to 0.

Without TEXT, you'd have to work out a way to read a font file (or would you want to write your own and include it in the program?) and render the characters.  It sounds like a lot of work.  Since TEXT is available, I don't really see the point of that.

How about this?  Make sure you make your plotting window large enough to see the whole code.

> A:= [
"bjo9IG5vcHMoQSk6IHdpdGgocGxvdHMpOg==",
 "ZGlzcGxheShb",
 "dGV4dHBsb3QoWzAsMCwiQTo9IFsiXSxhbGlnbj1yaWdodCks",
 "c2VxKHRleHRwbG90KFswLC1pLA==",
 "Y2F0KCJcIiIsQVtpXSwiXCIsIildLA==",
 "YWxpZ249cmlnaHQpLCBpPTEuLm4tMSks",
 "dGV4dHBsb3QoWzAsLW4sIGNhdCgiXCIiLCBBW25dLA==",
 "ICJcIl07IildLGFsaWduPXJpZ2h0KSw=",
 "c2VxKHRleHRwbG90KFswLC1uLWks",
 "U3RyaW5nVG9vbHNbRGVjb2RlXShBW2ldLGJhc2U2NCldLA==",  
 "YWxpZ249cmlnaHQpLGk9MS4ubildLA==",
 "YXhlcz1ub25lLHZpZXc9WzAuLjEuNSxERUZBVUxUXSk7"]:
n:= nops(A): with(plots):
 display([
 textplot([0,0,"A:= ["],align=right),
 seq(textplot([0,-i,
 cat("\"",A[i],"\",")],
 align=right), i=1..n-1),
 textplot([0,-n, cat("\"", A[n],
 "\"];")],align=right),
 seq(textplot([0,-n-i,
 StringTools[Decode](A[i],base64)],
 align=right),i=1..n)],
axes=none, view=[0..1.5,DEFAULT]);

@stefano : Are you sure you have the right system?  Why do you write g(t) - 1.5*g(t) and h(t) - 1.5*h(t) instead of -g(t)/2 and -h(t)/2?

I tried finding series solutions.  You didn't provide initial conditions, so for convenience I assumed g(0) = h(0) = 0.  Then of course invg(0) = invh(0) = 0, and it should make sense to take all series around 0. Thus:

> eq1:= (3/2*f(t) - t)*(diff(g(t),t)*h(t) + diff(h(t),t)*g(t))-g(t)*h(t);
 eq2:= (3/2*g(t) - t)*(diff(f(t),t)*h(t) + diff(h(t),t)*f(t))-f(t)*h(t)
+ (g(t) - 3/2*g(t))*g(invh(g(t)))*f(invh(g(t)))*diff(f(t),t);
 eq3:= (3/2*h(t) - t)*(diff(g(t),t)*f(t) + diff(f(t),t)*g(t))-f(t)*g(t)
+ (h(t) - 3/2*h(t))*f(invg(h(t)))*(diff(g(t),t)*f(t)+diff(f(t),t)*g(t));
 eq4:= h(invh(t))-t;
 eq5:= g(invg(t))-t;
ord:= 8; # order for the series
eqs8:= series~(eval([eq1,eq2,eq3,eq4,eq5],
[f=unapply(add(F[i]*t^i,i=0..ord),t), g=unapply(add(G[i]*t^i,i=1..ord),t),
h=unapply(add(H[i]*t^i,i=1..ord),t), invg=unapply(add(GI[i]*t^i,i=1..ord),t),
invh=unapply(add(HI[i]*t^i,i=1..ord),t)]),t,ord);
solve(coeffs~(convert~(%,polynom),t));


Interestingly enough, although this does get me a solution to order 8, if I change ord to 9 there is no solution.  I found that rather surprising: usually if something like this works for low orders, it will also work to arbitrary orders.  But it appears this is a true result, not a flaw in the Maple code: there is no solution as formal power series in t with g(0) = h(0) = 0.

@stefano : Are you sure you have the right system?  Why do you write g(t) - 1.5*g(t) and h(t) - 1.5*h(t) instead of -g(t)/2 and -h(t)/2?

I tried finding series solutions.  You didn't provide initial conditions, so for convenience I assumed g(0) = h(0) = 0.  Then of course invg(0) = invh(0) = 0, and it should make sense to take all series around 0. Thus:

> eq1:= (3/2*f(t) - t)*(diff(g(t),t)*h(t) + diff(h(t),t)*g(t))-g(t)*h(t);
 eq2:= (3/2*g(t) - t)*(diff(f(t),t)*h(t) + diff(h(t),t)*f(t))-f(t)*h(t)
+ (g(t) - 3/2*g(t))*g(invh(g(t)))*f(invh(g(t)))*diff(f(t),t);
 eq3:= (3/2*h(t) - t)*(diff(g(t),t)*f(t) + diff(f(t),t)*g(t))-f(t)*g(t)
+ (h(t) - 3/2*h(t))*f(invg(h(t)))*(diff(g(t),t)*f(t)+diff(f(t),t)*g(t));
 eq4:= h(invh(t))-t;
 eq5:= g(invg(t))-t;
ord:= 8; # order for the series
eqs8:= series~(eval([eq1,eq2,eq3,eq4,eq5],
[f=unapply(add(F[i]*t^i,i=0..ord),t), g=unapply(add(G[i]*t^i,i=1..ord),t),
h=unapply(add(H[i]*t^i,i=1..ord),t), invg=unapply(add(GI[i]*t^i,i=1..ord),t),
invh=unapply(add(HI[i]*t^i,i=1..ord),t)]),t,ord);
solve(coeffs~(convert~(%,polynom),t));


Interestingly enough, although this does get me a solution to order 8, if I change ord to 9 there is no solution.  I found that rather surprising: usually if something like this works for low orders, it will also work to arbitrary orders.  But it appears this is a true result, not a flaw in the Maple code: there is no solution as formal power series in t with g(0) = h(0) = 0.

@Alejandro Jakubi : Curiouser and curiouser!  Are we using different versions of Maple 15.01?
Mine (in Standard GUI, 32 bit, under Windows 7) is build ID 636299, dated June 2, 2011, and showstat(`.`,1..3) produces

`.` := proc()
local s, t, i, const, const1, f, hasrtable, b1, b2, opts;
   1   if nargs = 2 and type([args],['rtable', 'rtable']) then
   2     try
   3       return LinearAlgebra:-Multiply(args,_fromdot)
         catch :
           ...
         end try
         ...
       end if;
       ...
end proc

In Maple 15.01, I get the same result with LinearAlgebra:-Multiply(A,B) as I get with A.B.  However,  strangely enough, LinearAlgebra:-Multiply(A,B, _fromdot) gives the correct answer.
This is really weird, because according to the code for `.`, A . B (with A and B being rtables) should call LinearAlgebra:-Multiply(A,B, _fromdot).  But apparently, the command A . B does not cause the code for `.` to be executed at all in Maple 15.01 (it would in Maple 14).

Certain automatic simplifications are applied to all output, and are very hard to avoid.  This includes taking out that factor of 1/3.  Well, in a pinch you could do this:

> (5*x+4)/Typesetting:-Typeset(3*x) = 7;

First 6 7 8 9 10 11 12 Last Page 8 of 187