John May

Dr. John May

2351 Reputation

17 Badges

12 years, 308 days
Maplesoft
Guru
Pasadena, California, United States

Social Networks and Content at Maplesoft.com

Maple Application Center

I am a Senior Developer in the Mathematical Software Group and have been with Maplesoft since 2007. I am also an Adjunct Assistant Professor in the School of Computer Science at the University of Waterloo.

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.

My main research interests in are computational linear and polynomial algebra, especially numerical polynomial algebra. I currently work on the exact algebraic solvers as well as other subsystems of Maple.

MaplePrimes Activity


These are answers submitted by John May

The first problem is your parentheses.  Get rid of the one before "int" and get rid of one of the four right before ",t=0.."and that should get rid of the error.

John

This sort of thing can happen if there are no unconditional solutions. That is, the solutions depend on the values of the parameters.  At some point we will get to beefing up solve for parameters, but that is fundamentally difficult and thus is a big project.

In the mean time, here are a couple of things to try:

# solve with variables in a list, variables first, followed by parameters.
sol := solve({eq},[p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11, 
    pi01, pi00, pi10, pi11, sigma10, sigma01, sigma11, sigma00]);


[
p0 = -p1 + pi01, 
p1 = p1,
p2 = pi11 + sigma11 + sigma10 - pi01 - sigma01 - sigma00, 
p3 = 0,
p4 = -sigma11 - sigma10 + sigma01 + sigma00, 
p5 = 0,
p6 = 1 - pi11 - sigma11 - sigma10,
p7 = 0,
p9 = -sigma01 + sigma11, 
p10 = -p11 + sigma10,
p11 = p11,
# omitting equations like parameter=parameter
pi00 = -pi01 + 1 - sigma01 - sigma00,
pi10 = 1 - pi11 - sigma11 - sigma10
]

Notice.  p1 and p11 are free variables, and there are only solutions if parameters satisfy the last two equations.

Or write as a matrix:

eq := {eq}; vars := {p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11};
M:=LinearAlgebra:-GenerateMatrix(eq, vars, augmented=true):
# The REF:
LinearAlgebra:-GaussianElimination(M);

# Looking at the last two columns:
# 1. the solutions space is two dimensional
# 2. there are solutions only if the parameters satisfy:
0 = sigma00 - sigma10 - pi11 + pi01 + pi00 - pi10 + sigma01 - sigma11

(notice this equation is the sum of the last two equations that solve gave)

John
Math Software, Maplesoft

 Double integrals can also be given with the ranges as a list in a single Int call:

Int(1,[rbar=0..RootOf(f(r,z)-1,r),z=0..1]);

evalf will sometimes be able to compute things in a smarter fashion when iterated integrals are given this way.

John

If you are trying to grab the kth decimal digit, you could use mod like this:

f := (x,k) -> ((x mod 10^k) - (x mod 10^(k-1)))/10^(k-1):
> f(345792,2);

                                       9

> f(345792,5);

                                       4

John

Are these homework problems from a course?  If so, in the future, you should be posting in the student forums.

John

This sort of thing is hard to do completely automatically in general.

algsubs seems like it should do the first one directly, but it does not:

> algsubs((x+1)^2 = a, (x+1)^6+2);      
                                        6
                                 (x + 1)  + 2

But, you can stack calls to do it:

> subs(z=x+1,algsubs((z)^2 = a, algsubs(x+1=z, (x+1)^6+2)));
                                         3
                                    2 + a

For the second one, it is trickier, you could probably use subsindets first:

> subsindets((x^2+2*x+1)^3+2, `^`, factor);
                                        6
                                 (x + 1)  + 2

 

John

Depending on how you call solve, you will either get RealRanges or inequalities for your output. If it would be more useful to have inequalities, enclose your variable in a set.  e.g.

> solve(x<2,x);      
                         RealRange(-infinity, Open(2))

> solve({x<2},{x});
                                    {x < 2}

If you really just need to process real ranges, you can use op:

> foo:=RealRange(Open(-1.345464), Open(4.5343535));
              foo := RealRange(Open(-1.345464), Open(4.5343535))

> op([1,1],foo);
                                   -1.345464
> op([2,1],foo);
                                   4.5343535

John

Indeed there is.  Take a look at: ?LinearAlgebra,SubMatrix

For simple examples you can also used ranged indexing.  For your examples: A[2..3,2..3];

John

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

First 6 7 8 9 10 Page 8 of 10