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 answers submitted by John May

This is a bit of a work around:  try calling radnormal.  It works on the first example:

> radnormal((-q0)^(-1/2)*(-k0)^(1/2)*q0*k0^2*(q0*k0)^(-1/2));
                                     1/2      2
                                (-k0)    q0 k0
                              -------------------
                                   1/2        1/2
                              (-q0)    (q0 k0)

> radnormal((-q0)^(-1/2)*(-k0)^(-3/2)*q0*k0^4*(q0*k0)^(-1/2));
                                     1/2      2
                                (-k0)    q0 k0
                              -------------------
                                   1/2        1/2
                              (-q0)    (q0 k0)

On the second example, I  called simplify(<foo>,'exp'); first (simplifies just the exp/log terms), then radnormal, then a full simplify:

> simplify(radnormal(simplify(exp(-1/2*ln(-q0)+1/2*ln(-k0))*q0*k0^2/(q0*k0)^(1/2),exp))) assuming q0<0,k0<0;
                                        2
                                     -k0

> simplify(radnormal(simplify(exp(-1/2*ln(-q0))*exp(1/2*ln(-k0))*q0*k0^2/(q0*k0)^(1/2),exp))) assuming q0<0,k0<0;
                                        2
                                     -k0

 

John

If you do want to get the polynomial as a vector or list of its coefficients, there is command to do that as well:

?PolynomialTools,CoefficientList

John

This is the system as rendered by my cut and paste and plain text conversion of sb572's post.

sys:=
{
x[1]^2+y[1]^2-1 ,
x[2]^2+y[2]^2-1 ,
x[3]^2+y[3]^2-1 ,
x[4]^2+y[4]^2-1 ,
x[5]^2+y[5]^2-1 ,

x[1]*x[2]+x[1]*x[3]+y[1]*y[3]+y[1]*y[5]+x[2]*x[5]+x[1]*x[5]+x[3]*x[4]
+y[3]*y[5]+y[2]*y[3]+x[3]+x[4]+x[1]*x[4]+x[3]*x[5]+y[1]*y[4]
+y[1]*y[2]+x[4]*x[5]+x[5]+y[3]*y[4]+y[2]*y[4]+x[1]+x[2]*x[4]+x[2]+
x[2]*x[3]+y[2]*y[5]+y[4]*y[5],

-y[3]*sqrt(3)*x[4]+2*x[1]*y[3]*sqrt(3)-2*y[1]*x[3]*sqrt(3)+x[3]*y[4]*sqrt(3)
+y[1]*x[5]*sqrt(3)-2*y[2]*sqrt(3)*x[5]+y[4]*sqrt(3)*x[5]+y[2]*sqrt(3)*x[3]
-x[2]*y[3]*sqrt(3)+y[1]*sqrt(3)+3*y[1]*y[5]+3*x[1]*x[5]
-x[4]*y[5]*sqrt(3)+3*x[3]*x[4]+3*y[2]*y[3]+2*x[2]*y[5]*sqrt(3)
-2*y[4]*sqrt(3)+3*x[4]*x[5]+3*y[3]*y[4]+3*x[1]-x[1]*y[5]*sqrt(3)
+3*x[2]+y[2]*sqrt(3)+3*x[2]*x[3]+3*y[4]*y[5],

-y[3]*sqrt(3)*x[4]+x[3]*y[4]*sqrt(3)+y[1]*x[5]*sqrt(3)-y[2]*sqrt(3)*x[5]
-y[4]*sqrt(3)*x[5]+y[2]*sqrt(3)*x[3]-y[1]*x[2]*sqrt(3)-x[3]*y[5]*sqrt(3)
+x[1]*y[2]*sqrt(3)-x[2]*y[3]*sqrt(3)+y[1]*sqrt(3)+x[4]*y[5]*sqrt(3)
-3*x[1]*x[3]+x[2]*y[5]*sqrt(3)-y[4]*sqrt(3)+y[3]*sqrt(3)-3*x[5]-3*y[1]*y[3]
-3*y[2]*y[4]+x[1]*y[4]*sqrt(3)-x[1]*y[5]*sqrt(3)-y[2]*sqrt(3)-3*x[2]*x[4]
+y[3]*sqrt(3)*x[5]-y[1]*x[4]*sqrt(3),

y[2]*sqrt(3)*x[4]-x[1]*y[3]*sqrt(3)+y[1]*x[3]*sqrt(3)-y[4]*sqrt(3)*x[5]
+y[2]*sqrt(3)*x[3]-2*x[3]*y[5]*sqrt(3)-x[2]*y[3]*sqrt(3)+y[1]*sqrt(3)
+y[5]*sqrt(3)+x[4]*y[5]*sqrt(3)-3*y[2]*y[3]-3*y[4]*y[5]-3*x[1]*x[3]
-3*x[2]*x[3]-x[2]*y[4]*sqrt(3)-3*x[5]-3*y[1]*y[3]-3*x[1]-3*y[2]*y[4]
+2*x[1]*y[4]*sqrt(3)-2*y[2]*sqrt(3)-3*x[4]*x[5]-3*x[2]*x[4]
+2*y[3]*sqrt(3)*x[5]-2*y[1]*x[4]*sqrt(3),

y[2]*sqrt(3)*x[4]+x[1]*y[3]*sqrt(3)-y[1]*x[3]*sqrt(3)+2*y[4]*sqrt(3)*x[5]
+y[5]*sqrt(3)-2*x[4]*y[5]*sqrt(3)-3*y[2]*y[3]-3*y[4]*y[5]-3*y[3]*y[5]
-3*y[2]*y[5]-3*x[3]-2*y[3]*sqrt(3)-3*x[4]-3*x[2]*x[3]-x[2]*y[4]*sqrt(3)
-2*x[1]*y[2]*sqrt(3)-3*y[1]*y[4]-3*x[1]-3*x[1]*x[2]+2*y[1]*x[2]*sqrt(3)
+y[3]*sqrt(3)*x[4]-3*x[3]*x[5]+y[2]*sqrt(3)-3*x[1]*x[4]-3*x[4]*x[5]
-3*x[2]*x[5]+x[1]*y[5]*sqrt(3)-y[1]*x[5]*sqrt(3)-x[3]*y[4]*sqrt(3)-3*y[1]*y[2]
};

