John May

Dr. John May

3026 Reputation

18 Badges

17 years, 320 days
Maplesoft
Pasadena, California, United States

Social Networks and Content at Maplesoft.com

Maple Application Center
I have been a part of the Mathematical Software Group at Maplesoft since 2007. I have a Ph.D in Mathematics from North Carolina State University as well as Masters and Bachelors degrees from the University of Oregon. I have been working on research in computational mathematics since 1997. I currently work on symbolic solvers and visualization as well as other subsystems of Maple.

MaplePrimes Activity


These are replies submitted by John May

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

As I understand it, subs will do a traversal of the full expression as a tree while eval is more careful to traverse the expression as a DAG.  With Maple 12:

foo:=expand((sin(v)+cos(u)+cos(w)+2*sin(w)+cos(v)+1)^20):
bar := add((foo+i)^2, i=1..200):
time(subs({u=0, v=0, w=0}, bar));
                                    30.717

memory used=893.6MB, alloc=809.4MB, time=51.10

And in a new session:

foo:=expand((sin(v)+cos(u)+cos(w)+2*sin(w)+cos(v)+1)^20):
bar := add((foo+i)^2, i=1..200):
time(eval(bar,{u=0, v=0, w=0}));
                                    22.685

memory used=588.2MB, alloc=573.5MB, time=41.66

John

As I understand it, subs will do a traversal of the full expression as a tree while eval is more careful to traverse the expression as a DAG.  With Maple 12:

foo:=expand((sin(v)+cos(u)+cos(w)+2*sin(w)+cos(v)+1)^20):
bar := add((foo+i)^2, i=1..200):
time(subs({u=0, v=0, w=0}, bar));
                                    30.717

memory used=893.6MB, alloc=809.4MB, time=51.10

And in a new session:

foo:=expand((sin(v)+cos(u)+cos(w)+2*sin(w)+cos(v)+1)^20):
bar := add((foo+i)^2, i=1..200):
time(eval(bar,{u=0, v=0, w=0}));
                                    22.685

memory used=588.2MB, alloc=573.5MB, time=41.66

John

This tripped me up not too long ago, so I feel like I should point out that the subs and eval calls Mariner gives are not exactly equivalent,   the eval will do similtaneous substitution while the subs call does sequential substitution:

subs( x=y, y=x, [x,y] );
                                    [x, x]

eval([x,y], {x=y, y=x});
                                    [y, x]


The following are equivalent:

# Simultaneous
subs( [x=y, y=x], [x,y] );
                                    [y, x]

eval([x,y], {x=y, y=x});
                                    [y, x]

and

# Sequential
subs( x=y, y=x, [x,y] );
                                    [x, x]

eval(eval([x,y], x=y), y=x);
                                    [x, x]

Also, on complicated expressions, the calls to eval will generally be faster.

John

This tripped me up not too long ago, so I feel like I should point out that the subs and eval calls Mariner gives are not exactly equivalent,   the eval will do similtaneous substitution while the subs call does sequential substitution:

subs( x=y, y=x, [x,y] );
                                    [x, x]

eval([x,y], {x=y, y=x});
                                    [y, x]


The following are equivalent:

# Simultaneous
subs( [x=y, y=x], [x,y] );
                                    [y, x]

eval([x,y], {x=y, y=x});
                                    [y, x]

and

# Sequential
subs( x=y, y=x, [x,y] );
                                    [x, x]

eval(eval([x,y], x=y), y=x);
                                    [x, x]

Also, on complicated expressions, the calls to eval will generally be faster.

John

Sorry, I spoke too soon.

f(x,a,b) assuming a::integer, b::integer;

does indeed return a result invalid for b>1. Using the other seems to give the expected value:

g := f(x,a,b) assuming a::positive, b::positive;

                        a
                       x  hypergeom([a, -b + 1], [1 + a], x)
                  g := -------------------------------------
                                         a
 (expand@simplify)(eval(g, [a=2, b=3]) ); 
                                2        3        4
                           1/2 x  - 2/3 x  + 1/4 x

Clearly something wierd is happening in both cases, however.

John

Sorry, I spoke too soon.

f(x,a,b) assuming a::integer, b::integer;

does indeed return a result invalid for b>1. Using the other seems to give the expected value:

g := f(x,a,b) assuming a::positive, b::positive;

                        a
                       x  hypergeom([a, -b + 1], [1 + a], x)
                  g := -------------------------------------
                                         a
 (expand@simplify)(eval(g, [a=2, b=3]) ); 
                                2        3        4
                           1/2 x  - 2/3 x  + 1/4 x

Clearly something wierd is happening in both cases, however.

John

value( Int( t^(a-1) * (1-t)^(b-1), t=0..x ) );

is exactly equivalent to calling

int( t^(a-1) * (1-t)^(b-1), t=0..x );

The point of setting

f := (x,a,b) -> Int( t^(a-1) * (1-t)^(b-1), t=0..x );
                                    x
                                   /
                                  |    (a - 1)        (b - 1)
               f := (x, a, b) ->  |   t        (1 - t)        dt
                                  |
                                 /
                                   0

would be work around the problem in int to avoid computing the integral with symbolic parameters a and b and instead to only compute after values were assigned:

v1:=value( f(2,3,2/5) ); v2:=evalf( f(2,3,.4) );

In Maple 12:

                                                 2/5
                                  125   365 (-1)
                            v1 := --- - -----------
                                  84        84

                      v2 := 0.1453428221 - 4.132567005 I

If you are only interested in floating point values this is probably faster (v2 is computed very differently than v1).

John

value( Int( t^(a-1) * (1-t)^(b-1), t=0..x ) );

is exactly equivalent to calling

int( t^(a-1) * (1-t)^(b-1), t=0..x );

The point of setting

f := (x,a,b) -> Int( t^(a-1) * (1-t)^(b-1), t=0..x );
                                    x
                                   /
                                  |    (a - 1)        (b - 1)
               f := (x, a, b) ->  |   t        (1 - t)        dt
                                  |
                                 /
                                   0

would be work around the problem in int to avoid computing the integral with symbolic parameters a and b and instead to only compute after values were assigned:

v1:=value( f(2,3,2/5) ); v2:=evalf( f(2,3,.4) );

In Maple 12:

                                                 2/5
                                  125   365 (-1)
                            v1 := --- - -----------
                                  84        84

                      v2 := 0.1453428221 - 4.132567005 I

If you are only interested in floating point values this is probably faster (v2 is computed very differently than v1).

John

First 13 14 15 16 17 18 19 Page 15 of 19