Advanced Search

vv

13143 Reputation

20 Badges

9 years, 79 days
Contact vv

MaplePrimes Activity


These are Posts that have been published by vv

A limit from an undergraduate competitio...

vv 13143 Maple
limit int
May 17 2018
1 8
 

At a recent undegraduate competition the students had to compute the following limit

 

> 

Limit( n * Diff( (exp(x)-1)/x, x$n), n=infinity ) assuming x<>0;

Limit(n*(Diff((exp(x)-1)/x, `$`(x, n))), n = infinity)

 (1)

 

Maple is able to compute the symbolic n-fold derivative and I hoped that the limit will be computed at once.

Unfortunately it is not so easy.
Maybe someone finds a more more straightforward way.

 

> 

restart;
> 

f := n * diff( (exp(x)-1)/x, x$n );

n*(-1/x)^n*(-GAMMA(n+1)+GAMMA(n+1, -x))/x

 (2)
> 

limit(%, n=infinity);

limit(n*(-1/x)^n*(-GAMMA(n+1)+GAMMA(n+1, -x))/x, n = infinity)

 (3)
> 

simplify(%) assuming x>0;

limit(-(-1)^n*x^(-n-1)*n*(GAMMA(n+1)-GAMMA(n+1, -x)), n = infinity)

 (4)

 

So, Maple cannot compute directly the limit.

 

> 

convert(f, Int) assuming n::posint;

-n*(-1/x)^n*(-x)^(n+1)*GAMMA(2+n)*(Int(exp(_t1*x)*_t1^n, _t1 = 0 .. 1))/((n+1)*(Int(_k1^n*exp(-_k1), _k1 = 0 .. infinity))*x)

 (5)
> 

J:=simplify(%)  assuming n::posint;

n*(Int(exp(x*_k1)*_k1^n, _k1 = 0 .. 1))*GAMMA(n+1)/(Int(_k1^n*exp(-_k1), _k1 = 0 .. infinity))

 (6)
> 

L:=convert(J, Int) assuming n::posint;

n*(Int(exp(x*_k1)*_k1^n, _k1 = 0 .. 1))

 (7)
> 

L:=subs(_k1=u, L);

n*(Int(exp(x*u)*u^n, u = 0 .. 1))

 (8)

 

Now it should be easy, but Maple needs help.

 

> 

with(IntegrationTools):
> 

L1:=Change(L, u^n = t, t) assuming n::posint;

Int(exp((x*t^(1/n)*n+ln(t))/n), t = 0 .. 1)

 (9)
> 

limit(L1, n=infinity);  # OK

exp(x)

 (10)
> 

####################################################################

Note that the limit can also be computed using an integration by parts, but Maple refuses to finalize:

> 

Parts(L, exp(u*x)) assuming n::posint;

n*(exp(x)/(n+1)-(Int(u^(n+1)*x*exp(x*u)/(n+1), u = 0 .. 1)))

 (11)
> 

simplify(%);

n*(-x*(Int(u^(n+1)*exp(x*u), u = 0 .. 1))+exp(x))/(n+1)

 (12)
> 

limit(%, n=infinity);

limit(n*(-x*(Int(u^(n+1)*exp(x*u), u = 0 .. 1))+exp(x))/(n+1), n = infinity)

 (13)
> 

value(%);  # we are almost back!

limit(n*((-x)^(-n)*(-(n+1)*n*GAMMA(n)/x-(-x)^n*(x-n-1)*exp(x)/x+(n+1)*n*GAMMA(n, -x)/x)+exp(x))/(n+1), n = infinity)

 (14)
> 

 

add, floats, and Kahan sum...

vv 13143 Maple
evalf add
April 27 2018
3 2

add, floats, and Kahan sum

 

I found an intresting fact about the Maple command add for floating point values.
It seems that add in this case uses a summation algorithm in order to reduce the numerical error.
It is probably the Kahan summation algorithm (see wiki), but I wonder why this fact is not documented.

Here is a simple Maple procedure describing and implementing the algorithm.

 

 

> 

restart;
> 

Digits:=15;

15

 (1)
> 

KahanSum := proc(f::procedure, ab::range)  
local S,c,y,t, i;      # https://en.wikipedia.org/wiki/Kahan_summation_algorithm
S := 0.0;              # S = result (final sum: add(f(n), n=a..b))
c := 0.0;              # c = compensation for lost low-order bits.
for i from lhs(ab) to rhs(ab) do
    y := f(i) - c;     
    t := S + y;              
    c := (t - S) - y;        
    S := t;                  
od;                         
return S
end proc:

 

Now, a numerical example.

 

 

> 

f:= n ->  evalf(1/(n+1/n^3+1) - 1/(n+1+1/(n+1)^3+1));

proc (n) options operator, arrow; evalf(1/(n+1/n^3+1)-1/(n+2+1/(n+1)^3)) end proc

 (2)
> 

n := 50000;
K := KahanSum(f, 1..n);

50000

  

.333313334133301

 (3)
> 

A := add(f(k),k=1..n);

.333313334133302

 (4)
> 

s:=0.0:  for i to n do s:=s+f(i) od:
's' = s;

s = .333313334133413

 (5)
> 

exact:=( 1/3 - 1/(n+1+1/(n+1)^3+1) );

6250249999999900000/18751875067501050009

 (6)
> 

evalf( [errK = K-exact, errA = A-exact, err_for=s-exact] );

[errK = 0., errA = 0.1e-14, err_for = 0.112e-12]

 (7)
> 

