@spark1631 

I'll get to your example later (honest).

I made the throwaway remark

whilst this approach might work for the trivial problem which you gave, it is doomed to failure whenever you have more than one abs() term in your expression(s).

You stated that you wanted to work with

multiple absolute term in the equations? Is there a general way to get what I need?

and this is where you are doomed to failure. Suppose you have a system of three equations in three unknowns, with one of the equations being |x-2|+|y-3|=z: it is simply not possible to remove the abs() terms by "squaring".

The only types of equations where your approach is possible are

1. the equation contains a single isolatable abs() term
2. the equation has two abs terms, and no others, ie can be written as abs(f(x))=abs(f(y))

so the squaring approach works for an incredibly restricted set of trivial equations and fails everywhere else

Back to your example

sys := {x = abs(y-1), 5^(2*x+2) = 25^(2*x-3), 5^(3*x+1)/5^(x-1) = 25^(2*x-3)};

This has several peculiarities: my first observation is that the second and third equations are identical, which you can confirm by

Psys:=map(simplify, sys)

which will return a set of two equations (because sets cannot contain two identical elements). The next observation is that these equations are not coupled, because the x-value can be obtained by

Qsys:=isolate(Psys[2],x);

and then the y-value by

Rsys:=subs(Qsys, Psys[1]);
solve(Rsys);

Since the solution is trivial to find, I tried to isolate what might be going wrong. My thought process is outlined in the attached worksheet

teq.mw

My tentative conclusion is that you have inadvertently uncovered a bug, in that the expression

testeq(5^(2*abs(y-1)+2) = (1/15625)*625^abs(y-1))

errors catastrophically. Throw your equations away and try this on its own.

I want to think about this a little more before I maybe reprot it as a bug, but I have to go do something else right now

Call me cynical but I think if you ask this question on a Maple forum, then the answer will be that Maple is better.

On the other hand if you ask this queation on a Matlab forum then the answer will be that Matlab is better

Realistically both Maple and Matlab appear to use a suite of linear algebra routines called lapack/blas. These are open-sourcer/free whatever the current terminology is, so no paid-for software (maple/Matlab) wants to admit that they are charging customers for including packages which is essentially free. (They will claim that they are charging for the "interface" to such software, or some similar weaselly comment).

Given that both products use the same underlying code, why should ther be any significant differences???? Beats me!

I suspect that the only way you can get a realisitic answer to this question is the old-fashioned method. Come up with a suite of linear algebra problems, and try the in your cnadidate software packages. Probably a good idea to submit tem to a matlab forum and separately to a Maple forum, then just sit back and watch who wins!!!

@Alas 

For the equations and bcs you supplied, answers are
0., 0.3320573, 15.02141
0.02, 0.3588413, 14.79601
0.04, 0.3839690, 14.60015
0.06, 0.4077149, 14.42746

Try the attached if you don't believe me

tp3.mw

@Markiyan Hirnyk 

I was just working from the manual - I never actually tried my suggestion. The only two things I can think to try now would be

1. (temporarily) remove the look-up method (the B[x[j]] thing) and go back to your original way of doing this part of the calculation - see if it makes any difference
2. GlobalOptima accepts an initialpoint as one of the options - so try giving it an explicit initalpoint which satisfies the problem condition add(x[j]^20, j = 1 .. 20) >= add(x[j]*10^(j-1), j = 1 .. 20). Can't be hard to find one by inspection!!

Is there any way of determining the initialpoint which was used when you had the nonnegint option????

@Markiyan Hirnyk 

Well Sergeii Moiseev is one smart guy - unfortunately I loaded this package a while ago  -and I probably forgot the author's name about 10mins after I loaded it. It happens when you are as old as I am.

Slightly puzzled by bad index into Array problem.

B[] is an array with indices 0,1,2,..9

x[j] is a digit between 0 and 9

therefore B[x[j]] always ought to return the 20th power of the digit x[j]. I can only think of two ways that one can fail with this error

1. although x[j] is constrained to be nonnegint, nothing stops it being 11 or 12 or whatever - and that will break the call to B[]. Doesn't actually seem too unlikely :-(
2. possibly the nonnegint constraint is applied at the start/end of each iteration, but *possibly* at an intermediate point in the iteration it is allowed to be non-integer, which would also break the call to B[]: this one seems a lot less likely than (1) above

According to the help, GlobalOptima() variables are among the options to be passed to the Search command, and the Search() command allows variables to be defined as discrete - as in

The k-dimensional discrete variables can be given as [x1, x2, ..., xk] = M within constr, where x1, x2, ..., xk is variable names, M is one of the  listlist, Matrix containing feasible k-dimensional points for the multi-dimensional discrete variables

My understanding of this option is that rather than using nonnegint you would set up a 20*10 matrix where each row would be 0,1,2,...9: this simple matrix therefore contains all possible values for all x[j], j=1,..20. Using this option ought to restrict the search space a lot more than the the simple nonnegint constraint, and also handle item (1) in the above list, so that lookups could be used. Both of thes should speed the problem up considerably

