kill(all)$

maxima:Euler-Mascheroni constant,limit,difference,version 1.0
peter.vlasschaert@gmail.com,15/03/2022

ref :
https://nl.wikipedia.org/wiki/Constante_van_Euler-Mascheroni
https://en.wikipedia.org/wiki/P_series

load packages
------------------

load(lrats)$

definition : Euler-Mascheroni constant = γ
---------------------------------------------------

p1:'integrate(exp(-x)*log(x), x,0,inf)=integrate(exp(-x)*log(x), x,0,inf);

→ numerical value : -γ

p1a:(-1)*(rhs(p1)=float(rhs(p1)));/* approximation of :γ */

definition : 'exp(x) @  x=1' = e
-------------------------------------

→ limit : inf,find = %e

p2: 'limit ((1+1/y)^y,y,inf) =limit ((1+1/y)^y,y,inf);

p2x:args(lhs(p2));

→ numerical value : e

p2a:float(p2);

p2b:rhs(p2) = float(rhs(p2) );

→ limit : 0, find : %e

p2c:psubst([1/y=z,y=1/z,∞=0],p2x);

(p)subst : parallel substitution.

listofvars(p2c[1]);

p2c1:limit(p2c[1],%[1],p2c[3])='limit(p2c[1],%[1],p2c[3]);

what about : limit below

p3:'limit ((1+a/y)^y,y,inf) =limit ((1+a/y)^y,y,inf);

p3a:subst([a=1],p3);/*see p2*/

find : exp(-x),use p3,needit for integration

p4:subst([a=-x],p3);

p4a:'limit ((1+a/y)^(y-1),y,inf) =limit ((1+a/y)^(y-1),y,inf);

why : ∞ = ∞-1 

p4b: reverse(subst([a=-x],p4a));

p4b1:args(lhs(p4b))=args(rhs(p4b));

first step : p4b into lhs(p1)
-------------------------------

p5:rhs(p4b1)[1];

p5a:psubst([lhs(p4b)=p5,∞=y],lhs(p1));

second step : (1-x/y)=u
----------------------------

menu : Calculus →Change variable
Integral : p5a
old   variable = x
new variable = u
equation        : u=1-x/y

p5b:changevar(p5a, u=1-x/y, u, x);

p5b1:expand(p5b);

p5b2:args(p5b1);

assume(y-1>0);

first part integral :

p5b2a1:'integrate(part(p5b2[1],1),p5b2[2],p5b2[3],p5b2[4])=
integrate(part(p5b2[1],1),p5b2[2],p5b2[3],p5b2[4]);

second part integral :

→ first way : maxima (direct computation)

p5b2a2:'integrate(part(p5b2[1],2),p5b2[2],p5b2[3],p5b2[4])=
integrate(part(p5b2[1],2),p5b2[2],p5b2[3],p5b2[4]);

→ second way : use taylor expansion

menu : Calculus → Get series → series
Expression: log(1-u)
Variable    : u
Point         : 0
select :powerseries

p5b2b2:log(1-u)=niceindices(powerseries(log(1-u), u, 0));

p5b2b2a:subst([p5b2b2],lhs(p5b2a2));

p5b2b2a1:intosum(p5b2b2a);

p5b2b2a2:args(args(p5b2b2a1)[1]);

p5b2b2a3:args(args(args(p5b2b2a1)[1]))[1];

p5b2b2a4:args(args(args(args(p5b2b2a1)[1]))[1]);

rem : intosum

p5b2b2a5:intosum(p5b2b2a4[1]);

p5b2b2a6:subst([∞=y],%);

p5b2b2a7:factor(integrate(p5b2b2a6,u,0,1));

partfrac(sum(1/(i*(y+i)),i,1,y), i);

partfrac(1/(i*(y+i)), i);

p5b2b2a8:p5b2b2a1=intosum(p5b2b2a7*p5b2b2a2[2]*(-1));

p5b2b2a9:partfrac(part(rhs(p5b2b2a8),1),i);

declare(sum,linear)$

p5b2b2a10:(sum(p5b2b2a9,i,1,inf));

p5b2b2a11:subst([y=∞],part(p5b2b2a10,2,1));

simp:false;

rem : otherwise direct evaluation.

p5b2b2a12:p5b2b2a11=
sum(part(p5b2b2a11,1),i,1,v)+sum(part(p5b2b2a11,1),i,v+1,∞);

p5b2b2a13:part(rhs(p5b2b2a12),2)=sum(1/(v+i),i,1,inf);

p5b2b2a14:subst([p5b2b2a13],p5b2b2a12);

p5b2b2a15:lhs(p5b2b2a14)=psubst([∞=v,v=y],rhs(p5b2b2a14));

p5b2b2a16:subst([p5b2b2a15,v=inf],p5b2b2a10);

because : before use : simp:false

simp:true;

p5b2b2a17:p5b2b2a16;

final result :integral "p1" = first integral+second integral

p6:p1=limit(rhs(p5b2a1)+p5b2b2a17,y,inf);

p6a:(-1)*(rhs(p1)=rhs(p6));

γ = limit    (H(i,y)-ln(y))
       y→∞

H(i,y) = Harmonic series always, divergence , p series p =1

result : maxima output
first integral :

p6b1:p5b2a1;

second integral:

p6b2:p5b2a2;

p series:
----------

p7:sum(1/i^p,i,1,y);

p7a:'limit(p7,y,inf) =subst([y=inf],p7);

simpsum;

simpsum:true;

p7b:divg;/*subst([p=1],rhs(p7a))*/;

error : -- an error. To debug this try: debugmode(true);
sum : sum divergent

example :p=2, convergence
---------------------------------

p7b1:'subst([p=2],rhs(p7a))=subst([p=2],rhs(p7a));

example :p=1/2, convergence

p7b2:'subst([p=1/2],rhs(p7a))=subst([p=1/2],rhs(p7a));

 load(simplify_sum)$

p7b3:'subst([p=1/2],rhs(p7a))=simplify_sum(subst([p=1/2],rhs(p7a)));

result : p series
------------------

rem : p > 1     : cvg:convergence
          0<p ≤ 1 : dvg:divergence

---------------------------------------------------------------------------------

p8a:0<'integrate(1/x^p,x,1,inf);

p8a1: M1:rhs(p8a);

p8a2:M1< sum(1/n^p,n,1,inf);

pa3:M2:rhs(p8a2);

pa4:pa3< 1+ rhs(p8a);

info : about the p:series