Subject: numerical summation of series, rational case.
From: Michel Van den Bergh
Date: Tue, 26 Dec 2006 23:08:21 +0100
Hi,
Is what you are doing comparing with series like
1=sum(1/n^k-1/(n+1)^k,n,1,inf)=k*sum(1/n^(k+1)+...,n,1,inf)?
Then I understand your method. It is very nice!
I think it is indeed easy to write it purely in maxima code. Also I
don't see
why the magic number 20 is needed.
Regards,
Michel
>/* sumninf(a(n),p) computes an approximate value of
>sum(a(n),n,1,inf) with a(n) a rational function.
>
>It accelerate the convergence rate of a(n):
>there is b(n) such that lim a(n)/b(n) = 1 then
>
>a(n) - b(n) converges to zero more quickly than a(n) and
>sum(b(n),n,1,inf) known.
>
>Repeating this procedure we speed up the covergence rate,
> p is the number of times we speed up a(n).
>
>Finally sum(a(n),1,inf) = sum(a(n),1,20) + sum(a(n),20,inf)
> aprox
>sum(a(n),n,1,20)+ sum(b(n),n,20,inf)
>because r(n) = a(n) - b(n) converges very quickly to zero
> this procedure gives a good approximation. */
>
>--source code -------------------------------------------
>
>:lisp (defun principal(ex)
>(let* ((num (cadr ex)) (den (cddr ex)) (c1) (c2 (third den))
> (deg) (g2 (second den)))
>(if (atom num) (setq c1 num g1 0)
>(setq c1 (third num) g1 (second num)))
>(setq deg (- g2 g1)) (list '(MLIST SIMP) c1 c2 deg))) ;
>
>g1(ex):= block([l],l:?principal(rat(ex)), ex -
>l[1]/l[2]*product(1/(n+i),i,0,l[3]-1));
>
>matchdeclare(x_,integerp);
>
>defrule(uno,1/(n+x_),sum(-1/k,k,1,x_+20));
>
>
>sumninf(ex0,p):=block([ex],ex:rat(ex0),
> for i:1 thru p do ex:g1(ex),
> sum(ex0,n,1,20)+
> apply1(partfrac(ex0 - ex,n),uno));
>---------------------------------------------------
>
>
>