Salvo, no problem, over all the years I get used to that, but only almost :-)

Salvo,

double exponential integration is not my method, it is due to Mori and known to work well even for higher precision. All I did was translating/adapting existing code to Maple, as stated.

Roughly it computes Int(f(x), x = -1 .. 1) by transforming it through x = tanh( Pi/2 * sinh(t) ) to the full real line and using the trapez rule there.

PS: Alex is not quite the same as Axel ...

I lost a bit the overview for the variants you have in mind, but hope the
attached sheet does it for your sub-task Int(t^a/S(t), t = b ... c).

A bit strange, that Maple's on-board-method double exponential integration
seems to be on strike for that. So I modified what I once had translated
for my needs, hope it helps.

For the jittering between levels: it is a 'zooming' or scaling problem and
does not stand for errors - it more or less says that the 'resulution' is
not fine enough (and one can be trapped by graphics as well). May be you
also consider the advice 'make it dimensionless' (but as a former Math guy
I am not used really to that recipe).

The last thing i tried (not today) was towards approximating the integrand
by a Hermite expansion. But your stuff as pdf (i.e. with a=0) is quite far
away from a normal distribution (even on that finite range and certainly
not over the Reals). But as 'zooming' is your problem one would need quite
a good one (even if symbolic solutions exist). So I did not: time for a beer.

Download 102_double_exponential_integration_with_varaible_digits.mws
View file details

Edited 29 Jan 2010

The routine has a typo: at the beginning Pi/2 is computed as
pi2 := evalhf(1/2*Pi) and that has to be pi2 := evalf(1/2*Pi) to
have that with the correct precision.

seems to be too high ...

In your case however it is easy:

((x^7+1)^4/x^30); expand(%);

                     1     4     6     4      1
                    --- + --- + --- + ---- + ----
                     30    23    16     9      2
                    x     x     x      x      x

from which you can read it off or continue by

[op(%)];
map('t -> residue(t, x=0)', %);
convert(%, `+`);

                                  0

I do not know an actual answer to your question and I use M12.

a) Maple can be convinced to give a 'general' anwer

  Int(1/(x^4+3*x^2+1)^(p),x=0..infinity);
  Change(%, x=1/2*(-6+4*t)^(1/2),t); 
  # motivated by the symmetry of the *real* partial fraction decomposition
  # x^2+3/2=t, certainyl invertible on Positives
  value(%); combine(%);
  R:=%;

I omit to print the result, but Maple seems not to use reflection formula
to go further (I think it should be possible)

b) Now with Digits:=14 compare

  pTst:=1 + 1e-10;
  'eval(R, p=pTst)'; evalf(%);
  'Int(1/(x^4+3*x^2+1)^(p),x=0..infinity)': 'eval(%, p=pTst)'; evalf(%);
 
                           0.70258767772836

                           0.70248147304481

Is it really only correct for 3 leading decimals?

c) What happens here? Certainly a bug ...

  convert(R, StandardFunctions) assuming p::integer;

                                  0

I do not understand this today ...

It also may depend on your 'expression'. If it is not to complicated then it is better to post it as well.

You do not provide the situation, where the error occurs (the linked sheet works for me).

From the error message I would guess: you have negative values for x.

In your other post it seems that it can not be solved for beta ...

Let me write it in a more readable form, using x as the variable:

sqrt(Pi)*beta=-1/h*SteL/(exp(h^2*beta^2)*erf(h*beta))+
  R1*SteV/(exp((beta/R1)^2)*erfc(beta/R1));
subs(SteL=a,SteV=b, R1 = r, beta=x, %): 
B0:=combine(%, exp);
B0/sqrt(Pi)/x; 
B1:=collect(%, exp);

            1 = (some complicated stuff in exp and erf)

theAssumption:=a<0, b<0, h>0, r>0, x>0;

Now staring at it one can guess the pattern:

1 = '-b/2*D(erfc)(u)/erfc(u)/u - a/2*D(erf)(v)/erf(v)/v'; 
subs(v=h*x, u=x/r, %);
% = (B1);
is(%);
                                true

You can put the factor 1/2 into your a and b to have it nicer,
but that equation does not cry for an explicite solution in x.
I do not even see a series approximation.

So my feeling says: you need numerical solutions. 

Computing that will need some care: in x=0 you have a pole,

  B1; rhs(%);
  series(%, x=0, 5); because of D(erf)/erf

    (-1/2*a/h^2)*x^(-2)+b/Pi^(1/2)*r*x^(-1)+(1/3*a+2*b/Pi)+O(x)

