Axel Vogt

5018 Reputation

19 Badges

15 years, 305 days
Munich, Germany

MaplePrimes Activity


These are answers submitted by Axel Vogt

Strange ... I switched to exact values for S[j] and if I use w1 = 1/sqrt(1-x^2)
as weight and change to x = arcsin(t) then it plots pretty smooth (formally there
is a singularity in x=0), now over -Pi/2 .. Pi/2.

I would expect Maple to find it, at least after splitting in 0.

You could convert to a string and take its length
A:=x+2*y;
convert(A, string); length(%);
                                  5
To get rid of the length of variable names you can use LeafCount
B:=subs(x=theOneVariable, y=theOtherVariable, A);
               B := theOneVariable + 2 theOtherVariable
length(A), length(B);
                                9, 37
MmaTranslator[Mma][LeafCount](A), MmaTranslator[Mma][LeafCount](B);;
                                 5, 5

You can export through the menu bar File / Export as in various ways, HMTL, text formats (plain text or RTF) or similar through Save As (using 'Maple Input' generates a file *.mpl which is pure text) 

 

PS: Delete output within the sheet in case you want to see the input commands. That can be done through the menu bar Edit / Remove Output

There are two things.
In your integration over MI you mixed up the variables, try it step by step:
  H:= z -> evalf( Int( x -> MI(x,z), 0..3,method=_d01akc, epsilon=1e-8) );
  evalf( Int(H, 0 .. 2, method=_d01akc, epsilon=1e-8) );
                           1.43166543880904
Which differs a lot from the true value 24/Pi^2 ~ 2.413...
But the approximation is quite bad
  plot3d(sin(Pi*x/Lx)*sin(Pi*z/Lz) - 'MI'(x,z), 
    x=0..3,z=0..2, axes=boxed, title="absolute errors");
I used 'quotes' for MI to enforce simple plotting over a range which
is more convenient than feeding the approximation points.

PS: this is 2-dim integration, you integrate over 2 variables
tmp:=x^2+y^2 - 1:
indets(tmp, symbol); # ok, gives only w
indets(tmp);         # so you need 0 < w because of roots
plot(abs(tmp));      # w should be close to zero
plot(tmp, w=0 .. 1, 0 .. 10 );
Now try
fsolve(tmp, w=0 .. 1);
                       0.0240992450396199
Here is a very simple version for the bug: 1 / linear (please also see Carl's answer,
http://www.mapleprimes.com/questions/212935-Computing-A-Complex-Residue#answer228699).
So I guess the bug is not due to 'residue' itself, but some simplifications
or transformations.
 restart:
 g:=1/(z-sqrt(-1+I));
                                     1
                         g := ---------------
                                          1/2
                              z - (-1 + I)
 z0:=2^(1/4)*exp(3*Pi*I/8); 
 z1:= simplify( convert(%, radical) );
                             (1/4)
                      z0 := 2      exp(3/8 I Pi)

                    1/4        1/2 1/2           1/2 1/2
                   2    ((2 + 2   )    I + (2 - 2   )   )
             z1 := --------------------------------------
                                     2
 #eval(g, z =z1); radnormal(%); # div by zero -> z1 = pole
 z2:=sqrt(-1+I);
                                        1/2
                          z2 := (-1 + I)
 '0=z2 - z1'; evalc(%): simplify(%): is(%);
 '0=z2 - z0'; evalc(%): simplify(%): is(%);
                             0 = z2 - z1

                                 true

                             0 = z2 - z0

                                 true
 'residue(g, z = z0)': '%'=%;
 'residue(g, z = z1)': '%'=%;
 'residue(g, z = z2)': '%'=%;
                        residue(g, z = z0) = 0
                        residue(g, z = z1) = 0
                        residue(g, z = z2) = 1
 'residue(1/(z-a), z = a)': '%'=%;
                                1
                      residue(-----, z = a) = 1
                              z - a
 
PS: one can find that through looking at the partial fractions
and checking for which summand it fails
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
4 5 6 7 8 9 10 Last Page 6 of 86