# Kapitel 20: Potenzreihen
# ================= 
# Eine der ersten Aufgaben f"ur Potenzreihen der Form
# sum(a[k]*z^k,k=0..infinity);  ist die Bestimmung
# des Konvergenzradius rho. Dazu werden als Standardmethoden die Formel
# von Cauchy-Hadamard
#  limit((abs(a[n])^(1/n),n=infinity)
# und die Quotientenformel 
#  rho=limit(a[k]/a[k+1],k=infinity) 
# benutzt.
# Weitere Aufgaben besch"aftigen sich mit der Bestimmung von Grenzwerten
# und der Berechnung von Cauchy-Produkten.
#  
# Aufgabe 20.1
# =========
# Man beweise f"ur komplexes z , |z|<1 :
#  sum(n^2*z^(n-1),n=1..infinity)=(z+1)/((1-z)^3 .
> restart;
> Sum(n^2*z^(n-1),n=1..infinity)=sum(n^2*z^(n-1),n=1..infinity);
> simplify(");
# Konvergenzradius nach Cauchy-Hadamard:
> 1/limit((n^2)^(1/n),n=infinity);
# Konvergenzradius nach Quotientenformel: 
> limit(n^2/(n+1)^2, n=infinity);
# Aufgabe 20.2
# =========
# Man zeige:
# i) sum(k!*z^k,k=1..infinity)  hat Konvergenzradius rho = 0 .
# ii) sum(z^k/product((1+1/j),j=1..k), k=1..infinity) hat
# Konvergenzradius rho = 1 .
> a:=k->k!;
> b:=k->1/product(1+1/j,j=1..k);
# Konvergenzradius nach Hadamard:
> 1/limit((a(k))^(1/k),k=infinity);
> 1/limit((b(k))^(1/k),k=infinity);
# Konvergenzradius mit Quotientenformel:
> limit(a(k)/a(k+1), k=infinity);
> limit(b(k)/b(k+1), k=infinity);
# Maple kann zur zweiten Reihe sogar den Grenzwert bestimmen:
> simplify(sum(b(k)*z^k,k=1..infinity));
# Aufgabe 20.3
# ==========
# Man bestimme die Konvergenzradien der folgenden Potenzreihen:
# i)  sum((k!/product(2*j+1,j=1..k))^2*z^k, k=1..infinity)   f"ur
# komplexes z,
# ii) sum(a^(k!)*z^k,k=1..infinity) f"ur a>0,
# iii) sum((-1)^k*z^k/2^(k+1),k=1..infinity) ,
# iv) sum((1+1/k)^(k^2)*z^k,k=1..infinity) .
# Zu i):
> a:=k->(k!/product(2*j+1,j=1..k))^2;
> 1/limit((a(k))^(1/k),k=infinity);limit(a(k)/a(k+1), k=infinity);
# Wenn nach dem Grenzwert gefragt wird, liefert Maple lediglich eine
# Umwandung in eine hypergeometrische Funktion:
> sum(a(k)*z^k,k=1..infinity);
# ----------------------------------------------------------------------
# -------------
# Zu ii):
> a:=k->alpha^(k!);
> 1/limit((abs(a(k)))^(1/k),k=infinity);limit(abs(a(k)/a(k+1)),k=infinit
> y);
> assume(alpha>0);
> 1/limit((abs(a(k)))^(1/k),k=infinity);limit(abs(a(k)/a(k+1)),k=infinit
> y);
> simplify("");
# Offenbar ist eine Fallunterscheidung n"otig:
> assume(0<alpha,alpha<1);
> about(alpha);
> 1/limit((abs(a(k)))^(1/k),k=infinity);limit(abs(a(k)/a(k+1)),k=infinit
> y);
> simplify("");
> alpha:=1;
> 1/limit((abs(a(k)))^(1/k),k=infinity);limit(abs(a(k)/a(k+1)),k=infinit
> y); alpha:='alpha':
> assume(1<alpha);
> 1/limit((abs(a(k)))^(1/k),k=infinity);limit(abs(a(k)/a(k+1)),k=infinit
> y); alpha:='alpha':
# ----------------------------------------------------------------------
# -------------
# Zu iii):
> a:=k->(-1)^k/2^(k+1);
> 1/limit((abs(a(k)))(1/k),k=infinity);limit(abs(a(k)/a(k+1)),k=infinity
> );
> simplify("");simplify("");
# Grenzwert:
> simplify(sum(a(k)*z^k,k=1..infinity));
# ----------------------------------------------------------------------
# ------------
# Zu iv):
> a:=k->(1+1/k)^(k^2);
> 1/limit((abs(a(k)))^(1/k),k=infinity);
> limit(abs(a(k)/a(k+1)),k=infinity);
> a:='a';
# ==============================================
# Cauchy-Produkte
# ---------------------
# Oft benutzt man den Satz "uber Cauchy-Produkte von Potenzreihen, um
# aus Reihen mit bekannten Grenzwerten eine neue Reihe zu bilden und
# deren Grenzwert zu bestimmen. Zwar kann die "sum"-Funktion f"ur 
# viele Reihen den Grenzwert direkt bestimmen, doch gibt es auch F"alle,
# in denen die Benutzung von Cauchy-Produkten unter Maple sinnvoll ist.
# Wir schreiben ein Programm, das aus zwei gegebenen Reihen S, T die
# Produktreihe bildet und deren Grenzwert als Produkt der Grenzwerte 
# von S und T ausgibt. Der gr"osste Teil des Programms besteht aus
# Abfragen, ob die Daten der Reihen zueinander passen, bzw. der
# Angleichung der Laufparameter. Wir beschr"anken uns auf den Fall, dass
# beide Reihen den Entwicklungspunkt Null haben und "uberlassen dem
# Leser weitere Verbesserungen.
# Das Programm soll mit dem Befehl "cp(S,T)" aufgerufen werden,
# wobei S und T zwei unausgewertete Reihen sind, also von der Form
> S:=Sum(a[j]*z^j,j=m..infinity);
> T:=Sum(b[i]*z^i,i=m..infinity);
# Dabei soll m eine nichtnegative ganze Zahl sein. Die Summationsindizes
# d"urfen verschieden sein, und der Fall a[i]=1 f"ur alle i ist
# zugelassen.
# Das Programm berechnet die Koeffizienten
> c:=n->sum(a[k]*b[n-k],k=m..n);
# des Cauchy-Produkts, die Grenzwerte der Reihen S und T, und gibt als
# Ergebnis aus:
#  sum(c(k)*z^k,k=m..infinity) = S*T .
> cp:=proc() local a, b, c, S, T, n, Anfang, KoeffS, KoeffT, stringS, 
> stringT; S:=args[1]; T:=args[2]; stringS:=`first argument should be an
> expansion in terms of powers of z only`; stringT:= `second argument
> should be an expansion in terms of powers of z only`; if nargs<>2then
> RETURN(`ERROR, must have two arguments`) fi; Anfang:=
> op([-1,2,1],S); if evalb(op(0,S)=Sum)=false or evalb(op(0,T)=Sum)=
> false then RETURN(`ERROR, both entries must be unevaluated Sums`)
> elif evalb(op([-1,1],S)=k)=false then S:=subs(op([-1,1],S)=k,S); fi; 
> if evalb(op([-1,1],T)=k)=false then T:=subs(op([-1,1],T)=k,T) fi; if 
> has(op(1,S),z^k)=false then ERROR(stringS) else KoeffS:=
> simplify(op(1,S)/z^k); if has(op(1,T),z^k)=false then ERROR(stringT)
> else KoeffT:=simplify(op(1,T)/z^k) fi; if
> evalb(op([-1,2,1],S)=op([-1,2,1],T))=false then RETURN(`ERROR,
> starting points must be equal in both sums`) elif
> evalb(op([-1,2,2],S)=infinity)=false or
> evalb(op([-1,2,2],T)=infinity)=false then RETURN(`ERROR, sums must
> beinfinite`) fi; a:=unapply(KoeffS,k); if has(a,z) then ERROR(stringS)
> fi;  
> b:=unapply(KoeffT,k); if has(b,z) then ERROR(stringT) fi; 
> c:=unapply(sum(a(k)*b(n-k), 'k'=Anfang..n),n); Sum(simplify(c(n))*z^n,
> 
> n=Anfang..infinity)=eval(subsop(0=sum,S))*eval(subsop(0=sum,T)); fi;
> end;
# Zun"achst testen wir das Programm. Kommentare folgen anschliessend.
# Ein einfaches Beispiel:
> S:=Sum(z^j,j=1..infinity);T:=Sum(m*z^m,m=1..infinity);
> cp11(S,T);
# Ein schwierigeres Beispiel:
> S:=Sum(z^k/k^2,k=1..infinity);
> T:=Sum(k*z^k,k=1..infinity);
> cp(S,T);
# Dabei ist die Funktion polylog(a,z) definiert durch 
#  polylog(a,z):=sum(z^n/n^a,n=1..infinity0 ,
# siehe auch "?polylog".
# Maple kann die Reihe auf der linken Seite nicht direkt berechnen:
> H:=z->sum((-n*Psi(1,n+1)-Psi(n+1)+n*Pi^2/6-gamma)*z^n,n=1..infinity);
# In der Tat ist sogar die Berechnung von einzelnen Werten wie etwa
# H(1/2) wegen zu grossen Speicherbedarfs unm"oglich. Um zu pr"ufen, ob
# unser Programm wirklich den Reihengrenzwert gefunden hat, verwenden
# wir deshalb die Teilsummen der Ordnung 100 von H. Definiere die m-te
# Teilsumme
> HH:=(m,z)->evalf(sum((-n*Psi(1,n+1)-Psi(n+1)+n*Pi^2/6-gamma)*z^n,n=1..
> m));
# Definiere die rechte Seite als eine Funktion "PL":
> PL:=z->evalf(polylog(2,z)*z/(z-1)^2);
# Nun vergleiche deren Werte:
> HH(100,1/2);
> PL(1/2);
# =============================================
# Kommentare zum Programm "cp".
# -----------------------------------------
#  Zun"achst zu dem hier verwendeten Befehl "op"
# und zur Z"ahlweise von Operanden in Maple:
# Der Befehl "op(k,expr)" veranlasst Maple, den k-ten
# Operanden des Ausdrucks "expr" auszugeben. Dazu muss man
# wissen, welches die Operanden von "expr" sind. Das kann
# man mittels "op" und "nops" erfahren. 
# Zum Beispiel hat der Ausdruck
> S:=Sum(a[j]*z^j,j=m..infinity);
# zwei Operanden:
> nops(S);
# und zwar:
> op(S);
# Zus"atzlich gibt es einen nullten Operanden: Er enth"alt in der
# Regel (siehe dazu "?op") den "surface datatype" oder 
# "basic type" des 
# Ausdrucks, hier ist das der Programm-Name:
> op(0,S);
# Jeder Operand ist wiederum aus Unteroperanden aufgebaut. Der zweite
# Operand von S etwa hat die Teile
> op(op(2,S));
# so dass man mit dem Befehl
> op(2,op(2,S));
# den Laufbereich und mit
> op(1,op(2,op(2,S)));
# den Summationsanfang herausholen kann. Dasselbe geht einfacher mit dem
# Befehl (man beachte die umgekehrte Reihenfolge):
> op([2,2,1],S);
# Je nach Kompliziertheit des Ausdrucks "expr" kann die 
# Stelle, an der ein bestimmter Operand (uns wird z. B. z^k
# interessieren) steht, verschieden sein. Deshalb ist es n"utzlich,
# gelegentlich 
# "von hinten" zu z"ahlen: -1, -2, ... stehen f"ur den letzten,
# vorletzten, ... (Unter-) Operanden:
> op([-1,-1,1],S);
# ============================================
# Zum Programmablauf:
# ---------------------------
# Zun"achst werden die beiden Variablen f"ur den internen Gebrauch
# in "S" und "T" umbenannt (falls sie nicht schon so heissen).
# Falls das Programm nicht mit zwei Argumenten
# aufgerufen worden ist, wird es mit "RETURN" abgebrochen,
# wobei der string `ERROR, must have two arguments` ausgegeben
# wird. Dann wird in die lokale Variable "Anfang" der untere
# Summationsindex der Summe S eingetragen und "uberpr"uft,
# ob beide Variable vom Typ "Sum" mit grossem S sind.
# Dabei haben wir den Befehl "evalb" verwendet, der zur Auswertung von
# Booleschen Ausdr"ucken dient: evalb(op(0,S)=Sum) liefert
# "true", wenn Sum der basic type von S ist, anderenfalls "false".
# Weiter wird nach dem Laufparameter der Summe S gesucht, und dieser
# wird ggfs. mit dem Befehl "subs" in "k" umgewandelt. 
# 
# Insgesamt ist dabei eine Schleife der Form
# 
# if--or--then--elif--then--fi
# 
# entstanden: Man beachte, dass am Ende nur ein einziges "fi" steht,
# obwohl in elif eigentlich ein zweites if enthalten ist.
# Dann wird auch in T der Laufparameter in "k" umgewandelt und mittels
# des "has"-Befehls "uberpr"uft, ob im ersten Operanden von S 
# "uberhaupt der Ausdruck z^k vorkommt. Wenn nicht, dann ist S keine
# Potenzreihe in z (m"oglicherweise ist S eine Entwicklung in Potenzen
# von (z-z[0]), und das Programm bricht mit der Fehlermeldung ab, die
# zuvor in "stringS" gespeichert worden war.
# Ein Beispiel hierzu:
> T:=Sum(z^j,j=1..infinity);S:=Sum(m*(z-1)^m,m=
> 1..infinity); 
> cp(S,T);
Error, (in cp) first argument should be an expansion in terms of
powers of z only
# Anderenfalls wird aus dem ersten Term von S der Faktor z^k
# herausdividiert und der Rest als "KoeffS" gespeichert. 
# Entsprechendes geschieht mit der Summe T. Nachdem gepr"uft ist,
# ob beide Summen den gleichen Startwert haben und bis Unendlich
# laufen, werden aus den Ausdr"ucken KoeffS und KoeffT mittels
# "unapply" Funktionen von k gemacht.
# 
# Nun k"onnten die Koeffizienten immer noch von z abh"angig sein,
# zum Beispiel, wenn zu dem zuletzt definierten T ein S von der 
# Form 
> S:=Sum(z^(k+1/2),k=1..infinity);
# gew"ahlt wird: Dann bricht das Programm ab mit der Meldung
> cp(S,T);
Error, (in cp) first argument should be an expansion in terms of
powers of z only
# Damit ist der Abfrageteil des Programms beendet, und das Berechnen 
# des Cauchy--Produkts geschieht mittels der beiden Befehle
# 
# c:=unapply(sum(a(k)*b(n-k), 'k'=Anfang..n),n);
# 
# Sum(simplify(c(n))*z^n,n=Anfang..infinity)=eval(subsop(0=sum,S))*eval(
# subsop(0=sum,T));
# Der erste bildet die Koeffizienten c[k] der Produktreihe aus der 
# "Schr"aglinienformel", und der zweite schreibt eine Gleichung,
# deren linke Seite die Produktreihe und deren rechte Seite das Produkt 
# der Grenzwerte der Einzelreihen ist. Dabei wurde der Befehl
# "subsop" verwendet, um in den Summen den jeweils nullten Operanden
# ("Sum") durch "sum" zu ersetzen. Durch eval(subsop(0=sum,S) wird 
# also das Ausrechnen der Summe S bewirkt. Damit ist das Programm
# beendet.
# ==============================================
# Ein n"utzliches Hilfsmittel
# --------------------------------
# zum "Uberpr"ufen des Programmablaufs ist die "trace"-Funktion,
# insbesondere bei der Fehlersuche.
# Man unterwirft das Programm "cp" mit dem Befehl
> trace(cp);
# der "Uberwachung. Wenn man das Programm dann aufruft, etwa f"ur 
> S:=Sum(z^j,j=1..infinity);T:=Sum(m*z^m,m=1..infinity);
# dann wird der interne Ablauf sichtbar. Hier sehen wir zum Beispiel, 
# dass die beiden Summen als Summen "uber k neu geschrieben wurden:
> cp(S,T);
# Die trace-Funktion bleibt so lange aktiv, bis sie mittels 
# "untrace" abgeschaltet wird.
> untrace(cp);
# ===============================================
# Hinweis
# ----------
# Die Maple--Befehle "powseries" und
# "powseries,multiply" sind f"ur die hier verfolgten Ziele
# nicht brauchbar. Auf sie werden wir in Kapitel 23 eingehen.
# ===============================================
# Eine Variante des Programms "cp":
# -------------------------------------------
# Man kann den "subsop"-Befehl am Ende durch "value" ersetzen, da 
# value, angewendet auf einen "inert"-Befehl, daruas einen sofort
# auf"uhrbaren Befehl macht:
> cp11:=proc() local a, b, c, S, T, n, Anfang, KoeffS, KoeffT, stringS, 
> stringT; S:=args[1]; T:=args[2]; stringS:=`first argument should be an
> expansion in terms of powers of z only`; stringT:= `second argument
> should be an expansion in terms of powers of z only`; if nargs<>2then
> RETURN(`ERROR, must have two arguments`) fi; Anfang:=
> op([-1,2,1],S); if evalb(op(0,S)=Sum)=false or evalb(op(0,T)=Sum)=
> false then RETURN(`ERROR, both entries must be unevaluated Sums`)
> elif evalb(op([-1,1],S)=k)=false then S:=subs(op([-1,1],S)=k,S); fi; 
> if evalb(op([-1,1],T)=k)=false then T:=subs(op([-1,1],T)=k,T) fi; if 
> has(op(1,S),z^k)=false then ERROR(stringS) else KoeffS:=
> simplify(op(1,S)/z^k); if has(op(1,T),z^k)=false then ERROR(stringT)
> else KoeffT:=simplify(op(1,T)/z^k) fi; if
> evalb(op([-1,2,1],S)=op([-1,2,1],T))=false then RETURN(`ERROR,
> starting points must be equal in both sums`) elif
> evalb(op([-1,2,2],S)=infinity)=false or
> evalb(op([-1,2,2],T)=infinity)=false then RETURN(`ERROR, sums must
> beinfinite`) fi; a:=unapply(KoeffS,k); if has(a,z) then ERROR(stringS)
> fi;  
> b:=unapply(KoeffT,k); if has(b,z) then ERROR(stringT) fi; 
> c:=unapply(sum(a(k)*b(n-k), 'k'=Anfang..n),n); Sum(simplify(c(n))*z^n,
> n=Anfang..infinity)=value(S)*value(T); fi; end;