d:=401/2 - 1880/(3*Pi):
op([1,1..2],d);
op([2,1,1..2],d);
                                401, 2
                               -1880, 3

Applying the definition you are using directly the primitive/antiderivative function:

F:=int(1/(4*sin(x)-3*cos(x)), x);

F := 1/5*ln(3*tan(1/2*x)-1)-1/5*ln(tan(1/2*x)+3);

evaluated at the limits of each subinterval. But the integration "constant" has a jump in its imaginary part at arctan(3/4):

plot(Im(F),x = -(1/2)*Pi ..(1/2)*Pi);

I guess that the code for 'CauchyPrincipalValue' somehow is able to handle this discontinuity.

I wonder whether the code for definite integration includes a check on whether the integrand does not change sign within the interval of integration.

For  an integrand  that  is an 'abs'  call  and some simple variations it  is trivial.  But  for more complex  functions  it could be  difficult  (if not  impossible).  For this reason it might be  that there is no such check.

 

You may look at the thread Viewing function as a table.

The command

INTERFACE_WORKSHEET(display,file="path/file.mw");

opens the worksheet file.mw in a new tab. Presumably, you could add several such lines in your maple.ini file and get those worksheets opened in their tabs.

This integral can be calculated using a double regularization (I think it can only be done through a double limit). That is, I choose to regularize the integral

J:=Int((1/(exp(x)-1)-1/x+1/2*exp(-x))/x,x=0..infinity);

through the introduction of a generalized integral K that depends on a pair of parameters (alpha,beta) in this way:

K:=Int((1/(exp(x)-1)-1/x+1/2*exp(-x))*exp(-alpha*x)/x^beta,x=0..infinity);

K := Int((1/(exp(x)-1)-1/x+1/2*exp(-x))*exp(-alpha*x)/(x^beta),x = 0 .. infinity); The integral J is recovered for alpha=0, and beta=1. Then I expand in the first term 1/(exp(x)-1) in powers of exp(-x), which converges for x>0:

with(IntegrationTools):
f:=GetIntegrand(J);
fn:=op(1,f);
fn1:=op([1,1],f);
convert(subs(exp(x)=u,fn1),FormalPowerSeries,u=infinity);
subs(u=exp(x),%);
fne:=subs(fn1=%,fn);
ge:=fne*exp(-alpha*x)/x^beta;

Then combine the exponentials:

expand(ge);
gec:=combine(%,exp);
suma:=numer(select(has,gec,Sum));
Tk:=combine(op([1,1],%)*op(2,%),exp);algsubs(suma=Sum(Tk,k=0..infinity),gec);
gecc:=expand(%,exp);

Now I integrate term by term in a region of the parameters where the integral converges, and form the partial sum for k=0..M

op([1,1],gecc)*op([1,2,1],gecc);
Jk:=int(%,x=0..infinity) assuming alpha>0,beta>0,k::posint;
UM:=sum(%,k=0..M);

UM := sum((k+alpha+1)^(-1+beta)*GAMMA(1-beta),k = 0 .. M); and take the limit M-> infinity :

U:=limit(UM,M=infinity) assuming alpha>0,beta<0;

U := (-GAMMA(-beta)*beta*Zeta(0,1-beta,2+alpha)*alpha-Zeta(0,1-beta,2+alpha)*GAMMA(-beta)*beta-(alpha+1)^beta*GAMMA(-beta)*beta)/(alpha+1); Now I integrate the second and third terms in the same region:

Int(op(2,gecc)+op(3,gecc),x=0..infinity);
V:=value(%) assuming alpha>0,beta<0 ;

V := -alpha^beta*GAMMA(-beta)+1/2*(alpha+1)^(-1+beta)*GAMMA(1-beta); And finally expand the sum of terms in series about beta=1, and take the limit for alpha -> 0:

series(U+V,beta=1,2):
 convert(%,polynom);
 limit(%,alpha=0);
 combine(%,ln);

-1/2*ln(2*Pi)

to a "property expression" type, extracted from the code of 'assume' in Maple 11.02:

`type/propexp`:={'`or`', '`not`', '`and`', '`in`', 
'`::`', 'relation', ('specfunc')('anything',{'Or', 
'And', 'Non', 'Not', 'AndProp', 'OrProp'})};

type(a>0,propexp);
                                 true
type(Not(x::posint),propexp);
                                 true
type(n::posint,propexp);
                                 true
type(4,propexp);
                                false
type(sin(x),propexp);
                                false

But not completely right, as eg:

type(4>0,propexp);
                                 true

The editor has options to edit size  through the Insert/Edit icon.

