It is doing a completely different computation:

sqrt(5*L[47]-L[48]): lprint(%);
3*8930^(1/2)

evalf(sqrt(5*L[47]-L[48]));
                                  283.4960317

John

 It is doing a completely different computation:

sqrt(5*L[47]-L[48]): lprint(%);
3*8930^(1/2)

evalf(sqrt(5*L[47]-L[48]));
                                  283.4960317

John

That a a great example of how to do the numerical computation very fast, but I should point out that it is not really a fair comparison.
For the exact operation this is about the same speed:

L:=[seq](2*i+1, i=10000..100000):
L:=Array(L):
L1:=ArrayTools:-Alias(L,[1..90000]):
L2:=ArrayTools:-Alias(L,1,[1..90000]):
time(map[inline](sqrt,5*L1+L2));
                                    2.708

John

That a a great example of how to do the numerical computation very fast, but I should point out that it is not really a fair comparison.
For the exact operation this is about the same speed:

L:=[seq](2*i+1, i=10000..100000):
L:=Array(L):
L1:=ArrayTools:-Alias(L,[1..90000]):
L2:=ArrayTools:-Alias(L,1,[1..90000]):
time(map[inline](sqrt,5*L1+L2));
                                    2.708

John

> L:=[seq](2*i+1, i=10000..100000):
> time(zip((x,y)->sqrt(5*x+y),L[1..-2],L[2..-1]));
                                     2.632
> restart
> L:=[seq](2*i+1, i=10000..100000):
> time(seq( sqrt(5*L[i]- L[i+1]),i=1..nops(L)-1) );
                                     1.676

(timings from a different machine than my earlier post)

John

> L:=[seq](2*i+1, i=10000..100000):
> time(zip((x,y)->sqrt(5*x+y),L[1..-2],L[2..-1]));
                                     2.632
> restart
> L:=[seq](2*i+1, i=10000..100000):
> time(seq( sqrt(5*L[i]- L[i+1]),i=1..nops(L)-1) );
                                     1.676

(timings from a different machine than my earlier post)

John

This is an automatic simplification that you can't do much about. Is it any comfort that it prettyprints nicely everywhere except prettyprint=0?

> a/(b*c);
                                       a
                                      ---
                                      b c
> lprint(%); # prettyprint=0 output 
a/b/c

> dismantle(%); # underlying DAG a^1*b^(-1)*c^(-1)

PROD(7)
   NAME(4): a
   INTPOS(2): 1
   NAME(4): b
   INTNEG(2): -1
   NAME(4): c
   INTNEG(2): -1

John

Interestingly, neither Maple 8 nor Maple 9 can solve this either. 

John

Interestingly, neither Maple 8 nor Maple 9 can solve this either. 

John

 Since _F1 is a function, you can substitute (almost) any function for it an have a solution.  The identity function is the simplest choice:
 

# sol is the output of pdsolve
id:=x->x:
sol1 := eval(sol, _F1=id):
expr := rhs(sol1);
                                          1
                          expr := --------------------
                                          3
                                 (y - 2 x)  (y + 2 x)

simplify(diff(expr,x)*(x+ y)+diff(expr,y)*(4*x + y));
                                       0

You could use use 'exp', or 'sin' instead of id, and it will also give a solution.  Anything sufficiently differentiable will work.  Your choice depends on your application I suppose.

John

 Since _F1 is a function, you can substitute (almost) any function for it an have a solution.  The identity function is the simplest choice:
 

# sol is the output of pdsolve
id:=x->x:
sol1 := eval(sol, _F1=id):
expr := rhs(sol1);
                                          1
                          expr := --------------------
                                          3
                                 (y - 2 x)  (y + 2 x)

simplify(diff(expr,x)*(x+ y)+diff(expr,y)*(4*x + y));
                                       0

You could use use 'exp', or 'sin' instead of id, and it will also give a solution.  Anything sufficiently differentiable will work.  Your choice depends on your application I suppose.

John

Very nice.

John

*edit: or that ---^

True, a one-liner would do it:

subsindets(sys, 'indexed', x->cat(op(0,x),'_',op(1,x),'_',op(2,x)) );

John

In fact it is a very subtle problem.  Indexed variables a[i,j] can cause problems in some solve subroutines, so one of the first things solve does is rename these variables to new symbols like x_nn.  In the case of this system, the order in which the polynomial system solver handles the variables matters a lot, and when the variables get renamed, they get ordered differently.  It seems more or less chance that the original variable names get ordered in a way so the polynomial system solver returns quickly.

Here is another work around (use symbols instead of indexed names):

vars := indets(sys);
vars2 := map(x->evaln(cat(op(0,x),'_',op(1,x),'_',op(2,x))), vars);
sys := subs([seq(vars[i]=vars2[i], i=1..nops(vars))], sys):
time(solve(sys));
                                    275.825

John

The solve problem reported in penkavmr's  original post in has been fixed.  In Maple 12, solve finds solutions for that system in less than one second. Try it:

time( solve(sys) );
                                     0.404

As far as the new system given by Alex, it looks like Maple versions 10-12 cannot solve it  (I stopped 10 and 12 each at 3000 seconds) while Maple 9.5 solves it in 200 seconds or so.  Even though it is a system similar to penkavmr's, the underlying problem appears to be completely different.

Also, in Maple 12, SolveTools[PolynomialSystems] does solve Alex's system (in about 250 seconds), and since solve pretty much immediately calls SolveTools[PolynomialSystems] in this case, that means the problem here is much weirder than the problem penkavmr's system was encountering (I suspect it is subtle interaction due to the renaming of the indexed variables). 

As for returning only one solution, that is the default, but SolveTools[PolynomialSystems] will return more if asked.  Take a look at its help page.

John May
Math Developer
Maplesoft