and for x ---> oo that converges to b:

  B1; rhs(%); 
  eval(%, erf = erf(infinity)); # it gives 1
  asympt(%, x, 4) assuming 0 < r, 0 < h;
  lprint(%);

                     b+1/2*r^2*b/x^2+O(1/x^4)

but I have not looked inbetween (you should: if it is decreasing,
then there can not be a solution, if 1 < b=SteV ...).

But as alraedy said: I doubt there is a symbolic solution and it
is not even clear, whether there always exists a (real) one.

A:=Pi^(1/2)*beta = 1.878297101/exp(20.00000001*beta^2)/erf(4.472135956*beta)-
    2.100000000*R1/exp(beta^2/R1^2)/erfc(beta/R1);
identify(%): convert(%, rational):
A:=lhs(%) - rhs(%);

plots[implicitplot](A, R1=0..1, beta=0..1);


Now determine the 'extreme' points you are interested in (guessing from the graphic):

0=limit(A, R1=0, right) assuming 0 < beta;
expand(%);
beta1:=fsolve(%, beta=0.15);

                        beta1 := 0.1978948870

eval(A, R1=1);
beta0:=fsolve(%, beta=0.15);

                        beta0 := 0.1184172673


Then solve the equation 'formally' for R1 (not for beta, that gives troubles)
and consider that as a function of beta (should work due to the graphic):

RootOf(A, R1);
g:=unapply(%, beta);

Your graphic becomes much more smooth (it is the reflected one, correct that),
thus the original graphic would give quite error prone results:

plot(Re(evalf('g'(t))), t=beta0 .. beta1, numpoints=5);
plottools[reflect](%, [[0,0],[1,1]]);


If you want values then for example you can use

myStepSize:=0.01;
[seq( [t, g(t)], t=beta0 .. beta1, myStepSize)]:
evalf(%); 
                          myStepSize := 0.01


  [[0.1184172673, 1.000000000], [0.1284172673, 0.8441382672],

        [0.1384172673, 0.7054733505], [0.1484172673, 0.5804226649],

        [0.1584172673, 0.4661017791], [0.1684172673, 0.3598910320],

        [0.1784172673, 0.2587389755], [0.1884172673, 0.1567487861]]

You can find a recipe here: www.math.rwth-aachen.de/mapleAnswers/html/313.html

PS: it is worth to look at this pages beyond that question (even if it a bit old, some stuff may be outdated)

you can look at that in the online help (at the top of this page you find a search, activate 'help' and enter "DiscreteTransforms[FourierTransform]"

www.maplesoft.com/support/help/AddOns/view.aspx

The following might do, if you work with 'Standard' interface

plot3d(-x^3-y^2, x=-1..1, y=-1..1, axes=boxed, style=point);

If you work with the classical interface, then the graphics are less fine,
points can not be 'sized' (at least for me), but other symbol should do:

plot3d(-x^3-y^2, x=-1..1, y=-1..1, axes=boxed, style=point, 
  symbol=circle,  symbolsize=15, numpoints=20^2);

After simplifications this is a very special case, which Maple does not answer:

Sum((k+j)^j*r^j/j!,j = 0 .. infinity); eval(%,k=1);

                         infinity
                          -----          j  j
                           \      (j + 1)  r
                            )     -----------
                           /          j!
                          -----
                          j = 0

I googled for the keyword, seems somewhat special, both in combinatorics and other fields (epidemiology, queing ...)
and now think, that you have to invest some theory for that topic: search for "arch98v127.pdf" gives me the impression,
that your formula stems from a Lagrange inversion. And that cries to be treated before using it.

a) why?

b) what do all elements have in common and what can you operationally say for that?

There is no minimum. Even if I suppose it is homework:

The last terms are linear, where x,y,z live in a cube, thus this gives
something bounded (IIRC Analysis correctly).

So your 1st term matters. Transform to positive x,y,z:

  subs(z=-z,%);  subs(y=y+1, %);

    a*b^x*c^(y+1)*d^(-z)-x-2*y-2+3*z

now with 0 < z < 1, 0 < y < 2

This suggests you want b,c,d to be positive to avoid imaginary values
(may be you forgot to say that).

Now choose: subs(x=1/2, y=1, z=1/2, %);

This gives a*somePositiveConstant - anotherConstant

So for a ---> - infinity ...

Or similar ...