#                                          MAPLE  für ANALYSIS II
     #                                          =====================
     # LEHRSTUHL A FÜR MATHEMATIK                                                      [15.txt]
     #
     # Prof. Dr. E. Görlich
     #
     #                                    Sommersemester 2000
     #
     #                                    18.05.2000
     #
     #
     #
     # 14) Uneigentliche Integrale
     # (limit, unapply, middlesum, verallg. Fresnel-Integral, changevar,
     # Clausen-Integral)
     # ==================================================================
     #
     # Zu Uebung 5, Aufgabe 1:
     # ===================
     # a)
     > Int(1/(5-x)^3,x=1..7)=int(1/(5-x)^3,x=1..7);
     > Limit(Int(1/(5-x)^3,x=1..X),X=5,
     > left)=limit(int(1/(5-x)^3,x=1..X),X=5, left);
     # -> divergent
     # =============================================================
     # b)
     > Int(log(1/(1-x)),x=0..1)=int(log(1/(1-x)),x=0..1);
     # -> konvergent
     # =============================================================
     # c)
     > Int(cos(1/x)/x^(8/7),x=0..2)=int(cos(1/x)/x^(8/7),x=0..2);
     # Ist diese Antwort richtig? Definiere eine Funktion von epsilon;
     # mittels "unapply":
     > test1:=unapply(int(cos(1/x)/x^(8/7),x=epsilon..2),epsilon);
     # Der Grenzwert limit(test1(epsilon),epsilon = 0,right); ist gesucht.
     > plot(test1(epsilon),epsilon=0.1..1, VIEW=-1000..1000);
     > kernelopts(cputime);
     > evalf(test1(0.1));
     > evalf(test1(0.11));
     > evalf(test1(0.09));
     > evalf(test1(0.08));
     > evalf(test1(0.07));
     > evalf(test1(0.06));
     > evalf(test1(0.05));
     > seq(evalf(test1(0.01+j*0.01)),j=1..10);
     > plot(test1(epsilon),epsilon=0.01..0.1);
     # Es hat den Anschein, dass limit(Int(cos(1/x)/(x^(8/7)),x = epsilon ..
     # 2),epsilon = 0,right); =0.1 ist, ungefaehr.
     # Ein Versuch mit Riemann-Summen zeigt, dass obige Divergenzbehauptung
     # sich auf die Divergenz des *gewoehnlichen*
     # Riemann-Integrals bezieht:
     > kernelopts(cputime);
     > %/60*minuten;
     > with(student);
     > middlebox(cos(1/x)/x^(8/7),x=0..2,10);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,10);
     > evalf(%);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,100):evalf(%);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,200):evalf(%);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,300):evalf(%);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,400):evalf(%);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,500):evalf(%);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,600):evalf(%);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,700):evalf(%);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,800):evalf(%);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,900):evalf(%);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,1000):evalf(%);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,1100):evalf(%);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,1200):evalf(%);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,1300):evalf(%);
     > middlesum(cos(1/x)/x^(8/7),x=0..2,1400):evalf(%);
     > rightsum(cos(1/x)/x^(8/7),x=0..2,1400):evalf(%);
     > leftsum(cos(1/x)/x^(8/7),x=0..2,1400):evalf(%);
     > evalf(1/700*Sum(700*cos(700*1/i)*700^(1/7)/(i^(8/7)),i = 1 .. 1399));
     # ->Dies deutet an, dass das gewoehnliche Riemann-Integral, also der
     # Grenzwert der Riemann-Summe
     # middlesum(cos(1/x)/x^(8/7),x=0..2,n) fuer n->infinity; , in der Tat
     # nicht existiert.
     # Untersuche nochmals das uneigentliche Integral:
     # ======================================================================
     # ===
     > restart:
     > with(student):
     > ?changevar
     > changevar(1/x=u,Int(cos(1/x)/x^(8/7),x=0..2));
     > evalf(%);
     # Dieser Hinweis bezieht sich wiederum auf das gewoehnliche R-Integral.
     # Nun betrachte das Vorzeichen des Interanden:
     > plot(cos(x)/x^(6/7),x=1/2..15);
     # und teile das Integral in Teilstuecke auf:
     > k:=1:Int(cos(u)/(u^(6/7)),
     > u=Pi/2+2*k*Pi..Pi/2+(2*k+1)*Pi)=evalf(int(cos(u)/(u^(6/7)),
     > u=Pi/2+2*k*Pi..Pi/2+(2*k+1)*Pi));Int(cos(u)/(u^(6/7)),
     > u=Pi/2+(2*k+1)*Pi..Pi/2+(2*k+2)*Pi)=evalf(int(cos(u)/(u^(6/7)),
     > u=Pi/2+(2*k+1)*Pi..Pi/2+(2*k+2)*Pi));
     > k:=2:Int(cos(u)/(u^(6/7)),
     > u=Pi/2+2*k*Pi..Pi/2+(2*k+1)*Pi)=evalf(int(cos(u)/(u^(6/7)),
     > u=Pi/2+2*k*Pi..Pi/2+(2*k+1)*Pi));Int(cos(u)/(u^(6/7)),
     > u=Pi/2+(2*k+1)*Pi..Pi/2+(2*k+2)*Pi)=evalf(int(cos(u)/(u^(6/7)),
     > u=Pi/2+(2*k+1)*Pi..Pi/2+(2*k+2)*Pi));
     > k:=3:Int(cos(u)/(u^(6/7)),
     > u=Pi/2+2*k*Pi..Pi/2+(2*k+1)*Pi)=evalf(int(cos(u)/(u^(6/7)),
     > u=Pi/2+2*k*Pi..Pi/2+(2*k+1)*Pi));Int(cos(u)/(u^(6/7)),
     > u=Pi/2+(2*k+1)*Pi..Pi/2+(2*k+2)*Pi)=evalf(int(cos(u)/(u^(6/7)),
     > u=Pi/2+(2*k+1)*Pi..Pi/2+(2*k+2)*Pi));
     # -> Das Integral sollte sich nach oben und unten abschaetzen lassen
     # durch alternierende Reihen mit
     # monoton fallenden Koeffizienten, also konvergieren aufgrund des
     # Leibniz-Kriteriums (Satz
     # IV.1.10 der Vorlesung).
     # Berechne int(cos(u)/(u^(6/7)),u = 1/2 .. k*Pi);  fuer k= 5, 50, 60,
     # 70, 80, 100, 200, 300:
     > k:='k';
     > test2:=k->evalf(int(cos(u)/(u^(6/7)),u = 1/2 .. k*Pi));
     > test2(5);
     > test2(50);
     > test2(60);
     > test2(70);
     > test2(80);
     > test2(100);
     > test2(200);
     > test2(300);
     > limit(test2(k),k=infinity);
     > evalf(%);
     > kernelopts(cputime);
     # -->> Maple hat berechnet, dass das uneigentliche Integral existiert,
     # aber keinen Hinweis darauf gegeben,
     # in welchen mathematischen Zusammenhang es gehoert. Ein Blick in eine
     # Formelsammlung
     # (Abramowitz-Stegun, Handbook of Mathematical Tables, Dover 1968, S.
     # 262, # 6.5.7)
     # zeigt, dass es sich hier um das "verallgemeinerte Fresnel-Integral"
     # handelt:
     > C(x,a) = Int(t^(a-1)*cos(t),t = x .. infinity);#
     # Es konvergiert fuer alle komplexen a mit Re(a) < 1;. In unserem Fall
     # ist a=1/7.
     # Die Hilfe zum Stichwort "Fresnel" bringt jedoch nur Material zu den
     # gewoehnlichen Fresnel-Integralen:
     #
     > ?Fresnel
     # ======================================================================
     # ========
     # Zu Uebung 5, Aufgabe 1 d):
     > Int(1/sqrt(x-x^2),x=0..1)=int(1/sqrt(x-x^2),x=0..1);
     # Fuehre die Rechnung im Einzelnen durch:
     # Nach dem Skript "Integrationsmethoden", Abschnitt 8 A, bietet sich die
     # Substitution
     # t=sqrt(x-x^2) an:
     > changevar(t=sqrt(x-x^2),Int(1/(sqrt(x-x^2)),x = 0 .. 1), t);
     # Die Integralgrenzen zeigen, dass hier ein Fehler vorliegt: Maple
     # hindert den Benutzer nicht
     # daran, eine Substitutionsfunktion einzusetzen, die nicht durchweg
     # monoton wachsend ist.
     # Offenbar ist der Integrand im Ausgangsintegral symmetrisch um den
     # Punkt x=1/2, d.h.
     # invariant gegenueber der Substitution x=1-u:
     > subs(x=1-u,1/(sqrt(x-x^2)));
     > normal(%);
     # Also ist das Ausgangsintegral gleich
     > 2*Int(1/sqrt(x-x^2),x=0..1/2);
     # und hier darf t=sqrt(x-x^2) substituiert werden:
     > changevar(t=sqrt(x-x^2),2*Int(1/sqrt(x-x^2),x=0..1/2),t);
     # Nun substituiere noch u=2*t:
     > changevar(u=2*t, 2*Int(2*1/(sqrt(1-4*t^2)),t = 0 .. 1/4*sqrt(4)),u);
     # Nach dem Skript, "Integrationsmethoden", Abschnitt 8 F, sollte nun
     # u=cos(t) substituiert werden:
     > changevar(u=cos(t), 2*Int(1/(sqrt(1-u^2)),u = 0 .. 1), t);
     # (Man beachte, dass die obere Grenze nicht 1/(2*Pi) sondern Pi/2
     # bedeutet). Ergebnis:
     > 2*int(1,t=0..Pi/2);
     # Das ist die oben erhaltene Loesung.
     # ================================================================
     # Zu Uebung 5, Aufgabe 1 e):
     > Int(x/sqrt(exp(2*x)-1),x=0..5)=int(x/sqrt(exp(2*x)-1),x=0..5);
     # Hier ist _k1 eine von Maple eingefuehrte Variable.
     > evalf(%);
     # Das bedeutet insbesondere, dass das Integral existiert. Auch hier
     # wollen wir versuchen, die Maple-Loesung in
     # Einzelschritten nachzuvollziehen:
     > with(student):
     > changevar(x=ln(1/sin(t)),Int(x/sqrt(exp(2*x)-1),x=0..5),t);
     > ?simplify
     > simplify(%,real);
     # Offenbar kann der Integrand vereinfacht werden zu ln(1/(sin(t)).
     > int(ln(1/sin(t)),t=arcsin(1/exp(5))..Pi/2);
     # Probe:
     > evalf(%);
     # Der Realteil stimmt mit dem anfangs erhaltenen numerischen Wert
     # ueberein.
     # Das fuehrt uns auf das Stichwort:
     > ?dilog
     # Im Zusammenhang mit dem Dilogarithmus findet man in Abramowitz-Stegun
     # S. 1005, # 27.8.1 das Clausensche Integral :
     # Fuer 0 < theta;, theta <= Pi; gilt:
     > -Int(ln(2*sin(t/2)),t = 0 .. theta) = Sum(sin(k*theta)/(k^2),k = 1 ..
     > infinity);#
     # Daraus ergibt sich fuer das gesuchte Integral die Darstellung:
     # ============================================================
     > (1/2)*(sum(sin(k*Pi)/k^2,k=1..infinity)-sum(sin(k*2*arcsin(exp(-5)))/k
     > ^2,k=1..infinity)-Pi*ln(1/2)+2*arcsin(exp(-5))*ln(1/2));
     > evalf(%);
     # Ersetze die unendliche Reihe durch die Summe bis k=1000, 3000, 8000:
     > evalf((1/2)*(sum(sin(k*Pi)/k^2,k=1..1000)-sum(sin(k*2*arcsin(exp(-5)))
     > /k^2,k=1..1000)-Pi*ln(1/2)+2*arcsin(exp(-5))*ln(1/2)));
     > evalf((1/2)*(sum(sin(k*Pi)/k^2,k=1..3000)-sum(sin(k*2*arcsin(exp(-5)))
     > /k^2,k=1..3000)-Pi*ln(1/2)+2*arcsin(exp(-5))*ln(1/2)));
     > evalf((1/2)*(sum(sin(k*Pi)/k^2,k=1..8000)-sum(sin(k*2*arcsin(exp(-5)))
     > /k^2,k=1..8000)-Pi*ln(1/2)+2*arcsin(exp(-5))*ln(1/2)));
     # -> Gute Uebereinstimmung mit dem anfangs erhaltenen Wert.
     # =======================================================