Maple 12's solve throws FGb at this system and it runs out of memory trying to get a tdeg Gröbner basis.  Maple 11's solve uses different techniques and hasn't crashed on it after about an hour and looks doomed to churn for as long as I let it.   So, the obvious things don't work...

John

Of course, identify works beautifully if you cheat:

flts:=[seq]( (evalf@eval)(k/Pi^4,k=i),i=100..110):
map(identify, flts, BasisPolyConst=[Pi^4]);

        100  101  102  103  104  105  106  107  108  109  110
       [---, ---, ---, ---, ---, ---, ---, ---, ---, ---, ---]
          4    4    4    4    4    4    4    4    4    4    4
        Pi   Pi   Pi   Pi   Pi   Pi   Pi   Pi   Pi   Pi   Pi

 

John

map and seq can be a lot faster that doing a for-do loop, but not in this case:

> L:=[seq](2*i+1, i=10000..100000):
> time(seq( sqrt(5*L[i]- L[i+1]),i=1..nops(L)-1) );
                                     2.316
> restart;
> L:=[seq](2*i+1, i=10000..100000):
> t:=time(): for i from 1 to nops(L)-1 do          
>    sqrt(5*L[i]-L[i+1]):                            
> end do: time() - t;
                                       1.880

seq and map are best in the case when you are creating a list as output. So, the following is extremely slow (I stopped it at 400 seconds)

> restart;
> L:=[seq](2*i+1, i=10000..100000):
> M:=[];
> t:=time(): for i from 1 to nops(L)-1 do
> M := [op(M), sqrt(5*L[i]-L[i+1])]:                            
> end do: time() - t;

But that is mostly due to the bad list construction technique, not the for loop. The following isn't too bad:

> restart;
> L:=[seq](2*i+1, i=10000..100000):
> M := Array(1..(100000-10000), 0);
> t:=time(): for i from 1 to nops(L)-1 do   
>    M[i]:=sqrt(5*L[i]-L[i+1]):                
> end do: M := convert(M, list): time() - t;
                                     2.792

John

