6237 Reputation

17 Badges

10 years, 105 days

MaplePrimes Activity

These are replies submitted by tomleslie

@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".


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


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.

@Markiyan Hirnyk 

Using my original code, if I increase the initial point count using the 'number' option, as in

DirectSearch[GlobalSearch](abs(eqn1), assume=posint, number=500):

I get 249 solutions - I suspect there are many, many more :-(

@Rouben Rostamian  

As my original response said, isolve() provides an accurate solution in terms of eight parameters and I actually reached the same point that you did, ie

5 _Z1 + 10 _Z2 + 25 _Z3 + 35 _Z4 + 25 _Z5 + 40 _Z6 + 5 _Z7 + 10 _Z8 <666400

in which _Z1, _Z2, _Z3, _Z4, _Z5, _Z8  are required to be positive (not arbitrary, as you stated) integers:  _Z6, and  _Z7 may be negative, but lower-bounded by the conditions that f,g must be positive, ie

 _Z6 >  -(1/3)*(2 + _Z3 + 2 _Z4)

 _Z7 > -(1/4)*(2 _Z1 + 3 _Z2 + 2 _Z3 + _Z4 + 3 _Z5)

Since one now has lower bounds on _Z1.._Z8, and upper bounds are trivial, eight for -loops to test the inequality and all solutions can be generated. However is this noticeably "better" than running nine for loops for the variables in the original equation???

The only part of my original post which I *really*wish I had stated differently is the sentence

Would I guarantee that these are *all* the posint solutions - well no probably not

which I would now rewrite as

Would I guarantee that these are *all* the posint solutions - almost certainly not

In addition to the well known problems of local/global optimisation, I suspect that the DirectSearch package has an upper limit on the number of solutions it returns (although I cannot find this documented anywhere)

Theoretically, isolve() should provide the  answer which you want  - and it does if you remove the constraint that all variables are positive integers.

However even the constraint-free solution is given in terms of eight parameters, which is accurate, but probably does not represent significant progress on the original problem. I personally failed to make much fiurther progress with the isolve() command: not saying it can't be done, just observing that I got nowhere.

My next attempt was to use the DirectSearch package, which is not technically part of Maple, but is a package which you can download (from maple's website and install). A little experimentation with this and I came up with

# rewrite the equation in a way which allows maxima and minima to be found (DirectSearch requirement)
# Perform a global search which *should* return all combinations of posint variables which
# produce either maxima or minima
   ans:=DirectSearch[GlobalSearch](abs(eqn1), assume=posint):
# Just because I'm paranoid about max/min/opt routines, perform a quick dirty check
# on whether the returned solutions are valid (and are minima), return only those which are
   validAns:=[seq(`if`(eval(eqn1, ans[j][2])=0,ans[j][2],NULL), j=1..ArrayTools[Size](ans)[1])];

The result of this code was to produce 44 combinations {a..i} of positive integers which satisfy the original equation.

Would I guarantee that these are *all* the posint solutions - well no probably not. This is not a reflection on the author(s) of the DirectSearch package, just a reflection of my experience being often caught out by local/global optimisation problems.

It would be possible to check for further solutionsby feeding the DirectSearch ones to the parameterised isolve() solution mentioned above, to see if any more crop up


Checked this in Maple 18 and Maple2015, and it produces 4.0000000 in both.


  1. Put a restart command immediately prior to the g:=Groe.....  and check whether the error still occurs
  2. Based on my tests the error will not occur - suggests that you have 'a' and/or 'b' defined to something a bit weird earlier in your worksheet



'sum' is for symbolic sums and it will do things like try to compute closed form solutions, (if it can) rather than just listing all the terms in the sum

'add' doesn't do any thing clever - it justs adds all the terms you give it

Respectfully suggest RTFM

@Thomas Richard 

As has been suggested


works (and keeps on working for more terms, higher powers). Didn't really attempt to find the limit doing it this way

Given Preben's comment I thought it might be possible to use the continuation method on the exp terms (see ?dsolve,numeric_bvp,advanced (bvp))

Having tried a couple of obvious "continuations", (including adjusting the continuation variable "manually") the *best* I managed to do was to convert the "Newton iteration is not converging" to "Error, (in dsolve/numeric/bvp) singularity encountered"

It may be possible to come with a continuation method which works, although I'm beginning to doubt it. I suggest you take a long hard look at the terms exp(beta/theta(eta)), particularly given the boundary condition theta(N) = 0

First 143 144 145 146 147 148 149 Last Page 145 of 157