Alec Mihailovs

Dr. Aleksandrs Mihailovs

4415 Reputation

20 Badges

17 years, 32 days
Mihailovs, Inc.
Owner, President, and CEO
Tyngsboro, Massachusetts, United States

Social Networks and Content at Maplesoft.com

Maple Application Center

I received my Ph.D. from the University of Pennsylvania in 1998 and I have been teaching since then at SUNY Oneonta for 1 year, at Shepherd University for 5 years, at Tennessee Tech for 2 years, at Lane College for 1 year, and this year I taught at the University of Massachusetts Lowell. My research interests include Representation Theory and Combinatorics.

MaplePrimes Activity


These are answers submitted by Alec Mihailovs

Assuming doesn't work like that, but if you just want to get rid of b, that could be done by assigning it to 0 in this example:

eval(eqn1,b=0);

                               y = a t

However, I wonder what result you would like to get, say for y=a*b, or for y=ab, or for y=sqrt(a^2+1)+b-a.

You might try some version of asympt, too (the error term would be wrong though, probably).

Alec

For such expressions with square roots, simplify shouldn't be used with option symbolic, just because of the following,

simplify(sqrt(x^2),symbolic);

                                  x

It should be used with assumptions instead,

simplify(sqrt(x^2)) assuming x<0;

                                  -x

which looks better than

simplify(sqrt(x^2));

                              csgn(x) x

Option size is useful sometimes, but sometimes - not, so you can try with it and without and see what is better. Sometimes a combination of them is useful - say simplify first and simplify with option size after.

Sometimes simplify is not useful at all - any combinations of it, with option size and without, make the original expression more complicated.

Alec

MultiSeries is a part of the algolib library. It doesn't seem to have Maple 9 version, but you might try Maple 8 version (or, maybe, Maple 10).

~ can be replaced with map for op and convert, so that op~(a) is replaced with map(op,a); and similarly with map2 for MultiSeries:-series and residues, moving the command to the first argument of map2.

I skipped some complicated looking output, but since you don't have a relatively recent Maple version, I am posting it here, in the input form,

poles:=
[u = -I, u = I, 
u = 2*beta*Q*Pi*(-beta+2*I*Pi)/(4*Pi^2+beta^2)-1/(4*Pi^2+beta^2)*
(4*beta^4*Q^2*Pi^2-16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1/2), 
u = 2*beta*Q*Pi*(-beta+2*I*Pi)/(4*Pi^2+beta^2)+1/(4*Pi^2+beta^2)*
(4*beta^4*Q^2*Pi^2-16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1/2), 
u = -2*beta*Q*Pi*(beta+2*I*Pi)/(4*Pi^2+beta^2)-1/(4*Pi^2+beta^2)*
(4*beta^4*Q^2*Pi^2+16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1/2), 
u = -2*beta*Q*Pi*(beta+2*I*Pi)/(4*Pi^2+beta^2)+1/(4*Pi^2+beta^2)*
(4*beta^4*Q^2*Pi^2+16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1/2)];