Even with k constant, this type of system is hard for solvers (Maple just returns a RootOf that allvalues doesn't expand as Axel points out and other computer algebra systems return similar things).    Solvers can give symbolic answers to some exp-trig expressions but generally only if the answer can be found by replacing the exp and trig terms with symbols (i.e. if the solution is fundamentally algebraic and not analytic). 

John

If you look closely at the output of pdesolve,

> diff(F(x,y),x)*(x+ y)+diff(F(x,y),y)*(4*x + y):
> pdsolve(%);                                    
                                             1
                      F(x, y) = _F1(--------------------)
                                             3
                                    (y - 2 x)  (y + 2 x)

> lprint(%);
F(x,y) = _F1(1/(y-2*x)^3/(y+2*x))

you'll see that _F1 is a function and not a symbolc multiplied by the rest of the expression.  This is why solve can't do anything for you. A look at the help page for pdsolve suggest that _F1 is actually an arbitrary function:

All functions introduced by pdsolve beginning with _F and followed by a number are assumed to be arbitrary, sufficiently differentiable functions of their arguments.

 

 

John

 Take a look at the ImageTools package, in particular, the Read function.  For an example of using ImageTools, you can also take a look at this blog post.

John

 Looks like these messages are being printed on the 'hints' infolevel. 

infolevel[hints] := 0;

should turn them off.

John

Since you are looking for a floating point answer, and all your expressions have floating point numbers, you should really use  fsolve instead of solve.  That is, replace T= solve(f=0,T); with T = fsolve(f,T);

If you use solve, your f is converted to (notice no floating point values):

f:=-1+657894737/500000000000*10^(40356350000-173063/100/(T+116713/500));

                                  /                   173063     \
                                  |40356350000 - ----------------|
                                  |                  /    116713\|
                                  |              100 |T + ------||
                                  \                  \     500  //
                      657894737 10
            f := -1 + --------------------------------------------
                                      500000000000

which Maple tries to solve exactly (it then calls evalf on the answer in the end so you don't notice), and that appears to cause big expression blow up.

John

You can of course do this with ListTools as well:

with(ListTools):

L:=[a,b,c,d,c,d,e,f,a]:
MakeUnique(FindRepetitions(L));
	[c,d,a]

FindRepetitions removes one copy of each unique element which is why you get the ordering of the elements you see here.

Lots of useful stuff in ListTools.

John

Strangely, even the output of Groebner[Solve] does not get solve() to return a solution even though Groebner clearly says the system is two dimensional.

I was able to use triangular decomposition to get an answer however:
 

with(RegularChains):
R:=PolynomialRing([x2,x3,x5,x8,x6,x7,x9,y2,y1,x1,x4]);
sol:=Triangularize(sys, R);
solve(Equations(sol[1], R));

This gave me the following solution with x8 and y2 free:

      15751       7813
{x3 = -----, x2 = ----, x8 = x8, y2 = y2,
       125        125

         -7812               15876
    x1 = -----, y1 = 0, x4 = -----, x9 = RootOf(
          125                 125

               2             2
    15625 y2 _Z  + (-15625 y2  - 61027344) _Z

                       15876
     - 15625 y2), x7 = -----, x5 =
                        125

             2                               2
    125 y2 x8  - 15876 y2 + 15876 x8 - 125 y2  x8
    ---------------------------------------------,
                       125 x8

               2
         125 x8  - 15876
    x6 = ---------------}
             125 x8

John May
Developer, Maplesoft

Setting the environment variable _MaxSols can force solve to return only one solution:
_MaxSols:=1;

This can also make solve run faster internally in some, but certainly not all, cases.

John

One posibility is to avoid simplifying entirely.  In this case, you can compute the symbolic matrix "Du" by turning off normalization with:

Normalizer:=x->x; # Tells LinearAlgebra not to call normal() on intermediate results
Testzero:=testeq; # Tells LinearAlgebra to use the testeq(foo) instead of evalb(Normalizer(foo)=0)

Du := -(G21.(G11-Pt)^(-1).G12-G22)^(-1): # <- do not try to print this result

But now, Du is extremely big:

> map(length, Du)[1,1];

53574525

It would be a very bad idea to try to print or simplify Du, but it is diagonal:

> map(testeq, Du);

	[ false  true   true  ]
	[ true   false  true  ]
        [ true   true   false ]

If you are careful, you might be able to do something with this result, but since it is so large, it is going to be tough.

John

Assuming eq1 is expanded out into sums of the form c*y^n you could do

pp, c := selectremove(has, [op(eq1)], y); # c is a list of constant terms
pw, pd := selectremove(type, pp, `^`); # pd is products, pw is powers
pd, lin := selectremove(type, map2(op, -1, pd), `^`); # lin is the linear terms in y
[op(map(0,c)), op(map(1, lin)), op(map2(op, 2, pw)), op(map2(op, 2, pd))];

Of course, unlike Robert's solution, this won't preserve the order of the exponents.


First 7 8 9 10 11 Page 9 of 11