as that example code worked well with Maple 10.00 and produces this error message since, at least, Maple 10.03.

With Maple 9.52:

sol := {T(y,z) = _C1*exp(_c[1]*y)*_C2/z*WhittakerM(1/4*I*_c[1]^(1/2),0,_c[1]^(1/2)*z^2*I)+_C1*exp(_c[1]*y)*_C3/z*WhittakerW(1/4*I*_c[1]^(1/2),0,_c[1]^(1/2)*z^2*I)}; Not checked.

which is free, you can get the other one:

So, you may want to check results.

The packages 'LinearAlgebra' (or 'linalg') only do matrix linear algebra for linear vector spaces of finite dimension. But, apparently, there is no package to do just linear algebra.

Note, the expression "abstract linear algebra" is sometimes used to mean just this plain linear algebra, as opposed to "concrete lineal algebra" that is a matrix representation for a linear space of a given finite dimension. So, in this context, do not be afraid of the expression...

But I think that it should not be very difficult to make up something with a module, exporting suitable `diff`, `*`, etc.
 

A remembrance.

Reading the signature in the original post, it reminds me that it was at the CBPF, in the inaugural lecture of a conference for relativists, that I have seen for the first time a CAS in action. It was Reduce in mid 1987.

Instead of the original integral, I will consider the integral:

J:=Int(exp(-a*x^2)*sech(x/2)^2, x=0..infinity);

where a>=0. Due to the even symmetry of the integrand, the integrand on (-infinity,infinity) is twice this one. I am not aware of a closed form for this integral, except for a=0:

eval(J,a=0);
value(%);

2

So, I will seek for a series expansion. Clearly, different approaches are possible. For this one I will use commands from some packages:

with(IntegrationTools):
with(MultiSeries):
with(SeriesInfo):
with(gfun):
with(plots):

First, I convert the 'sech' part to exponentials and extract the integrand:

convert(J, exp):
expand(%):
Je:=Combine(%);
fe:=GetIntegrand(Je);

fe := 4/(exp(x)+1)^2*exp(x)/exp(a*x^2); On one hand, this integrand is bounded, rapidly convering to 0 for x->infinity. On the other hand, gaussian integrals with just an exponential factor have a closed form. Hence I will expand the 'sech' part in powers of exp(-x). That is I will use the power expansion of (1+u)^(-2), where u=exp(-x)0 (hence convergent).

subs(exp(x)=u,fe);
s:=multiseries(%,u=infinity);

s := 4/exp(a*x^2)/u-8/exp(a*x^2)/u^2+12/exp(a*x^2)/u^3-16/exp(a*x^2)/u^4+O(1/(u^5)); What I need now is the general term of this expansion. One way to do it is:

L:=ListCoefficients(s);
rec:=listtorec(L,v(n));
rsolve(op(1,rec),v(n));
A:=unapply(%,n);

So that I get the generic term 'Tn' in the expansion of the integrand:

Tn:=A(n)/exp((n+1)*x);

Tn := 4*(-1)^n*(n+1)/exp(a*x^2)/exp((n+1)*x); and I integrate term by term (presumably, as everything is so smooth and convergent, it should be fine) to obtain the generic term 'Jn' of the expansion of the integral:

Int(Tn, x=0..infinity):
value(%) assuming a>0:
factor(%):
Jn:=convert(%,erfc);

Jn := -2*(-1)^(n+1)*Pi^(1/2)*exp(1/4*(n+1)^2/a)*(n+1)*erfc(1/2*(n+1)/a^(1/2))/a^(1/2); Here, I have converted to the complementary error function because with the 1-erf form the numeric errors soon became too large as the argument alpha=1/2*(n+1)/a^(1/2) grows at about 8. Well, this expression is not so nice in regards to convergence as I would have expected. Moreover, does this series expansion converge? I want to see the behavior of abs(Jn) for large n at fixed a. Or equivalently, its behavior for alpha -> infinity:

abs(Jn) assuming a>0,n::posint;
Ja:=algsubs(1/2*(n+1)/a^(1/2)=alpha,%);

Ja:=4*Pi^(1/2)*exp(alpha^2)*erfc(alpha)*alpha;

Jm:=series(%,alpha=infinity,4);

