Alec Mihailovs

Dr. Aleksandrs Mihailovs

4485 Reputation

21 Badges

20 years, 99 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

I would just use Mayavi.

Alec

Did you try evalhf instead of evalf, or using evalf with 14 digits, i.e. evalf(..., 14) ?

Alec

I think, I understand that now. It says - "Left 10 comments", so more than 10 (as well as less than 10) don't qualify.

Alec

assign((eq->lhs(eq)=unapply(rhs(eq),indets(rhs(eq))[]))~(resultset));

Alec

Substitution z=b*w reduces that to the case with b=1. Now, taking the contour close to the real axis, (x+I*y), sqrt((x^2+I*y)^2+1), and sqrt(1-1/(x+I*y)^2) with small y have the same sign of the imaginary part as the sign of y, and the first 2 of these imaginary parts are small, so that product also will have the same sign of the imaginary part as the sign of y.

The sign of the imaginary part of sqrt(z^4-1) is positive for real z on the interval between -1 and 1. That means that the integral over the bottom part of the contour (close to the real axis), oriented counterclockwise, will be close to the -int(sqrt(x^4-1),x=-1..1), and the same for the integral over the top part of that contour. Taking the limit, our contour integral (with b=1) equals

-2*int(sqrt(x^4-1),x=-1..1)=-4*int(sqrt(x^4-1),x=0..1);

                                  1/2
                   1/2           2
           -4/3 I 2    EllipticK(----) = -I Beta(1/4, 3/2)
                                  2

which we can also express in terms of Γ,

convert(rhs(%),GAMMA);

                                   3/2  1/2
                          -2/3 I Pi    2
                          -----------------
                                       2
                             GAMMA(3/4)

And in the general case the answer will be that expression multiplied by b^3 = E^(3/4).

If the originally chosen function had the opposite sign, i.e. if it was -z*sqrt(b^2+z^2)*sqrt(1-z^2/b^2), the sign of the (imaginary part of the) contour integral would be opposite as well.

Alec

PS By the way, this expression of EllipticK(1/sqrt(2)) through GAMMA(3/4) seems to be unknown to simplify,

simplify(2*EllipticK(sqrt(2)/2)*GAMMA(3/4)^2/Pi^(3/2));

                                 1/2
                                2               2
                    2 EllipticK(----) GAMMA(3/4)
                                 2
                    -----------------------------
                                  3/2
                                Pi

evalf(%);

                             0.9999999994

The conversion to hypergeom works though,

simplify(convert(EllipticK(1/sqrt(2)),hypergeom));

                                    3/2
                                  Pi
                           1/2 -----------
                                         2
                               GAMMA(3/4)

but simplify doesn't use this conversion itself, even if the option hypergeom is added to it.

I'm not actually suggesting to add that to simplify - but, perhaps, a new FullSimplify command could be added, working longer, perhaps, than simplify, but trying various combinations of possible simplifications.

Alec

In Standard Maple on Windows,

system("taskkill /im maple* /f /t");

I guess, that can be done similarly on other systems. It's not hard to find the code on the web - I just can't test it (at least for Mac - never had these beasts at home).

Alec

simplify(t^(1/n),symbolic)^n;

                                     n
                                (2/5)

Alec

MultiSeries:-limit(((x^3-2*x^2+sin(x))^(1/3))/(6*x+2),x=infinity);

                                 1/6

Alec

Copy OrthogonalExpansions.mla and OrthogonalExpansions.hdb from your Maple 14/lib folder and paste to Maple 15/lib, restart Maple 15 and it will work.

Alec

While Maple is not able to do the 2-dimensional case directly, it has no problems with evaluating twice a half of the integral,

8*int(int(sqrt(x^2+y^2)*(1-x)*(1-y),y=0..x),x=0..1);

                    1/2              1/2        1/2         1/2
           58 + 41 2    + 85 ln(1 + 2   ) + 60 2    ln(1 + 2   )
      1/15 -----------------------------------------------------
                                        1/2
                               17 + 12 2

radnormal(%);

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

