Axel Vogt

5821 Reputation

20 Badges

20 years, 228 days
Munich, Bavaria, Germany

MaplePrimes Activity


These are answers submitted by Axel Vogt


  assume(0<A, 0 <=t, 0 < xi, xi < 1); 
  map(getassumptions, [A,t,xi]);

    [{A::RealRange(Open(0), infinity)}, {t::RealRange(0, infinity)}, {xi::RealRange(Open(0), Open(1))}]


  evalc(A*exp(-sqrt(0.1-1)*t));

          A cos(0.9486832981 t) - A sin(0.9486832981 t) I


  evalc(A*exp(-sqrt(xi^2-1)*t));

                        2 1/2                   2 1/2
           A cos((1 - xi )    t) - A sin((1 - xi )    t) I

Be aware that 'evalc' alraedy assumes that variables are Reals, except you
assume something else (highlight that command (with the mouse) and press
F1 key in classical or F2 key in the other interface to see the help).
Here is (roughly) what I mean. It is easy to code it (with some care) and
only for x ~ 1 would need some extra twist (it works in 'any' precision).
And it is reasonable fast and flexible, if one is not willing to use what
Maple already provides.

theNewton:=proc (f,x0,n)
# searches for zero of the function f 
# x0 = initial guess, n = steps
 local x,x_new,fx,i;
 i:= 0; 
 x:= evalf(x0);
 while (true) do 
   fx:= evalf(f(x));
   x_new:= evalf(x-fx/D(f)(x));
   # one should check ^^^^^^^ first for (almost) zero
   if( abs(evalf(x-x_new)) < 1e-14  or n <= i) then
     #                       ^^^^^ or what one needs
     break
   end if;
   x:=x_new;
   i:= i+1;
 end do; 
 print(`steps performed`=i); 
 return x;
end proc; #maplemint(%);

Now use it (with 10 Digits):

### for 0 <= x < 1
  x:=1.0-1e-9; # 
  f0:= y -> y*exp(y) - x;

  someRoot:=theNewton(f0,x/2,20);
  LambertW(x);
  f0(someRoot);

                         steps performed = 5
                       someRoot := 0.5671432900
                             0.5671432900
                                     -9
                               0.1 10


### for 1 <= x < infinity
  x:=1e201;
  f1:= y -> ln(y)+y  - ln(x);

  someRoot:=theNewton(f1,ln(x),20);
  LambertW(x);
  f1(someRoot);

                         steps performed = 2
                       someRoot := 456.6955866
                             456.6955866
                                  0.

In x=1 a short Taylor polynom then is enough, i.e.

.205247033774894651+.361896256634889222*z-.107270323141071853*(z-1)^2

That's all I think (except one wants speed or else).

First you should ASCII for your file names.

And second: if you open that file with an text editor the you can guess it is an html file - it says "Page not found"

May be you want to post the task before writing down the path integral along the ellipse

PS: I think the correct restrictions (with a = 1) are

a := 1; # assume(a, real);
assume ( 0 < b, b < 1); #assume(b, real); additionally(1 > b); 
getassumptions(b);
assume(0<=p,p<2*Pi); assume(0<=q,q<2*Pi);
map(getassumptions, [p,q]);

since the integrand is periodic in p and q

And it may be seen as a real part (because of the cos(q) in front )

It sounds confusing and somewhat desperate ...

It is only because "... and then reaction stops", from which I guess it is
from Chemestry. But it really would best you state your actual task in some
*precise* way (even if it may be homework). Including the range of all your
variables (positive values?). Up to now it is a bit messy.

And - by the way - it is not a lumbered W :-)


For others: s = LambertW(exp) is crying for Wrightomega ...

If one does not need a fast solution (closed form approximation), then LambertW
can be done through Newton's method easily I think:

x = a*exp(a) looks good enough for 0 <= a <= 1 and beyond ln(x) = ln(a)+a
is fine enough to solve for a = LambertW(x). One only needs some reasonable
initial guess, either just brute (x = a/2 resp. x = exp(a)) or not so bad by
using the approximation shown at Wikipedia.

Have not looked at it, but a 2nd order Newton will be fast (seems Brigg's
code is towards that).

If you want to work with LambertW on 'paper' you may will have to learn something about that function (and series in 0 will not tell you that much).

Otherwise - if you want numerical values - you should say it (and then the series will cover only a small part).

So: what to you intend?

Let us look at just one of your integrals and if my guessing is correct
then you want the value of 

  Int( 1/abs(r-r3), r3=a..R);

Please compare that against what you have written down. 

Here I used the UpperCase, which Maple essentially lets just print the
integral, so you can look at it and how Maple understands your input.

And after that I use the command

  value(%);

which then is the same as using lowercase and results in the warning
you mentioned.

To tell Maple about your assumptions you can do the following:

  assume(a < r, r < R);
  Int(1/abs(r-r3),r3=a..R);
  value(%);