It should have occurred to me that the easiest way to generate a "smooth" curve is not to try interpolating values from the five points in your original list for L. Rather just generate a lot more L-values and solve the problem for all of them.

The code in the attached worksheet (minor revisions from my previous version) solves your problem for L=2.5..10 in steps of 0.1, and becuase one now has ~75 points, the resulting curves look "smooth"

trickPlot.mw

Of course because the differential equations are now being solved about 75 times rather 5 times, it is a little slower, but still runs in about 1sec on my machine

@Markiyan Hirnyk 

Sorry Markiyan, don't understand that last comment at all. Who is Sergei Moiseev????

I was just trying to make the point that you used the command

DirectSearch:-GlobalOptima(add(x[j]*10^(j-1), j = 1 .. 20), {add(x[j]^20, j = 1 .. 20) >= add(x[j]*10^(j-1), j = 1 .. 20), seq(x[k] <= 9, k = 1 .. 20)}, assume = nonnegint, maximize);

and that evaluations *might* be speeded up if you had used something like (untested code for illustration only)

B:=Array(0..9, [seq(j^20, j=0..9)])
DirectSearch:-GlobalOptima(add(x[j]*10^(j-1), j = 1 .. 20), {add(B[x[j]], j = 1 .. 20) >= add(x[j]*10^(j-1), j = 1 .. 20), seq(x[k] <= 9, k = 1 .. 20)}, assume = nonnegint, maximize);

ie substitute a lookup for an evaluation

My comment was not based on any criticism of the DirectSearch package, just a method of speeding up the function evaluations whcih the DirectSearch package needs.

@Carl Love 

I'm with Carl - if the last entry represents the number of function evaluations then it can't take "a few hours" to do 2921 evaluations - just can't!

The only other speed-up which springs to mind is not to compute the ^20 in the function evaluations. The obvious alternative is to create an array like

B=Array( 0..9, [seq(j^20, j=0..9)])

and use lookups in the function evaluations - after all x[j] is always a member of  [0..9] so x[j]^20 would simply be B[x[j]], and B[x[j]] is going to be a lot faster that x[j]^20!!  If the 2921 function evaluations really are taking this long, then using lookups ought to speed things up considerably

I'm waaaaay out of my comfort zone here, but a lot of the DirectSeacrh commands use a "mesh" of (I think) 500 points, and I think the "mesh gets refined at each iteration - so maybe the 2921 refers to the number of mesh computations, which would jack the number of calculations up to 2921*meshpoints=1.5Million (roughly)- although I still wouldn't have expected 1.5M evaluation to take "hours".

@momoklala 

Well if you change the last but one command in the worksheet I sent previously to

ps:= seq( plot( (L, convert(resArr[j,2..6],list)),
                  color=colors[j+1],
                  style=line,
                  legend=typeset("eta=", resArr[j,1]),
                  legendstyle=[font=[Times, bold, 16]] ), j=0..0):

it will produce a plot of theta(0) versus L.

You only have five points to plot - it will never look "smooth". You could I suppose try using CurveFitting[ArrayInterpolation] to "fake" a smooth curve

@Markiyan Hirnyk 

If you have 20 digits, shouldn't the number be constructed as

add(x[j]*10^(j-1), j=1..20)

in order that x[1] represents 0...9, not 10,20,etc - or am I missing something??

@Markiyan Hirnyk 

Might also (instead?) want to try the DirectSearch package. I've had some surprising successes with this (somewhat offset by a few irritating buglets)

So to summarise, for a 40-element vector, with two storage options, I have verified that

                                                     storage=sparse                 storage=rectangular

add(p),                                              doesn't work                               works                
add(j, j in p),                                     doesn't work                               works
add(p[j], j=1..40),                                  works                                   works
add(p[1..40]),                                   doesn't work                              works
add(p[1..-1]),                                    doesn't work                              works
add(p[..]),                                         doesn't work                              works
add(seq(j, j in p )),                                  works                                  works
add(seq(p[j],j=1..40)),                            works                                  works

Results for mul() are identical.

It is *almost* true that problems are restricted to the 1-element add(): the exception being the 2-element
add(j, j in p)

Am I the only one who finds this unacceptable????

@Markiyan Hirnyk 

My speed-up point was addressed specifically to the issue of solving this problem using enumeration.

I completely accept that if the exponents are big enough, then enumeration is a bad idea.

Considering your near-identical previous question was

30a+75b+110c+85d+255e+160f+15g+12h+120i=800000

I think you should clarify: are you trying to win an award for stupidity or idleness?

I'm not exactly sure that I understand what you are trying to do, but the attached code produces the "correct" answers for your test case, and ought to handle more than two equations, and equations of higher order, although I haven't checked these two cases.

polyChk.mw