John May

Dr. John May

2351 Reputation

17 Badges

12 years, 72 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

If A is going to be between 0 and 1, you can get a better formula for the integral with:

E := int(1/(1+A*cos(2*Pi*x))^3, x = 0 .. 1, allsolutions) assuming A>0, A

Then you can generate the values with the ?seq command:

F := [seq(E, A = 0.1e-2.. 1, .1)];

 

John

Edit your file to be in the form:

edges := {{1,2},{1,3},...,{99999,99987},{99999,100000}};

then since it is now in Maple syntax, you just need to read it in with the ?read command.

read("C://Desktop//list.dat");

 

John

Check the help page: ?Groebner,Basis to see you aren't calling the Basis command exatly right.  The 'tord' argument directs to command to choose its own basis and then assign it to the variable 'tord' after it returns:

G1 := Basis({expand(w*I2[1]),expand(w*I2[2]), expand((1-w)*J2[1])}, 'tord'):
tord;
                         plex(y, x, w)

The last argument is actually ignored in that case.  I suspect that you want the following:

G2 := Basis({expand(w*I2[1]),expand(w*I2[2]), expand((1-w)*J2[1])}, lexdeg([w],[x,y]));

 

John

A RootOf() represents a generic/unspecified root of its argument.  So it does satisfy the expression in its argument.  For example:

p:=x^6+randpoly(x,degree=5,terms=3)+1;
r := RootOf(p,x);
eval(p, x=r);
evala(Normal(%)); #=0

Various other commands under ?evala allow you to compute with generic polynomial RootOfs: i.e. calculations that would be valid for any specific root of the equation.   You only need to use ?allvalues if/when you need specific values of the roots (for example, when you want to find all the positive roots).

John

Sometimes Optimization is able to do something with overconstrained systems. One way to try this is:

   read("Cvector_with_a.m"):
   infolevel[Optimization]:=4;
   Optimization:-Minimize(1,{seq(C[i]=0, i=1..9), c1>=1,c<=900,c2>=1,c2<=900,
c3>=1,c3<=900,c4>=1,c4<=900,c5>=1,c5<=900,c6>=1,c6<=900,c7>=1,
c7<=900,hb>=0.01,hb<=900, a1>=1.5, a1<=1.6, a2>=5,a2<=5.1},
feasibilitytolerance=10^(-3), iterationlimit=1000, initialpoint=[a1=1.52, a2=5.04,
c=1,hb=0.012,c1=46.9, c2=9.60, c3=7.15,c4=56.1,c5=18.9,c6=2.65,c7=7.74]);

Unfortunately, that does not appear to converge on a solution. You could play around with the various ?Optimization,Minimize options, the linear constraints on the variables, and the precision of the initial point.  You can also move the C[i]=0 constraints to the objective function and try that instead.  One of these might work for you.

Good Luck,

John

Use the ?seq command or a ?for loop.  I like seq() personally:

with(geom3d):
L:=[[-5, -5, 8], [-5, -1, 10], [-5, 3, 10], [-5, 7, 8], [-5, 8, -5], [-5, 8, 7], [-5, 10, -1], [-5, 10, 3], [-1, -5, 10], [-1, 7, 10], [-1, 10, -5],
[-1, 10, 7], [3, -5, 10], [3, 7, 10], [3, 10, -5], [3, 10, 7], [7, -5, 8], [7, -1, 10], [7, 3, 10], [7, 7, 8], [7, 8, -5], [7, 8, 7], [7, 10, -1],
[7, 10, 3], [8, -5, -5], [8, -5, 7], [8, 7, -5], [8, 7, 7], [10, -5, -1], [10, -5, 3], [10, -1, -5], [10, -1, 7],
[10, 3, -5], [10, 3, 7], [10, 7, -1], [10, 7, 3]]; seq(point(cat(`A_`,i),op(L[i])),i=1..nops(L)); eqS:=Equation(sphere(S,(x-1)^2 + (y-1)^2 +(z-1)^2 -121=0,[x,y,z],'centername'=T)); seq(TangentPlane(cat(`P_`,i),S,cat(`A_`,i)),i=1..nops(L)); seq(sort(Equation(cat(`P_`,i),[x,y,z])),i=1..nops(L));

 

John

Unfortutately, due to subtleties in how Maple represents symbolic products and quotients internally (see ?ProgrammingGuide,Chapter03 ), it is hard to refer to just part of a product using subs.  By using the whole product on the left hand side of the substitution, you can get the desired result:

f:=cos(sqrt(g)*sqrt(m+M)/(M^(3/2)*sqrt(l))*t);
subs(sqrt(g)*sqrt(m+M)/(M^(3/2)*sqrt(l))*t=omega*t,f); 

 

John

You can take a look at exactly what it is doing by using ?showstat

showstat(LinearAlgebra:-MinimalPolynomial);

It appears to be directly computing the linear dependencies on the powers of the matrix A by flattening A and its powers into rows of a larger n by n^2 matrix (systab in the code).

I actually wonder whether recent (ish) improvments in ?LinearAlgebra:-CharacteristicPolynomial would make it faster just call that and then ?PolynomialTools:-SquareFreePart

Adding to the suggestions in the comments, I would use 'add' instead of 'sum', 'mul' instead of 'product', use an explicit `*` for all multiplications (instead of implicit multiplications or `.`), and use tau_||i for variable names instead of tau(i). All of these help minimize problems. Doing that, I get the equations

eqn1 := (mul(1-tau_||i, i = 1 .. 4)*delta*Delta*tau_||0)*(1+(1/16)*(add((16-j+1)/(g_0)*(j), j = 1 .. 16))+(1-mul(1-tau_||i, i = 1 .. 4)*delta*Delta)*(1+(1/16)*(add((16-j+1)/(g_0)*(j), j = 1 .. 16)))+(1-mul(1-tau_||i, i = 1 .. 4)*delta*Delta)^2*(1+(1/32)*(add((32-j+1)/(g_0)*(j), j = 1 .. 32)))+(1-mul(1-tau_||i, i = 1 .. 4)*delta*Delta)^3*(1+(1/32)*(add((32-j+1)/(g_0)*(j), j = 1 .. 32)))+(1-mul(1-tau_||i, i = 1 .. 4)*delta*Delta)^4*(1+(1/32)*(add((32-j+1)/(g_0)*(j), j = 1 .. 32))))/(1-(1-mul(1-tau_||i, i = 1 .. 4)*delta*Delta)^3) = 1;

eqn2 := (mul(1-tau_||i, i = 2 .. 4)*(1-tau_||0)*delta*Delta*tau_||1)*(1+(1/8)*(add((8-j+1)/(g_1)*(j), j = 1 .. 8))+(1-mul(1-tau_||i, i = 2 .. 4)*(1-tau_||0)*delta*Delta)*(1+(1/8)*(add((8-j+1)/(g_1)*(j), j = 1 .. 8)))+(1-mul(1-tau_||i, i = 2 .. 4)*(1-tau_||0)*delta*Delta)^2*(1+(1/16)*(add((16-j+1)/(g_1)*(j), j = 1 .. 16)))+(1-mul(1-tau_||i, i = 2 .. 4)*(1-tau_||0)*delta*Delta)^3*(1+(1/16)*(add((16-j+1)/(g_1)*(j), j = 1 .. 16)))+(1-mul(1-tau_||i, i = 2 .. 4)*(1-tau_||0)*delta*Delta)^4*(1+(1/16)*(add((16-j+1)/(g_1)*(j), j = 1 .. 16))))/(1-(1-mul(1-tau_||i, i = 2 .. 4)*(1-tau_||0)*delta*Delta)^3) = 1;

eqn3 := ((1-tau_||0)*(1-tau_||1)*(1-tau_||3)*(1-tau_||4)*delta*Delta*tau_||2)*(1+(1/4)*(add((8-j+1)/(g_2)*(j), j = 1 .. 4))+(1-(1-tau_||0)*(1-tau_||1)*(1-tau_||3)*(1-tau_||4)*delta*Delta)*(1+(1/4)*(add((4-j+1)/(g_2)*(j), j = 1 .. 4)))+(1-(1-tau_||0)*(1-tau_||1)*(1-tau_||3)*(1-tau_||4)*delta*Delta)^2*(1+(1/8)*(add((8-j+1)/(g_2)*(j), j = 1 .. 8)))+(1-(1-tau_||0)*(1-tau_||1)*(1-tau_||3)*(1-tau_||4)*delta*Delta)^3*(1+(1/8)*(add((8-j+1)/(g_2)*(j), j = 1 .. 8)))+(1-(1-tau_||0)*(1-tau_||1)*(1-tau_||3)*(1-tau_||4)*delta*Delta)^4*(1+(1/8)*(add((8-j+1)/(g_2)*(j), j = 1 .. 8))))/(1-(1-f_2*delta*Delta)^3) = 1;