The answer is 'infinity' and to see what happens just try some data:

  plot(1/abs(1-r3),r3=0..2); # a=0, R=2 and r=1

You see a 'pole' at 1 and that integral does not converge: it is like
asking for Int( 1/x, x=0..1).

I do not understand all of your words (like "row" ...)

But it is my impression you want values and do not like the answer 'LambertW'.

If  it is so then however it is (in general) not the right way to ask for a Taylor series ...
and you should not hesitate to ask your actual question.

Sorry, if that is a wrong impression.

 

PS: that is an example, why it makes sense to provide a more complete profile,
often it is just the language (certainly in my case) or interest which may hint
others to give you an adequate reply

Using Maple 12 classic, Win XP SP2 32 bit with AMD (3500, 2.2 GHz,
1 MB) the timings are more or less the same, 0.063 and 0.031 (but
not constant in various runs). Whether as function or as proc I
always saw that speed up by choosing a fixed NAG method.

I did that with Digits:=14

Int(x^2/(exp(x)-1), x=0..infinity) - Int(x^2/(exp(x)-1), x=0..y);
is given in polylogs, which allows to check for y ~ 100.

The first integral as numerical constant = 2*zeta(3) can be used
to re-formulate the task in terms of the second integral, as it
was alraedy suggested, which speeds up by a factor 2.

I did not find more advantage by precomputing for 0 ... 4 and then
only using Int(..., x = 4 .. y ).

Also I did not find an advantage in defining the integrand as proc
and using option hfloat there.

 

Edited (after seeing Alec's edit ...): IIRC otherwise evalf@Int constructs a proc

Pondering about the task and the linked paper makes me feel, that it is
a kind of 'mixed' approach for valuating (call) options by binomial trees
but instaed of summing an integration is used.

The 'clustering' into groups after p.7 may have computational raesons.

If that is true, then there are at least two ways: pass to normal pdf as
the limit for binomial pdf or use the Beta function instead of binomial.

The correct pay-off for the 'option' is given on p. 9 (alternatives are
A and B, a fixed amount or a call option).

Then it reads like a (social) experiment to simulate the distribution
and calculate from there.

But it is just a thought and my feeling may be wrong.



Ps: acer's approaches made me to file this thread (and these techniques
reminds me of something ...), thank you - that is very fine!

It depends. But try GAMMA( x ) or GAMMA( ln(x) ) or ...  for x = some number with 7 or with 100 Digits.

I think you can not rely on 'order' (especially in a sum) - do you want to select a special summand?

http://www.mapleprimes.com/forum/helpplot3dx0

If I write E as Int with symbolic upper bound R and make change coordinates
by Change(  E(c0, c1, c2) , x = xi*R, xi) then I get an integral over unit
interval of the form Int( A(x)/B(x), x=0 ..1), xi = x

Here A(x) = A(-x) is an odd polynomial of degree=29 ,A(x) = x*A1(x^2)
and B(x) = B(-x) of form B(x) = B1(x^2)^(5/2) / sqrt(1-x^2).

This has the taste of a Chebyshev integration over half its range.

May be that helps for the integration routine.


For the specifically found & guessed configuration by acer:

  E(1, 0, 0):
  Change(%, x = xi*R, xi): subs(xi=x,%);
  value(%);

                             1
                            /
                           |        x
                        4  |   ----------- dx
                           |         2 1/2
                          /    (1 - x )
                            0


                                  4

Which does not depend on R = 3.91 ... hm ?
First: why should Maple find a symbolic solution?
Second: I always prefer 'exact' values over floats, so y:=1 and V:=2 and 
use functions (if not having specific reasons), hence

1+((abs(V))^2/W^2)*(CylinderD(z,y)^2)*(CylinderD(z,-y)^2-CylinderD(z,y)^2);
f:=unapply(%, z);

After staring and playing ... and converting to hypergeometrics and so
on I got aware, that you work on the cylinder function w.r.t. *degree*
and that is the reason, why Maple refuses to work:

Maple does not know the solution (to evaluate it numerically) and using
hypergeomtric sums it seems, that this would need series of Psi.

Without diving into theory (for which I am missing motivation) you can
use a pretty nice feature of Maple: numerical differentiation.

  g:=proc(t)
  fdiff(f(phi(T)),T=t)/f(phi(t)) * D(phi)(t);
  evalf(%);
  end proc;

This is a bit rude: g(0) and g(2*Pi) would not work, but g(Pi) and if
you plot that stuff then you can guess a bit more:

  plot( Re('g'(t)),t=0+0.01 .. 2*Pi-0.01);

This suggests symmetry at t=Pi and

  plot( Im('g'(t)),t=0+0.01 .. 2*Pi-0.01);

suggests point symmetry and that the integral over it is thus 0.

Finally I get

  Digits:=14;
  2*Int( Re@g, 0 .. Pi, method = _d01ajc, digits = 8); evalf(%);

                              2.7870136

which should be the result for your task with 8 decimal digits.
First 52 53 54 55 56 57 58 Last Page 54 of 92