Jm := 4-2/alpha^2+3/alpha^4+O(1/(alpha^5)); Danger! the leading term produces an undefined sum of alternating terms (-1)^n*4, with divergent subsums of even and odd terms. But this is not unusual for me. I do not know how mathematicians will think of what is going on next, but physicists are rather used to deal with diverging magnitudes and use tricks to extract finite parts: the "renormalization" procedure. In fact, I think that there is a connection with the transseries that Jacques has pointed out before (I have not finished reading that paper yet). But I could hardly found such a connection in the literature. So, I will put all the dirt within an additive constant C to be fixed using a suitable constraint. That is, I will say that in some sense: J = Sum(-2*(-1)^(n+1)*Pi^(1/2)*exp(1/4*(n+1)^2/a)*(n+1)*erfc(1/2*(n+1)/a^(1/2))/a^(1/2), n= 0 .. infinity)+C; The constraint that I will use is the value for a=0, ie J(0)=2. One way to use this constraint is to define the partial sums

Jsp:=unapply(Sum(Jn,n=0..M),a,M);

and "identify" the limit of these partial sums as a->0. But you see, the partial sums alternate sign, depending on whether M is even or odd. Eg for a=1, M from 10 to 20:

Digits:=10:
seq(evalf(value(Jsp(1,i))),i=10..20);

 2.767448625, -1.179163813, 2.775087537, -1.185286044, 2.780066841,
 -1.189388433, 2.783485632, -1.192266726, 2.785931105,
 -1.194361647, 2.787739149;

The limit a->0, fixed n, corresponds to the limit alpha->infinity in Ja. Hence, for even M, it is of the form 4+(-4+4)+...+(-4+4)=4, while for odd M it is 4+(-4+4)+...+(-4+4)-4=0. Numerically:

Digits:=20:
seq(evalf(value(Jsp(10^(-i),10))),i=1..8);

  3.5803451145890104598, 3.9380316837243077610, 3.9934351793273425299,
        3.9993394728060958994, 3.9999339064143410016,
        3.9999933902323510823, 3.9999993390186819947,
        3.9999999339018734719

seq(evalf(value(Jsp(10^(-i),11))),i=1..8);

  -0.4141223187961680176, -0.0614129920410232268,
        -0.0065092674317558932, -0.0006549716614965324,
        -0.0000655380303248472, -0.65542120956771 10  ,
                           -6                    -7
        -0.6554247547842 10  , -0.655425709726 10

Correspondingly, there are two values of the renormalization constant C depending on the parity of M: 4+Ceven=2, ie Ceven=-2; and Codd=2. Now we can compare these renormalized partial sums with the numerically calculated integral. I would expect that the partial sum for a given M provides a better approximation when alpha(M)=1/2*(M+1)/a^(1/2) is large, ie. a is small. And this seems to be what can be observed in this plot for M=10:

Digits:=20:
p1:=plot(evalf(value(Jsp(a,10))-2),a=0.1..10,color=red):
p2:=plot(evalf(J),a=0.1..10,color=blue):
display(p1,p2);

where red is the partial sum and blue the numerical integral and C=Ceven=-2. I find this coincidence quite reasuring that despite the "magics" this calculation makes sense. Similar good coincidence can be obtained for odd M partial sums. In regards to the opposite limit a->infinity, the leading term:

series(Jn,a=infinity,1);
Sum(convert(%,polynom),n=0..infinity);

Sum(-2*Pi^(1/2)*(-1)^(n+1)*(n+1)*(1/a)^(1/2),n = 0 .. infinity); could be "calculated" as in:

Sum((-1)^n*(n+1)/u^n,n = 0 .. infinity);
value(%);
eval(%,u=1);

1/4

This is, apparently, a transseries or antilimit calculation (I do not know which denomination is better). So, in this sense, yes, the limit of the series is 1/2*sqrt(Pi/a). The next step is relate to the other approaches. About handling with new primes editor. Beyond simple posts with a few lines of text, I find much more convenient to edit in my favorite editor, using the tags and paste in source mode.

I have got this series expansion for the integral on the interval (0,infinity):

Sum(-2*(-1)^(n+1)*Pi^(1/2)*exp(1/4*(n+1)^2/a)*(n+1)*erfc(1/2*(n+1)/a^(1/2))/a^(1/2),n = 0 .. infinity)

It holds for all values of a>0.

I do not have time now to give the details. Sometime tomorrow.

Its partial sum of the first 30 terms compares fine with the numerical integral for small a (up to about 10), using Digits:=20 or so.

Well almost. A factor 2 has to be substracted. Kind of "renormalization" of its value for a=0.

writing lambda and alpha in terms of beta.

eq1:= tan(beta)+2/alpha*(beta);
eq2:=beta= sqrt(lambda*alpha-(alpha/2)^2);
eq3:= eval( eq2, isolate(eq1,alpha) );
solve( eq3, lambda );
simplify( %, trig );
lambda1 :=combine(%,trig);
alpha1:=solve(eq1,alpha);
plot([alpha1,lambda1,beta=0..Pi/2],view=[-2..0,-10..0]);