evalf[20]( [errK = K-exact, errA = A-exact, err_for=s-exact] );

[errK = -0.33461e-15, errA = 0.66539e-15, err_for = 0.11166539e-12]

 (8)
> 

 


Download KahanSum.mw

Bug corrected in inttrans[hilbert]...

vv 13143 Maple
inttrans hilbert
April 10 2018
2 2

There is a bug in inttrans:-hilbert:

> 

restart;
> 

inttrans:-hilbert(sin(a)*sin(t+b), t, s);
# should be:
sin(a)*cos(s+b);   expand(%);

sin(a)*cos(s)

  

sin(a)*cos(s+b)

  

sin(a)*cos(s)*cos(b)-sin(a)*sin(s)*sin(b)

 (1)
> 

########## correction ##############
> 

`inttrans/expandc` := proc(expr, t)
local xpr, j, econst, op1, op2;
      xpr := expr;      
      for j in indets(xpr,specfunc(`+`,exp)) do
          econst := select(type,op(j),('freeof')(t));
          if 0 < nops(econst) and econst <> 0 then
              xpr := subs(j = ('exp')(econst)*combine(j/('exp')(econst),exp),xpr)
          end if
      end do;
      for j in indets(xpr,{('cos')(linear(t)), ('sin')(linear(t))}) do
          if type(op(j),`+`) then
              op1:=select(has, op(j),t); ##
              op2:=op(j)-op1;            ##
              #op1 := op(1,op(j));
              #op2 := op(2,op(j));
              if op(0,j) = sin then
                  xpr := subs(j = cos(op2)*sin(op1)+sin(op2)*cos(op1),xpr)
              else
                  xpr := subs(j = cos(op1)*cos(op2)-sin(op1)*sin(op2),xpr)
              end if
          end if
      end do;
      return xpr
end proc:
> 

#######################################
> 

inttrans:-hilbert(sin(a)*sin(t+b), t, s); expand(%);

-(1/2)*cos(a-b)*sin(s)+(1/2)*sin(a-b)*cos(s)+(1/2)*cos(a+b)*sin(s)+(1/2)*sin(a+b)*cos(s)

  

sin(a)*cos(s)*cos(b)-sin(a)*sin(s)*sin(b)

 (2)
> 

 


Download hilbert.mw

 

Multiple integrals over sets defined by ...

vv 13143 Maple 2017
October 28 2017
6 11

I have just posted an article with this title at Maplesoft Application Center here.
It was motivated by a question posed by  Markiyan Hirnyk  here and a test problem proposed there by Kitonum.

Now I just want to give the promissed complete solution to Kitonum's test:

Compute the plane area of the region defined by the inequalities:

R := [ (x-4)^2+y^2 <= 25, x^2+(y-3)^2 >= 9, (x+sqrt(7))^2+y^2 >= 16 ];

plots:-inequal(R, x=-7..10, y=-6..6, scaling=constrained);

The used procedures (for details see the mentioned article):

ranges:=proc(simpledom::list(relation), X::list(name))
local rez:=NULL, r,z,k,r1,r2;
if nops(simpledom)<>2*nops(X) then error "Domain not simple!" fi;
for k to nops(X) do    r1,r2:=simpledom[2*k-1..2*k][]; z:=X[k];
  if   rhs(r1)=z and lhs(r2)=z then rez:=z=lhs(r1)..rhs(r2),rez; #a<z,z<b
  elif lhs(r1)=z and rhs(r2)=z then rez:=z=lhs(r2)..rhs(r1),rez  #z<b,a<z
  else error "Strange order in a simple domain" fi
od;
rez
end proc:

MultiIntPoly:=proc(f, rels::list(relation(ratpoly)), X::list(name))
local r,rez,sol,irr,wirr, rels1, w;
irr:=[indets(rels,{function,realcons^realcons})[]];
wirr:=[seq(w[i],i=1..nops(irr))];
rels1:=eval(rels, irr=~wirr);
sol:=SolveTools:-SemiAlgebraic(rels1,X,parameters=wirr):
sol:=remove(hastype, eval(sol,wirr=~irr), `=`); 
add(Int(f,ranges(r,X)),r=sol)
end proc:

MeasApp:=proc(rel::{set,list}(relation), Q::list(name='range'(realcons)), N::posint)
local r, n:=0, X, t, frel:=evalf(rel)[];
if indets(rel,name) <> indets(Q,name)  then error "Non matching variables" fi;
r:=[seq(rand(evalf(rhs(t))), t=Q)];
X:=[seq(lhs(t),t=Q)];
to N do
  if evalb(eval(`and`(frel), X=~r())) then n:=n+1 fi;
od;
evalf( n/N*mul((rhs-lhs)(rhs(t)),t=Q) );
end proc:

Problem's solution:

MultiIntPoly(1, R, [x,y]):  # Unfortunately it's slow; patience needed!
radnormal(simplify(value(%)));

evalf(%) = MeasApp(R, [x=-7..10,y=-6..6], 10000); # A rough numerical check
           61.16217534 = 59.91480000

 

Riemann hypothesis is false! (simple pro...

vv 13143 Maple
October 01 2017
2 2

# Riemann hypothesis is false! (simple proof)
 

restart;
assume( s>0, s<1/2, t>0 );
coulditbe(abs(Zeta(s+I*t))=0);
                              true

# Q.E.D.

Unfortunately coulditbe(Zeta(s+I*t)=0) returns FAIL, but our assertion is already demonstrated!

The moral: the assume facility deserves a much more careful implementation.

2 3 4 5 6 7 Page 4 of 7