Alec

That could be done without rsolve as

f:=proc(k) option remember; <1,1;-1,1>.procname(k-1) end:
f(0):=<0,1>:

plot(Matrix(21,2,(i,j)->f(i-1)[j]));

Or, using rsolve, as

rec:=rsolve({x(k+1)=x(k)+y(k),y(k+1)=-x(k)+y(k),x(0)=0,y(0)=1},{x(k),y(k)}):

plot([seq(eval([x,y](i),eval(rec,k=i)),i=0..20)]);

Alec

Numerical approximation (up to 4 digits after the decimal dot) can be done in Maple as

f:=(x1,y1,x2,y2,x3,y3)->sqrt((x1-y1)^2+(x2-y2)^2+(x3-y3)^2):
evalf(Int(f,[(0..1)$6],method= _MonteCarlo,epsilon=0.5e-4));

                             0.6616946352

One-dimensional case can be done symbolically,

int(abs(x-y),[x=0..1,y=0..1]);

                                 1/3

Two-dimensional case can be done with 14 digits,

f:=unapply(simplify(int(sqrt((x-y)^2+c),[x=0..1,y=0..1])) ,c) assuming c>0;

                 (3/2)                    1/2
  f := c -> 2/3 c      - c ln(-1 + (c + 1)   ) + 1/2 c ln(c)

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

Digits:=14:
evalf(Int(f((x-y)^2),[x=0..1,y=0..1]));

                           0.52140543316473

or with higher precision after a slight modification of the integral,

Digits:=100:
4*eval(int(int(sqrt(x^2+y^2)*(1-x)*(1-y),x=0..1),y=0..1),csgn(y)=1);

       1
      /      4      3
     |      y      y          3     2         2     2
  4  |   - ---- + ---- + 1/4 y  ln(y ) - 1/4 y  ln(y )
     |      3      3
    /
      0

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

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

evalf(%);

  0.52140543316472067833098235660724397491403156777900834179621051\
        87505078933048158318679281329252614524

Maple has no problems with integrating symbolically the part of the integrand free of logarithms. It has problems with the following integral,

4*int(1/4*y^3*ln(y^2)-1/4*y^2*ln(y^2)-1/2*y^3*ln(1+(y^2+1)^(1/2))+1/2*y^2*ln(1+(y^2+1)^(1/2)),y = 0 .. 1)=1/2*2^(1/2)-1/3+1/3*ln(1+2^(1/2))+1/2*ln(2^(1/2)-1);

       1
      /
     |        3     2         2     2         3          2     1/2
  4  |   1/4 y  ln(y ) - 1/4 y  ln(y ) - 1/2 y  ln(1 + (y  + 1)   )
     |
    /
      0

                2          2     1/2
         + 1/2 y  ln(1 + (y  + 1)   ) dy =

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

evalf(%);

                     0.2268778500 = 0.2268778498

Rewriting similarly the integral for the 3-dimensional case, we can evaluate it numerically more precisely,

8*int(int(int(sqrt(x^2+y^2+z^2)*(1-x)*(1-y)*(1-z),x=0..1),y=0..1),z=0..1):
evalf(%);

                             0.6617071822

Similarly, it works in higher dimensions up to n=7,

f:=n->evalf(2^n*Int(sqrt(add((x||i)^2,i=1..n))*mul(1-x||i,i=1..n),
    [seq(x||i=0..1,i=1..n)]));
seq(f(n),n=1..7);

  0.3333333334, 0.5214054332, 0.6617071822, 0.7776656536,
                                                      
        0.8785309152, 0.9689420832, 1.051583873

Alec

evalc might work without additional rules, as in the original example.

Alec

evalc(r);

                   2      2
                  k  sigma                       1/2   1/2
            exp(- ---------) cos(k theta) sigma 2    Pi
                      2

Alec

It would be better to have more complete example of what you would like to get - an expression and the wanted result.

But in general, the function name can be obtained using op(0,f),

f:=y(n):
op(0,f),op(f);
                                 y, n

Alec

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