Axel Vogt

5018 Reputation

19 Badges

15 years, 315 days
Munich, Germany

MaplePrimes Activity


These are answers submitted by Axel Vogt

You will lose that additional digits if you continue to operate with the result, but you can try something like this (without going into details for different length)

Digits:=10;
a:=.123456789;
m:=proc(a,b) evalf[2*Digits](a*b) end proc;
m(a,a);
                         0.015241578750190521
m(a,a) + 1;
                             1.015241579

I only have some test version, but there is a (general) workaround:

Go to the very directory, ..\bin.win,
mark the program file, maplew.exe, using the mouse 
and now use the context menu of the mouse (= right click) to add the program to the start menu.
Dito for the classical interface cwmaple.exe

I used a change x = 1/3000*46^(1/2)*z * t, t) assuming 0 < z and then
z = w / ( 1/3000*46^(1/2) ) to get
w^3*Int(t^2*ln(1+exp(-w*(t^2-1)^(1/2))),t = 0 .. infinity)
Now real / imaginary is decided in t=1 and using eval(Re(..)), 0 < z I arrive at
Int(1/2*t^2*ln((1+cos(w*(-t^2+1)^(1/2)))^2+sin(w*(-t^2+1)^(1/2))^2), t=0..1, 
  digits = 10, method = _d01ajc) +
Int(t^2*ln(1+exp(-w*(t^2-1)^(1/2))), t=1 .. infinity,
  digits = 10, method = _d01amc);
w^3*%;
plot(%, w = 0 .. 30, numpoints=5);
That shows the rough locations of several extrema of a oscillating and
increasing function.
It does not directly answer your question, but for the above reformulation
you can diff for w formally and integrate and fsolve more "easily".
Hope that idea helps.
Here is mine, actually also handmade, but written in Maple (and perhaps
just some reformulation of the other answers):

  theAssumptions:= 0 < Omega, 0 < a, 0 < k, 0 < m;

Following Preben one can set Omega = 1 = m,

  subs(Omega=1, m=1, ex);
  P:=simplify(%) assuming theAssumptions;

    -2^(1/2)*(-a^2-2*k+a*(a^2+4*k)^(1/2))^(1/2) / ((a^2+4*k)^(1/2)-a)       

Maple has problems to find signs for the inputs to sqrt, but help it (and
have in mind strict convexity of sqrt):

  -a^2-2*k+a*(a^2+4*k)^(1/2):
  0 < (-1)*%;
  % +a*(a^2+4*k)^(1/2);
  lhs(%)^2 < rhs(%)^2; expand(%);
  is(%) assuming theAssumptions;

                         2             2       1/2
                    0 < a  + 2 k - a (a  + 4 k)

                                 true

Dito for the denominator. Now using sqrt(-p) = I*sqrt(p) for 0 < p we have

  P = -I *  2^(1/2)*(+a^2+2*k-a*(a^2+4*k)^(1/2))^(1/2)/((a^2+4*k)^(1/2)-a)

Maple still does not simplify the very term to 1, but knowing 'positive'
it now is enough to have it for the squares, like vv did.

  2^(1/2)*(+a^2+2*k-a*(a^2+4*k)^(1/2))^(1/2)/((a^2+4*k)^(1/2)-a);
  %^2; #map(expand, %);
  expand(numer(%)) / expand(denom(%)); # else denom is not expanded

                                  1

BTW: thx for the neat example

Often one way is to reduce to constant intervals.

In many cases order of integration does not matter (but you may check for your specific case), so interchange them and consider Int(Int(f(a,b),b=1-a..1+a),a=0..10). You may use Upper Case instead of lower case for integration untill the end.  

The inner integral can be mapped to any interval, in the following I map it to -1 ... 1. Then your problem has vanished.

  Int(f(a,b),b = 1-a .. 1+a);
  IntegrationTools:-Change(%, b = a*t+1, t): combine(%);
  Int(%,a=0..10);
                      10    1
                     /     /
                    |     |
                    |     |   f(a, a t + 1) a dt da
                    |     |
                   /     /
                     0     -1
Example:
  Int(Int(f(a,a*t+1)*a,t = -1 .. 1),a = 0 .. 10);
  eval(%, f = '(a,b) -> exp(-a^3-b)');
  value(%): subs(int=Int, %); # Maple does not know a symbolic solution
  evalf(%);
                   10    1
                  /     /
                 |     |         3
                 |     |   exp(-a  - a t - 1) a dt da
                 |     |
                /     /
                  0     -1

                10
               /
              |          3                  3
              |    exp(-a  + a - 1) - exp(-a  - a - 1) da
              |
             /
               0

                          0.370721316428104

Edited to add:

For that example good luck reduced the task even to dim = 1.

You can rewrite that final integrand as exp(-a^3)*exp(-1)*sinh(a)*2
to reduce (additive) cancellation errors if you want high precision
(though Maple should care for that on its own).

It is like integrating in a pole:

Int((1/2)/(s^(1/2)*GAMMA(1/2)*(t^(1/2)-s^(1/2)))*(s^6), s = 0 .. t);
IntegrationTools:-Change(%, s=x*t, x);
IntegrationTools:-GetIntegrand(%);
expand(%);
series(%, x=1, 3);
       2/(x-1)+23/2+219/8*(x-1)+O((x-1)^2)

So it is like int(1/(x-1), x=0..1), since the numerator in x=1 is finite.

You want to expand it (the ordering may change)

  x= (a+b*c)/d: 
  expand(%);
                                b c
                            x = --- + a/d
                                 d
It seems that 'simplify' writes (u-1)*(u+1) as u^2 - 1 and then does not factor again
restart: interface(version); # mine is Maple 18
c:= (1+sqrt(5))/2: #c:=2;
theAssumption:= 1<u,u<c:
g:=((u-1)*(u+1))^(1/2)*u/(u*(u-1))^(1/2);
simplify(g);
                              2     1/2
                            (u  - 1)    u
                            --------------
                                       1/2
                            (u (u - 1))
combine(g) assuming theAssumption;
simplify(%) assuming theAssumption;
#combine(%) assuming theAssumption;
                                  1/2  1/2
                           (u + 1)    u

In your case you can scale the matrix by it's determinant, solve and scale back. Find a file attached

MP_EigenValues_scaling.mw

MP_EigenValues_scaling.pdf

numer(%)/combine(denom(%)) and use Pi*x=arccos(r), Pi*y=arccos(s),
use cos(n*arccos(r)) = orthopoly[T](n,r) to come up with

Int(Int(
  2/(-s^2+1)^(1/2)/Pi^2*(T(n,r)*T(n,s)-1)/(-2+r+s)/(-r^2+1)^(1/2),
  r = -1 .. 1),s = -1 .. 1);

         1    1
        /    /
       |    |            2 (T(n, r) T(n, s) - 1)
       |    |   ------------------------------------------ dr ds
       |    |      2     1/2   2                 2     1/2
      /    /    (-s  + 1)    Pi  (-2 + r + s) (-r  + 1)
        -1   -1

Do it stepwise to get results for numerical values of n.
Here is an idea how to achieve it over the integers (I think there was some
similar question lately)

Generate an strictly upper triangle matrix with integer entries.

Generate an n-tuple from {+1,-1}. Take it, if it has an odd even number of "-1"s
and fill the diagonal.

The determinant is +1. Replace exactly one arbitrary +-1 by +-2. Then det = 2.

Making it more arbitrary now generate elementary transformations invertible
over Z (like changing rows or columns [or transpositions]) and apply it.

I am not the one to talk about how "random" the result will be.

That is what I see opening you sheet. Very funny.

 

Int(Int(f(x)*g(y),x),y); expand(%);

Int(Int(f(x)*g(y),x=-infinity..infinity),y=-infinity..infinity); expand(%);

The following is from the cited book: read cos as the real part of exp and
find a path which reduces oscillation (I always admire such).

f:= x -> exp(ln(x)/x*I)/x;

'eval(f(z(t))*D(z)(t), z = 't -> t+I*t*(1-t)')';
simplify(%): # do not simplify it too much
g:=unapply(%, t);

Int(Re@g, 0 .. 1, method = _Dexp); evalf[40](%); # quickly found

              0.3233674316777787613993700879521704466510

That also works for crazy stuff like cos(ln(x^1)/x^2)/x^3, which may
serve as another test for MMA.
By simplification one can ignore k1, k2, but has R to be arbitrarly real.

a := cos(x)*exp((R*sin(x)-4*x)*I)*sin(x);

Plotting it for some data "shows" that the imaginary part is symmetric,

a = eval(a, x=-x); map(Im, %); simplify(%) assuming R::real, x::real;
is(%);

                                 true
The real part integrates to 0:

evalc(Re(a)); combine(%);
int(%,x=-Pi/2 .. Pi/2) assuming R::real;

                  cos(x) cos(R sin(x) - 4 x) sin(x)

         -1/4 sin(-6 x + R sin(x)) + 1/4 sin(-2 x + R sin(x))

                                  0

And the imaginary part can be found:

evalc(Im(a)); #combine(%);
expand(%); #simplify(%);
2*Int(%,x=0 .. Pi/2) * I;

value(%) assuming R::real; #expand(%);
simplify(%, size);

It looks similar to MMA's automatic result, see at the end

For R = 2.3 I get -.693587157054838*I and for the integral

Int(a, x=-Pi/2..Pi/2);
subs(R=2.3, %); evalf(%);

                       0. - 0.693587157054844 I

with k1=1, k2 = 0 in the original formulation.

 

(4*R*Pi*(R^4-48*R^2+240)*BesselJ(1,R)+(36*Pi*R^4-480*Pi*R^2)*BesselJ(0,R)+(-2*R^5+224*R^3-1920*R)*cos(R)+(34*R^4-864*R^2+1920)*sin(R))/R^6

5 6 7 8 9 10 11 Last Page 7 of 86