residues:=
[1/2*I/beta^2, -1/2*I/beta^2, 1/2*(-8*beta^2*Q*Pi^2-I*beta*
(4*beta^4*Q^2*Pi^2-16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2)-2*I*beta^3*Q*Pi+8*I*Q*Pi^3*beta-
2*Pi*(4*beta^4*Q^2*Pi^2-16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2))*Pi/(-4*beta^4*Q^2*Pi^2+16*I*beta^3*Q^2*Pi^3-2*beta^2*Q*Pi*
(4*beta^4*Q^2*Pi^2-16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2)+16*beta^2*Q^2*Pi^4+4*I*Pi^2*beta*Q*
(4*beta^4*Q^2*Pi^2-16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2)+16*Pi^4+8*beta^2*Pi^2+beta^4)/beta^2, 1/2*(-8*beta^2*Q*Pi^2+I*beta*
(4*beta^4*Q^2*Pi^2-16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2)-2*I*beta^3*Q*Pi+8*I*Q*Pi^3*beta+2*Pi*
(4*beta^4*Q^2*Pi^2-16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2))*Pi/(-4*beta^4*Q^2*Pi^2+16*I*beta^3*Q^2*Pi^3+2*beta^2*Q*Pi*
(4*beta^4*Q^2*Pi^2-16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2)+16*beta^2*Q^2*Pi^4-4*I*Pi^2*beta*Q*
(4*beta^4*Q^2*Pi^2-16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2)+16*Pi^4+8*beta^2*Pi^2+beta^4)/beta^2, -1/2*(-8*beta^2*Q*Pi^2+I*
(4*beta^4*Q^2*Pi^2+16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2)*beta+2*I*beta^3*Q*Pi-8*I*Q*Pi^3*beta-
2*Pi*(4*beta^4*Q^2*Pi^2+16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2))*Pi/(4*beta^4*Q^2*Pi^2+16*I*beta^3*Q^2*Pi^3+2*beta^2*Q*Pi*
(4*beta^4*Q^2*Pi^2+16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2)-16*beta^2*Q^2*Pi^4+4*I*Pi^2*beta*Q*
(4*beta^4*Q^2*Pi^2+16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2)-16*Pi^4-8*beta^2*Pi^2-beta^4)/beta^2, -1/2*(8*beta^2*Q*Pi^2+I*
(4*beta^4*Q^2*Pi^2+16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2)*beta-2*I*beta^3*Q*Pi+8*I*Q*Pi^3*beta-
2*Pi*(4*beta^4*Q^2*Pi^2+16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2))*Pi/(-4*beta^4*Q^2*Pi^2-16*I*beta^3*Q^2*Pi^3+2*beta^2*Q*Pi*
(4*beta^4*Q^2*Pi^2+16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2)+16*beta^2*Q^2*Pi^4+4*I*Pi^2*beta*Q*
(4*beta^4*Q^2*Pi^2+16*I*beta^3*Q^2*Pi^3-16*beta^2*Q^2*Pi^4-16*Pi^4-8*beta^2*Pi^2-beta^4)^(1
/2)+16*Pi^4+8*beta^2*Pi^2+beta^4)/beta^2];

I plotted the imaginary parts of the 3rd, 4th, 5th, and 6th poles for various fixed values of Q, as a function of beta, and it looked as if (at least for those values that I plotted) the sign of imaginary part was depending only on beta and not on Q, so if beta>0 the poles that I mentioned had positive imaginary part, both for positive Q and for negative Q.

For both beta and Q positive, the second pole is I, so it has a positive imaginary part, the 3rd pole has the positive imaginary part (because the sign of the imaginary part of a square root of a non-real number is the same as the sign of the imaginary part of that number). The 5th has a negative imaginary part.

The 4th and the 6th are not that obvious. But they are conjugate, so one of them has a positive imaginary part, and another one - a negative imaginary part. Now, they are continuous functions of Q and beta for positive Q and beta, so their imaginary part is also continuous. That means that, if, say the 6th pole has imaginary parts of different signs for different values of Q and beta, there should be such Q and beta for which the imaginary part is 0. But looking at the function s2, it is obvious that it can't have real poles, so such a case is impossible, and the poles have the same sign on the entire domain with Q>0 and beta>0. So check that sign in one point (any) in that domain, and it should be the same everywhere there.

So the 2nd, 3rd, and 6th pole have positive imaginary parts for positive beta and Q, and other poles have negative imaginary part.

Alec

Seems to be working for me,

pd := VectorCalculus:-Laplacian(u(r,phi),polar[r,phi])=A*cos(phi):

simplify(pdsolve(pd));

  u(r, phi) =

                        2
        1/3 A cos(phi) r  + _F1(ln(r) I + phi) + _F2(-ln(r) I + phi)

simplify(eval(pd,%));

                       A cos(phi) = A cos(phi)

_F1 and _F2 arguments can be changed to their usual form through

convert(exp((-ln(r)*I+phi)*I),trig);

                      r (cos(phi) + sin(phi) I)

convert(exp((ln(r)*I+phi)*(-I)),trig);

                      r (cos(phi) - sin(phi) I)

Alec

Also, that could be done without a loop as

f:=Vector(2,i->(x,y)->x+y^i):

f[1](x,y);

                                x + y

f[2](x,y);

                                     2
                                x + y

That works also with the proc notation,

f:=Vector(2,i->proc(x,y) x+y^i end):

Instead of unapply, also the following can be used in a loop,

