# Kapitel 16.  Unendliche Reihen.
# ======================
# Aufgabe 16.1.
# Man verifiziere, da"s f"ur jede nat"urliche Zahl n gilt
#  sum(1/(n*(n+k)),n=1..infinity)=sum(1/j, j=1..k)/k;  .
# 
# Berechne die Differenz:
> Sum(1/(n*(n+k)),n=1..infinity)-Sum(1/j,
> j=1..k)/k=sum(1/(n*(n+k)),n=1..infinity)-sum(1/j, j=1..k)/k;
# Hier ist gamma die Eulersche Konstante,
# gamma=limit(sum(1/i,i=1..n)-ln(n),n=infinity) , siehe "?gamma", und
# Psi die sogenannte Digamma-Funktion, siehe "?Psi", eine der h"oheren
# transzendenten Funktionen.   
> simplify(");
# ==========================================
# Aufgabe 16.2
# ----------------
# Man bestimme den Grenzwert von sum(n*q^(n-1),n=1..infinity) .
> Sum(n*q^(n-1), n=1..infinity) =
> sum(n*q^(n-1), n=1..infinity);
# Der G"ultigkeitsbereich |q| < 1 ist vom Benutzer selbst zu
# bestimmen.
# ============================================
# Aufgabe 16.3.
# -----------------
# Man untersuche die folgenden Reihen auf Konvergenz:
# a) sum(k!^2/(2*k)!,k=1..infinity) .
> Sum(k!^2/(2*k)!, k=1..infinity) =
> sum(k!^2/(2*k)!, k=1..infinity);
# Maple verwendet die Gauss'sche hypergeometrische Funktion, siehe
# "?hypergeom".
> evalf(rhs("));
# ----------------------------------------------------------------------
# b) sum(binomial(3*k,k)/6^k, k=1..infinity) .
> Sum(binomial(3*k,k)/6^k,k=1..infinity)=sum(binomial(3*k,k)/6^k,k=1..in
> finity);
> evalf(rhs("));
# In der Tat ist nicht einmal die notwendige Bedingung f"ur Konvergenz
# erf"ullt:
> limit(binomial(3*k,k)/6^k,k=infinity);
# ----------------------------------------------------------------------
# c) sum(1/k^(3/2),k=1..infinity) .
> Sum(1/k^(3/2),k=1..infinity) =
> sum(1/k^(3/2),k=1..infinity);
# Hier ist die Reihe lediglich als Wert der Riemannschen 
# Zeta--Funktion an der Stelle 3/2 umgeschrieben worden. Unter
# "?zeta" findet man die Definition:
> zeta(z)=sum(1/k^z,k=1..infinity)#   f"ur Re(z)>1.
> evalf(rhs("));
# ----------------------------------------------------------------------
# d) sum(k^m*q^k,k=1..infinity)  f"ur nat"urliches m.
# Nebenbedingungen:
> assume(q, RealRange(-1,0) or
> RealRange(0,1));
> about(q);
> assume(m,integer);about(m);
# Wir erhalten kein Ergebnis f"ur den allgemeinen Fall:
> Sum(k^m*q^k, k= 1..infinity)=sum(k^m*q^k, k=1..infinity);
# Versuchsweise kann man einmal m=10 fixieren:
> Sum(k^10*q^k, k= 1..infinity)=sum(k^10*q^k,
> k= 1..infinity);
# Wir brauchen nicht den expliziten Grenzwert.
# Es h"atte gen"ugt, die Reihe sum(1/k^2,k=1..infinity) als Majorante
# heranzuziehen. Dazu h"atte man zeigen m"ussen, dass der folgende
# Grenzwert verschwindet, weil damit folgt: 
#  sum(k^m*q^k,k=1..infinity)<=sum(1/k^2,k=1..infinity) :
> limit(k^(m+2)*(q)^k, k=infinity);
> simplify(");
# "assume" !
# ----------------------------------------------------------------------
# e) sum((-1)^k*k/2^k,k=1..infinity .
> Sum((-1)^k*k/2^k, k=1..infinity) =
> sum((-1)^k*k/2^k, k=1..infinity);
# ----------------------------------------------------------------------
# f) sum(((k^2+(-1)^k)/(k^2+2))^(k^3),k=1..infinity) .
# Der folgende Befehl f"uhrt bei der zur Zeit verf"ugbaren Version von
# Maple V Release 4 zum Speicher"uberlauf. Selbst wenn das Programm von
# selbst anh"alt, ist beim anschlie"senden Versuch, das worksheet zu
# sichern, ein Absturz nicht auszuschliessen.
# Deshalb zeigen wir den Befehl nur in "auskommentierter" Form
# und empfehlen nicht, ihn auszuf"uhren.
# (Unter Release 3 wurde der Befehl nach kurzer Zeit unausgef"uhrt
# zur"uckgegeben.)
> #sum(((k^2+(-1)^k)/(k^2+2))^(k^3),
> k=1..infinity);
# Die Aufgabe ist keineswegs schwierig: Divergenz wegen
> limit(((k^2+(-1)^k)/(k^2+2))^(k^3),
> k=infinity);
# ----------------------------------------------------------------------
# g) sum((k+1)/(k^2+2*k+2),k=1..infinity) .
> Sum((k+1)/(k^2+2*k+2), k=0..infinity) =
> sum((k+1)/(k^2+2*k+2),k=0..infinity);
# ----------------------------------------------------------------------
# h) sum(k^(-m*k)*m^k f"ur nat"urliches m. 
> assume(m,posint);about(m);
> sum(k^(-m*k)*m^k, k=1..infinity);
# Da wir nur an der Konvergenzaussage interessiert sind, gen"ugt es, die
# Reihe sum(1/k^2,k=1..infinity)  als Majorantenreihe heranzuziehen: 
# Hinreichend f"ur die Beschr"anktheit des Quotienten 
# k^(-m*k)*m^k/k^2ist die Existenz des Grenzwerts
> limit(k^(-m*k)*m^k/k^2, k=infinity);
# ----------------------------------------------------------------------
# i) sum((2*k+1)/(4*k*(k+3)*(k+2)) .
> Sum((2*k+1)/(4*k*(k+3)*(k+2)), 
> k=1..infinity) =sum((2*k+1)/(4*k*(k+3)*(k+2)), k=1..infinity);
# ----------------------------------------------------------------------
# j)  sum((2+(-1)^k)/3^k,k=1..infinity) .
# Gefragt ist nach der Reihensumme. So geht es nicht:
> Sum((2+(-1)^k)/3^k, k=0..infinity) =
> sum((2+(-1)^k)/3^k, k=0..infinity);
# Da hier offensichtlich die Summe von zwei absolut konvergenten 
# (geometrischen) Reihen vorliegt, kann man diese auch einzeln zu
# berechnen versuchen, doch auch das f"uhrt nicht zum Ziel:
> Sum((2+(-1)^k)/3^k, k=0..infinity) =
> sum(2/3^k, k=0..infinity) +
> sum((-1)^k/3^k, k=0..infinity);
Error, (in sum/infinite) FAIL
# Offenbar hat das "sum}"-Programm eine eingebaute Sperre, die
# verhindert, dass die zweite geometrische Reihe ausgewertet wird. 
# Die Sperre kann man umgehen, indem man zuerst mit einem abstrakten q
# arbeitet und nachher 
# mittels "subs" f"ur q den Wert -1/3 einsetzt:
> Sum((2+(-1)^k)/3^k, k=0..infinity) =
> sum(2/3^k, k=0..infinity) + subs(q=-1/3, sum(q^k, k=0..infinity));
# ----------------------------------------------------------------------
# k) sum(1/(4*k^2-1),k=1..infinity) .
> Sum(1/(4*k^2-1), k=1..infinity) =
> sum(1/(4*k^2-1), k=1..infinity);
# ============================================
# ============================================
# Die folgende Aufgabe ist nicht mit Hilfe des Maple--Befehls 
# "sum" l"osbar.
# Aufgabe 16.4.
# -----------------
# Die Folge (a[n]) f"ur  n = 1,2, ... sei definiert durch: 
#  a[n]=0 , falls die Ziffer 5 in der Dezimaldarstellung von n 
#               auftritt,
#  a[n]=1/n , sonst.
# Man zeige: sum(a[n],n=1..infinity) ist konvergent.
# Dazu braucht man einen Befehl, der feststellt, ob in einer Dezimalzahl
# die Ziffer 5 vorkommt. Der Befehl "has" ist dazu nicht geeignet:
> has(567, '5');
> nops(567);
> op(1,567);
# Programm:
# "length(x)", "by", "trunc", ...
> hasfive:=proc(x) local j, rest, n;
> rest:=x; if type(x,integer)=false then FAIL else if
> length(x)>10 then RETURN(FAIL) else for j from length(x) by -1 to 1 do
> n:=trunc(rest/10^(j-1)); if  n=5 then RETURN('true') else
> rest:=rest-n*10^(j-1) fi od fi; RETURN('false');fi; end;
> hasfive(5);hasfive(5.5);hasfive(10^5);hasfive(5+10^5);hasfive(55555555
> 55);hasfive(55555555555);
# Damit kann man den Reihenterm a[k] berechnen:
> a:=proc(x) if hasfive(x)=true then 0 else 1/x fi end;
# Probe:
> seq(a(k), k=1..66);
# Es ist zuviel verlangt, die Konvergenz der Reihe von "sum"
# untersuchen zu lassen:
> sum('a(k)','k'=1..infinity);
# Aber man kann versuchen, die Summe bis 100, 1000, . . . auswerten
# zulassen. Dabei ist jedoch die langsame Konvergenz hinderlich.
> evalf(sum('a(k)','k'=1..100));
> evalf(sum('a(k)','k'=1..1000));
# F"ur die Summe bis 10000 ermitteln wir zun"achst die cpu - Zeit:
# (Vorsicht: Erhebliche Nachlaufzeit wegen evtl. Zugriffs auf den
# virtuellen Speicher. Siehe "Verbesserungsm"oglichkeiten". )
> #time(evalf(sum('a(k)','k'=1..10000)));

                                32.957

> #evalf(sum('a(k)','k'=1..10000));

                             7.865553848

# ============================================
# Verbesserungsm"oglichkeiten:
# ------------------------------------
# Benutze "evalhf" und "option remember".
# In "a" kann evalhf nicht eingebaut werden:
> b:=proc(x) if hasfive(x)=true then 0 else
> evalf(1/x) fi end;
> sum('b(k)','k'=1..1000);
Error, (in factor/factor) floats not handled
# Stattdessen verwende "evalhf" zum Aufsummieren der
# n-ten Teilsumme:
> TeilsummeH:=proc(n) local k, s; s:=0; if
> type(n,integer)=false then FAIL else for k from 1 to n do if
> hasfive(k)=false then s:=s + evalhf(1/k) fi od fi; s; end;
> time(TeilsummeH(10000));
Warning, computation interrupted
> TeilsummeH(10000);

                             7.865553834

# F"ur den letzten Befehl wird nochmals die gleiche Zeitspanne
# ben"otigt. Das kann verhindert werden durch Einf"ugen von "Option
# remember":
# Sicherheitshalber Neustart:
> restart;
> hasfive:=proc(x) local j, rest, n;
> rest:=x; if type(x,integer)=false then FAIL else if
> length(x)>10 then RETURN(FAIL) else for j from length(x) by -1 to 1 do
> n:=trunc(rest/10^(j-1)); if n=5 then RETURN('true') else
> rest:=rest-n*10^(j-1) fi od fi; RETURN('false');fi; end;
> TeilsummeRH:=proc(n) option remember; local k, s; s:=0; if
> type(n,integer)=false then FAIL else for k from 1 to n do if 
> hasfive(k)=false then s:=s+evalhf(1/k) fi od fi; s; end;
> time(TeilsummeRH(10000));

                                8.032

> TeilsummeRH(10000);

                             7.865553834

# F"ur TeilsummeRH(100 000) erh"alt man in etwa 96 Sekunden
# den Wert 9.262403941.
# ========================================
# Kommentar:
# Mit Papier und Stift geht es schneller:
# Teile die Summe in Abschnitte auf von k=10^n bis k=10^(n+1)-1
# f"ur n=0,1,2, ...
# Sch"atze in jedem Abschnitt die Anzahl der von Null verschiedenen
# Summanden und deren Maximalgr"osse nach oben ab. Erhalte 80 als obere
# Schranke f"ur die Summe.