Thanks acer, great!
This helps tremendously and will enable me to move forward with this simulation. I can also see how this will be useful in other problems I will be looking at soon.
There are 3 distinct things in your post:
- use eq as a constraint, rather than maximizing wrt the implicit solution
- initialize the problem with 'initialpoint'
- convert powers to the surd function
When stating the problem I described using eq as a constraint (point 1 above) but somehow got lost along the way trying to handle RootOf directly and missed it. On a better day I might have got there. After another day of searching through the help pages and the mapleprimes archive, I might (ought to) have been able to figure out 'initialpoint' (point 2 above). However, point 3 above would have stumped me, because without the conversion to 'surd', the Optimization fails. So converting powers to surd is a crucial step in solving the problem, it would appear.
restart; # acer's solution without surd fails
### original equation
eq := x = 11/500+6*u*(1/25)-u^2+(-174200000000*x^5+13045131*E^4*u^4)^(1/5)/(65225655*E^4*u^4):
Error, (in Optimization:-NLPSolve) complex value encountered
A natural follow-up question, therefore, is: How can I convert the powers to surd? These equations are the result of prior computations, so I can't just "input" them with surd written in, I need to convert them.
The powers appear in the middle of the expression, rather than around it, so I need to instruct convert(expr,surd) to search inside the expression for instances of `^`. How?
convert(eq,surd); # FAILS
convert((-174200000000*x^5+13045131*E^4*u^4)^(1/5),surd); # WORKS when the power wraps
hastype(eq,`^`); # Maple can "see" the powers inside the expression
Also, what if the powers are rational, e.g. 3/5 instead of 1/5?
The help pages suggest surd takes an integer:
x - any algebraic expression
n - any algebraic expression, understood to be an integer
Example with a rational power:
eq := x = 11/500+6*u*(1/25)-u^2+(-174200000000*x^5+13045131*E^4*u^4)^(3/5)/(65225655*E^4*u^4):
If this is a mess, well never mind, I'll restrict my parameter values accordingly.
Last remark, on accessing the solutions. Previously, I was using things like op([1,1],sol), but you took a different approach, so I reverse-engineered your code and came up with this:
### how I usually do this
xpts := [ seq( [ op([i,1],[sols]), op([i,2,1],[sols]) ], i = 1 .. nops([sols]) ) ];
upts := [ seq( [ op([i,1],[sols]), rhs(op([i,2,2,1],[sols])) ], i = 1 .. nops([sols]) ) ];
### alternative based on acer's code
xpts := [ seq( [ s, s[2,1] ], s in sols ) ];
upts := [ seq( [ s, rhs(s[2,2,1]) ], s in sols ) ];
The latter approach is clearly shorter and more elegant. Both approaches rely on the solution being ordered in a consistent way across sessions. To be "surer" that the ordering would be consistent I wrote 'initialpoint'=[u=0.5,x=1], so the u comes before the x, as it is with the arguments (u=0..1,x=0..1).