Axel Vogt

5821 Reputation

20 Badges

20 years, 222 days
Munich, Bavaria, Germany

MaplePrimes Activity


These are answers submitted by Axel Vogt

This is "by design" - consider the case that the integral can not be found.
You can use the following:
 
f:= x -> x;
F:= x -> int( f(X), X= 0 .. x);
F(x);
                                   2
                                  x
                                 ----
                                  2
# may be you want the following:
int( f(X), X= 0 .. x);
F1:=unapply(%, x);
                                   2
                                  x
                                 ----
                                  2

                                          2
                          F1 := x -> 1/2 x

"fsolve" gives a numerical result, depending on "Digits" - and it is a good idea to use Digits:=15 as default.

"solve" gives a symbolic result - but sometimes not all of them

You may use

solve( exp(z) = -1, allsolutions = true);

                          Pi I + 2 I Pi _Z1

Here the strange notation _Z1 is Maple's standard way to say that this is an arbitrary integer.

MP_compute_powers_loop.mws

# http://www.mapleprimes.com/questions/211577-Why-The-Evaluation-Of-Root-Of-A-Real

restart; interface(version); Digits:=15;

`Classic Worksheet Interface, Maple 2016.1, Windows, Apr 19 2016, Build ID 1132667`

Digits := 15

(1)

tst0:=proc(d::posint)
local i,j;
for i from 1 to d do
  for j from 1 to d do
    evalf(abs(i-j+1)^0.3-abs(i-j)^0.3):
  end do:
end do:
return 0;
end proc;

tst0 := proc (d::posint) local i, j; for i to d do for j to d do evalf(abs(i-j+1)^.3-abs(i-j)^.3) end do end do; return 0 end proc

(2)

tst1:=proc(d::posint)
local i,j;
for i from 1 to d do
  for j from 1 to d do
    evalhf(abs(i-j+1)^0.3-abs(i-j)^0.3):
  end do:
end do:
return 0;
end proc;

tst1 := proc (d::posint) local i, j; for i to d do for j to d do evalhf(abs(i-j+1)^.3-abs(i-j)^.3) end do end do; return 0 end proc

(3)

dim:=1000;

dim := 1000

(4)

#forget(tst0): gc():
CodeTools[Usage](evalhf(tst0(dim)));

memory used=0.59KiB, alloc change=0 bytes, cpu time=827.00ms, real time=826.00ms, gc time=0ns

 

0.

(5)

CodeTools[Usage](tst0(dim));

memory used=1.44GiB, alloc change=36.00MiB, cpu time=28.69s, real time=28.69s, gc time=795.61ms

 

0

(6)

CodeTools[Usage](tst1(dim));

memory used=232.60MiB, alloc change=0 bytes, cpu time=3.06s, real time=3.06s, gc time=109.20ms

 

0

(7)

forget(tst1): gc():
CodeTools[Usage](evalhf(tst1(dim)));

memory used=1.16KiB, alloc change=0 bytes, cpu time=764.00ms, real time=765.00ms, gc time=0ns

 

0.

(8)

 

cp_tst0:=Compiler:-Compile(tst0):

 

 

CodeTools[Usage](cp_tst0(dim));

memory used=1.21KiB, alloc change=0 bytes, cpu time=125.00ms, real time=125.00ms, gc time=0ns

 

0

(9)

 

 


Download MP_compute_powers_loop.mws

Maple has a "memory" for results, so times are not constant.

For computing it depends how you instruct Maple, for usual double precision that 1 million operations can be done from ~ 30 sec donw to 1/8 sec, see file attached

 

A symbolic solution is unlikely to exist, so you need to provide values
for your parameters to find the zeros numerically.

Here is a recipe how to do it for the real axis.

restart; Digits:=15;
 -(((beta*eta^2-(1/2)*beta+1)*p^2-(1/2)*beta^3)*KummerM((1/4)*((-beta+2)*
   p^2-beta^3)/p^2, 1, beta*eta^2)+     (1/2)*KummerM((1/4)*((-beta+6)*
   p^2-beta^3)/p^2, 1, beta*eta^2)*((beta-2)*p^2+beta^3))*
   exp(-(1/2)*beta*eta^2)/(p^2*beta);
 
 eval(%, [p = 1, eta = 1]); # some numerical values for parameters;
 
 f:=unapply(%, beta);
   beta -> 
     -((1/2*beta+1-1/2*beta^3)*KummerM(-1/4*beta^3-1/4*beta+1/2,1,beta)+
     1/2*KummerM(-1/4*beta^3-1/4*beta+3/2,1,beta)*(beta^3+beta-2))
     *exp(-1/2*beta)/beta
 # plot(f, -10 .. 10);
Maple does not see that there is a continuous extension to 0 being a zero.
 with(RootFinding);
  L:=NULL:
  t:=1e-3:               # starting point beyond beta=0
  for j from 1 to 3 do   # find 3 zeros, running to + infinity
    t:=NextZero(f, t);   # using that command, see the (crude) help
    L:= L, t;
  end do: t:='t':
  L:=[L];                # list of zeros
     L := [1.87588952396938, 2.58694685291078, 3.13796881735298]
 # check
 map(f, L);
                            -14                       -13
       [0.625993411295519 10   , -0.159057950051232 10   , -0.]
 # the negative values are zeros as well, a symmetry
 negL:=map(q -> -q, L);
  negL := [-1.87588952396938, -2.58694685291078, -3.13796881735298]
 # check that
 map(f, negL);
                        -14                      -13                       -13
  [-0.953311904713793 10   , 0.169099417324297 10   ,  0.336647449173820 10   ]
 # otherwise use g:=unapply(f(-x), x) and find zeros for g
 
Avoiding lists I would write that as

> convert(532, decimal, octal) -
> convert(235, decimal, octal);
> convert(%, octal);
                                 189

                                 275
 
> convert(275, decimal, octal) +
> convert(575, decimal, octal);
> convert(%, octal);
                                 570

                                 1072

> 2*8^0 + 7*8^1 + 0*8^2 + 1*8^3;
                                 570

 

f:=-144*z-44+12*sqrt(-12*z^3+96*z^2+24*z-15);
subs(z=x+I*y, f);
abs(%);
plot3d(%, x=-3..3,y=-3..3, axes=boxed);


Edited:

Formally: if it has a root then bring the sqrt on one side and square it:
0=f;
% +144*z+44;
lhs(%)^2=rhs(%)^2;
lhs(%)-rhs(%);
expand(%);
#%/4; # gives Christian Wolinski's polynom
R:=solve(%);
       R := -4/3, -4/3, -4/3
      
The solution must be z = -4/3. But that is no solution:
eval(f, z=-4/3); simplify(%);
                                 296
NB: that is similar to ask for solving I = - sqrt(z).
2*ln(x+(x+ln(x))^(1/2))+2*(x+ln(x))^(1/2)
That will not work for several reasons, a better way is
H22*x*(cos(epsilon*((1/2)*L-x)/R)+kappa*cosh(epsilon*((1/2)*L-x)/R)):
int(%, x=0 .. L);
                        -30      2                      -12
    -2.95229265242607 10    omega  + 1.93138121006040 10  

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
First 10 11 12 13 14 15 16 Last Page 12 of 92