eqn4 := ((1-tau_||0)*(1-tau_||1)*(1-tau_||2)*(1-tau_||4)*delta*Delta*tau_||3)*(1+(1/2)*(add((2-j+1)/(g_3)*(j), j = 1 .. 2))+(1-(1-tau_||0)*(1-tau_||1)*(1-tau_||2)*(1-tau_||4)*delta*Delta)*(1+(1/2)*(add((2-j+1)/(g_3)*(j), j = 1 .. 2)))+(1-(1-tau_||0)*(1-tau_||1)*(1-tau_||2)*(1-tau_||4)*delta*Delta)^2*(1+(1/4)*(add((4-j+1)/(g_3)*(j), j = 1 .. 4)))+(1-(1-tau_||0)*(1-tau_||1)*(1-tau_||2)*(1-tau_||4)*delta*Delta)^3*(1+(1/4)*(add((4-j+1)/(g_3)*(j), j = 1 .. 4)))+(1-(1-tau_||0)*(1-tau_||1)*(1-tau_||2)*(1-tau_||4)*delta*Delta)^4*(1+(1/8)*(add((8-j+1)/(g_3)*(j), j = 1 .. 8))))/(1-(1-(1-tau_||0)*(1-tau_||1)*(1-tau_||2)*(1-tau_||4)*delta*Delta)^3) = 1;

eqn5 := (f_4*delta*Delta*tau_||4)*(1+(1*1)*(add((1-j+1)/(g_4)*(j), j = 1 .. 1))+(1-f_4*delta*Delta)*(1+(1*1)*(add((1-j+1)/(g_4)*(j), j = 1 .. 1)))+(1-f_4*delta*Delta)^2*(1+(1/2)*(add((2-j+1)/(g_4)*(j), j = 1 .. 2)))+(1-f_4*delta*Delta)^3*(1+(1/2)*(add((2-j+1)/(g_4)*(j), j = 1 .. 2)))+(1-f_4*delta*Delta)^4*(1+(1/4)*(add((4-j+1)/(g_4)*(j), j = 1 .. 4))))/(1-(1-f_4*delta*Delta)^3) = 1;

These evaluate to a polynomial system in the tau variables, which is good, since Maple is generally very good at solving such systems. However, in this case, since these equations have many parameters in addition to the variables, solving using any of the methods in Maple runs out of memory (6GB in my case) fairly quickly. Your best hope, is to specialize some of your parameters, and solve in each case.

John

I am guessing you are using solve(%,x); on this equation?  In that case, Maple will not be able to get an answer unless (unlike this case) the integral has a closed form.  There is no machinery inside solve to get at variables inside definite integrals.

John

Be careful definiting variables as x[i], since this implicitly creates a table named 'x' and something like this can happen to you:

(**) x[1]:=1;
                                   x[1] := 1

(**) x := 10;
                                    x := 10

(**) x[1];
                                     10[1]

 I prefer to use concatenation using 'cat' or the `||` operator for this sort of thing:

(**) x_||1 := 1;
                                   x_1 := 1

(**) for i from 1 to 10 do x_||i := i; end do;

 

John

In (1), ModuleApply appear to be passing two arguments, but in fact it is passing all of the arguments it received via the special variable _passed (aka args).

(**) foo := proc() print(_passed); end proc: # foo passes all its arguments to print
(**) foo(1,2,3,4);
                                  1, 2, 3, 4
(**) foo(blah);
                                     blah

 

John

In any case that you want to treat expressions as names, frontend is the one for the job:

(**) frontend(collect, [hh, v^x] );
#                                  2             x
#                              (5 x  - 4 x + 4) v

Two things:

1. Since there are a lot of parameters and the "explicit" solution contains a large non-numeric RootOf() the solution is likely not to be valid for many values of the parameters.   It is almost certainly a better idea, in this case, to specialize parameters and then solve.  Even if you need to do this a lot of times, it will likely be faster and more reliable than substituting in the generic solution.

2.  Use: 1/2, 1/10, and 10^6 instead of .5, .1 and 0.1e7.  Even though they are "mathematically" equivalent, the former are exact and will lead to exact computations, while the latter are floating point numbers that will cause round off errors in subsequent computations.  If you want round-offs then use ?fsolve instead of solve.

RandomTools is using the MersenneTwister PRNG by default. You can set its seed with the ?randomize command, or by calling its SetState command directly (see ?RandomTools,MersenneTwister,SetState ). This should be convincing that that works:

randomize(0); 
setp1 := {$1..100};

to 10 do
    decision := RandomTools:-Generate(choose(setp1)); # chooses 45 every time
    randomize(0);
end do; 

If you are just looking for reproducibility, you should not need to set a seed, since Maple will set itself to the same state every time you start or restart a session.

2 3 4 5 6 7 8 Page 4 of 10