for j from 1 to 2 do
f[j]:=parse(sprintf("proc(x,y)
x+y^%d
end",j))
od;
                  f[1] := proc(x, y) x + y end proc
                 f[2] := proc(x, y) x + y^2 end proc

Another way is

for j from 1 to 2 do
f[j]:=(j->proc(x,y)
x+y^j
end proc)(j)
od;
                 f[1] := proc(x, y) x + y^1 end proc
                 f[2] := proc(x, y) x + y^2 end proc

which could be also written as

F:=j->proc(x,y) x+y^j end proc:
for j from 1 to 2 do f[j]:=F(j) od;

                 f[1] := proc(x, y) x + y^1 end proc
                 f[2] := proc(x, y) x + y^2 end proc

There are many other ways.

Alec

You could try ?RootFinding:-Analytic .

fsolve with option avoid is trying to find other roots than those specified in avoid, see ?fsolve/details . This can be put into a loop, with adding every new found root to avoid.

Alec

For example,

with(Ore_algebra):
A:=diff_algebra([Dx,x],[Dy,y]):
`&*`:=rcurry(skew_product,A):
d:=Dx-I*Dy:
p:=2*y^2-3*y+1:
d&*(p&*d);

                                         2              2
  (-4 I y + 3 I) Dx + (3 - 4 y) Dy + (2 y  - 3 y + 1) Dx

                  2                             2              2
         + (-4 I y  + 6 I y - 2 I) Dx Dy + (-2 y  + 3 y - 1) Dy

Alec

If beta≠0 and Q≠0, the poles and residues can be found as

_EnvExplicit:=true:
s2 := (1/2/Pi+2*Q*u/(1+u^2))^2/((1+beta^2*(1/2/Pi+2*Q*u/(1+u^2))^2)*(1+u^2)):
poles:=op~([singular(s2,u)]);
residues:=zip((x,y)->simplify(x*(lhs-rhs)(y)),
    convert~(MultiSeries:-series~(s2,poles,1),polynom),
    poles); 

If beta=0, then s2 is

s20:=simplify(eval(s2,beta=0));

                                  2            2
                            (1 + u  + 4 Q u Pi)
                     s20 := --------------------
                                  2       2 3
                              4 Pi  (1 + u )

so the poles are ±I, and the residues are

residues0:=simplify(residue~(s20,[u=-I,u=I]));

                               2  2                   2  2
                    1/8 I (2 Pi  Q  + 1)  -1/8 I (2 Pi  Q  + 1)
      residues0 := [--------------------, ---------------------]
                              2                      2
                            Pi                     Pi

Now, if there are no real poles, the integral is 2*Pi*I*(sum of residues at the poles with positive imaginary part).

So, if beta=0, the residue formula gives

2*Pi*I*residues0[2];

                                 2  2
                             2 Pi  Q  + 1
                             ------------
                                 4 Pi

and that coincides with the value calculated by Maple directly,

int(s20,u=-infinity..infinity);

                                 2  2
                             2 Pi  Q  + 1
                             ------------
                                 4 Pi

The next case is Q=0. In this case s2 is

s2Q:=simplify(eval(s2,Q=0));

                                     1
                   s2Q := ------------------------
                               2       2        2
                          (4 Pi  + beta ) (1 + u )

so the poles are again ±I, and the residues are

residuesQ:=simplify(residue~(s2Q,[u=-I,u=I]));

                               1/2 I         -1/2 I
             residuesQ := [-------------, -------------]
                               2       2      2       2
                           4 Pi  + beta   4 Pi  + beta

and the residue formula gives

2*Pi*I*residuesQ[2];

                                 Pi
                            -------------
                                2       2
                            4 Pi  + beta

which equals the integral calculated directly,

int(s2Q,u=-infinity..infinity);

                                 Pi
                            -------------
                                2       2
                            4 Pi  + beta

Now, if beta≠0 and Q≠0, it looks from the plots that assuming that beta and Q are real, if beta>0, the 2nd, 3rd, and 6th poles have positive imaginary part, while if beta<0, the 2nd, 4th, and 5th poles have positive imaginary part.

Assuming that this is true, the residues formula gives for beta>0 value

int1:=2*Pi*I*(residues[2]+residues[3]+residues[6]);

and for beta<0 value

int2:=2*Pi*I*(residues[2]+residues[4]+residues[5]);

We can check that numerically,

f:=(b,q)->evalf(eval(`if`(b>0,int1,int2),[beta=b,Q=q])):
g:=(b,q)->evalf(Int(eval(s2,[beta=b,Q=q]),u=-infinity..infinity)):

f(1,2)=g(1,2);
                                        -9
           1.735213452 - 0.1633628180 10   I = 1.735213454

f(-2,1)=g(-2,1);
                                        -9
          0.4327637188 + 0.3644247478 10   I = 0.4327637189

f(-1,-1)=g(-1,-1);

                  0.9273289786 + 0. I = 0.9273289789

f(2,2)=g(2,2);
                                        -9
          0.5944425012 - 0.5692565890 10   I = 0.5944424990

Alec

That was discussed in a few old threads. For example, here and here,

The cliff notes are that it is possible in Maple, but it is not as convenient as it is in Magma or Sage.

For example,

alias(a=RootOf(_Z^4+_Z+1)):
f:=a^3*x+(1+a^2)*y:
g:=a*x+y:

Expand(f*g) mod 2;

                     2            2  2    2      2
                    y  + x y a + y  a  + x  a + x

is not collected, so one has to use something like

h:=collect(Expand(f*g) mod 2,{x,y});

                              2                 2   2
                h := (a + 1) x  + x y a + (1 + a ) y

which has some terms (xya in this example) strangely ordered.

For division, one can use Normal,

collect(Normal(h/g) mod 2,{x,y});

                           3           2
                          a  x + (1 + a ) y

Factoring works as

Factor(h) mod 2;

                     2                   2
               (1 + a ) (a x + y) (y + (a  + a + 1) x)

Alec

Student:-Calculus1:-Asymptotes~([allvalues(RootOf(8*x^3+y^3-6*x*y-3,y))],x,y)[1];

                            [y = -2 x - 1]

Student:-Calculus1:-Asymptotes~([allvalues(RootOf(16*x^4-y^4-y^3-3*x^2+y+1,x))],y,x)[1];

                   [x = y/2 + 1/8, x = - y/2 - 1/8]

[seq](solve({i},y)[],i=%);

                   [y = 2 x - 1/4, y = -2 x - 1/4]

Alec

If the speed is not a concern, instead of adjusting Digits before using frem, the following procedure can be used for determining angles,

angle:=evalf@argument@exp@curry(`*`,I):

For example,

angle(10^23);

                             2.364596689

sin(%);

                             0.7011406397

sin(10.^23);

                             0.7011406399

It produces angles in the interval between -Pi and Pi, the same as frem,

angle(75);

                             -0.39822370

Alec

How can we tell? That depends on your calculations.

If you run an infinite loop, it may continue for a long time. If your calculations will end in an hour or two, they will.

Generally speaking, there is no a memory limit set for Maple usually (except, maybe, say 1GB for Java - but that's different, I think; or if you set such a limit starting Maple) - it will try to use all of computer's memory available if possible.

It is often a good idea to try your big calculations on a smaller example(s) first - with short loops and small amount of data. If they work - increase the loop length and the data size slightly, and compare the time and memory used - then you might be able to estimate the time and memory needed for your big calculations.

Also, you could try running it in the command line Maple - that frees resources used for the worksheet interface.

Typically, such a thing as you described, happens if you are outputting a large amount of data in the worksheet - even few thousand lines could hang Maple (especially Standard interface and especially in older Maple versions.)

Alec

It works without functional definitions as

alias(zeta=zeta(t),Yt=diff(y(t,x),t)):
A :=int(Yt^2, x = 0 .. zeta):
lprint(diff(A,t));

int(2*Yt*diff(Yt,t),x = 0 .. zeta)+diff(zeta,t)*D[1](y)(t,zeta)^2

Alec

@Alejandro Jakubi 

Or dividing the integration interval in 2 parts - from 0 to 1 and from 1 to infinity, with the first integral being

int(exp(-x)/(1-exp(-x)),x=0..1);

                               infinity

Alec

A workaround (not perfect) is

_EnvContinuous:=true:
Re(Mean(R));

                      1/2                 1/2 2
                     2              (1 + 2   )
                     ---- + 1/12 ln(-----------)
                      3               1/2     2
                                    (2    - 1)

expand(%);
               1/2
              2                  1/2            1/2
              ---- + 1/6 ln(1 + 2   ) - 1/6 ln(2    - 1)
               3

The problem seems to be in evaluating that integral - a workaround for that is dividing the integration interval in 2 parts - from 0 to 1 and from 1 to sqrt(2).

_EnvContinuous:=false:
int(rho*pdf, rho = 0..1) + int(rho*pdf, rho = 1..sqrt(2));

               1/2
              2                  1/2            1/2
              ---- + 1/6 ln(1 + 2   ) - 1/6 ln(2    - 1)
               3

eval(%,sqrt(2)-1=1/(sqrt(2)+1));

                        1/2
                       2                  1/2
                       ---- + 1/3 ln(1 + 2   )
                        3

evalf(%);

                             0.7651957162

Alec

2 3 4 5 6 7 8 Last Page 4 of 76