## Alec Mihailovs

Dr. Aleksandrs Mihailovs

## 4470 Reputation

20 years, 14 days
Mihailovs, Inc.
Owner, President, and CEO
Tyngsboro, Massachusetts, United States

## Social Networks and Content at Maplesoft.com

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.

## Mayavi...

I would just use Mayavi.

Alec

## evalhf or evalf with 14 digits...

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

## In one line...

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

Alec

## Completing the calculation...

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

## brutal...

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

n
(2/5)
```

Alec

## MultiSeries:-limit...

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

1/6
```

Alec

## Copy and paste...

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

## Two-dimensional case...

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

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

Alec

## Another way...

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...

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...

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

Alec

## evalc...

```evalc(r);

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

Alec

## op(0,f)...

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
﻿