# Lektion 9
# Integralrechnung
# © 2015  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(plots): with(plottools):
# Die Integration  (ein Ganzes ausmachend) ist die Umkehrung der Differenziation.  
# Differenzialrechnung:  y=f(x)      ->    y'=f '(x);
# Integralrechnung:        y'=f '(x)    ->     y=f (x);
# Im Zusammenhang mit der Integration nehmen wir noch folgende Umbenennung vor: f(x): Vorgegebene  1. Ableitung einer noch unbekannten Funktion. 
# Dann wird als Integral- oder Stammfunktion  F(x) einer gegebenen Funktion  f(x) eine Funktion bezeichnet, für die gilt:   F'(x)=f(x).;
# Bei der Integration einer Funktion  f(x)  wird eine  Stammfunktion  F(x) von  f(x) gesucht, deren Ableitung wieder die ursprüngliche Funktion  f(x) ergibt.;
# Schreibweise: F(x)+C  =  &int;f(x)&DifferentialD;x   mit  F'(x)=f(x). Die Funktion  f(x) wird als Integrand  bezeichnet und  x  heißt  Integrationsvariable. ;
# Zu jeder integrierbaren Funktion  f(x); gibt es unendlich viele Stammfunktionen  F(x) + C, die sich nur durch eine (additive)  Konstante  C  unterscheiden. 
# &int;f(x) &DifferentialD;x  =  &int;F'(x) &DifferentialD;x  =  &int;dF(x)=F(x)+C;
# Hinweis:  Alle Stammfunktionen  F(x)+C  besitzen an der festen Stelle  x =x[0]   dieselbe Steigung F'(x[0]). Sie gehen durch Parallelverschiebung ineinander über.;
# Ein Integral, dessen  Integrationskonstante  C  unbestimmt und damit nicht festgelegt ist, wird unbestimmtes Integral genannt. 
#  Beispiel: f(x);=2x   ’ F(x)+C;  =int(f(x), x);=  x^2; ,  da  F'(x) = f(x) = 2x  ist.  In Maple wird, ähnlich wie in vielen Tabellenbüchern der Ingenieurliteratur, die Konstante  C  nicht  ausgegeben.
> f:=2*x;
> Int(f,x)=int(f,x);
# Häufig ist es nicht möglich, die Stammfunktion einer Funktion  f(x)  in  geschlossener  Form anzugeben. Beispiele dazu sind:
> I1:=Int(exp(-x^2),x)=int(exp(-x^2),x):I2:=Int(sin(x)/x,x)=int(sin(x)/x,x):
> I3:=Int(cos(x)/x,x)=int(cos(x)/x,x):  I4:=Int(1/ln(x),x)=int(1/ln(x),x):
> print(I1,I2,I3,I4);
# Kann Maple keine geschlossene Form des Integrals berechnen, dann wird der Integralausdruck zurückgegeben:
> Int(exp(-x^2)*ln(x),x)=int(exp(-x^2)*ln(x),x);
# Das bestimmte Integral der Funktion  f(x)  lässt sich anschaulich deuten als Flächeninhalt A zwischen der Funktion f(x), der x-Achse und den zur y-Achse parallelen Geraden x = a und x = b. Das bestimmte Integral der Funktion  f(x)  wird symbolisch mit dem Zeichen   (&int;)[a]^(b)f(x) &DifferentialD;x   abgekürzt.;
> f:=x^2;
> Int(f,x=a..b)=int(f,x=a..b);
# Das inerte (untätige) Ojekt  Int  kann mit dem Befehl  value  berechnet werden:
> Int(f,x=a..b); value(Int(f,x=a..b));
> a:=1:b:=5:
> Int(f,x=a..b)=int(f,x=a..b);
> plot(f,x=a..b,title = "\nf(x) = x^2\n",titlefont = ["ROMAN", 15],color="NavyBlue",thickness=3,axes=boxed,gridlines,filled=[color="Blue",transparency =.5]);
# Zerlegt man die Fläche  A  in  n  Steifen der Breite  Delta*x = (b-a)/n;  und ersetzt  jeden Streifen durch ein Rechteck, dann ergibt sich für n = 5 folgendes Bild
> n:=5:
> plfn:=plot([seq(seq([a+(i+j)*(b-a)/n,(a+i*(b-a)/n)^2],j=0..1),i=0..n-1),[b,0]],color=red):
> plf:=plot(f,x=a..b,color=blue,thickness=2):
> plots[display]([plf,plfn],axes=frame,gridlines);
# Summieren wir nun über sämtliche Rechteckflächen, dann führt das auf die Untersumme
> i:='i':n:='n':
> USUM:=(b-a)/n*Sum((a+i*(b-a)/n)^2,i=0..n-1);
# Der Grenzübergang n ’ , der hier existiert, liefert 
> A:=simplify(limit(value(USUM),n = infinity));
# Die Untersumme konvergiert - allerdings sehr schlecht, wie die folgende Grafik zeigt - gegen einen Grenzwert, der geometrisch als Flächeninhalt  A  unter der Kurve  y = f (x)   gedeutet werden kann. Wir stellen den Fehler grafisch dar:
> s1:=seq(10*j,j=1..10):USUMJ:=subs(n=j*10,USUM): s2:=seq(A-value(USUMJ),j=1..10):
> p:=plot([s1],[s2],x=0..100,y=0..5,style=point,symbol=asterisk,color=blue):
> display(p,axes=frame,gridlines,labels=[n,Fehler]);
# Ist man am nummerischen Wert des Integrals interessiert, dann ist der Befehl
> Int(f,x=a..b)=int(f,x=a..b,numeric);
# zu verwenden. Ist eine der Grenzen des bestimmten Integrals eine Dezimalzahl, dann berechnet Maple das Integral nummerisch
> b:=5.;Int(f,x=a..b)=int(f,x=a..b);
# Um Maple anzuweisen, trotzdem eine symbolische Berechnung durchzuführen, ist der Parameter numeric = false zu setzen. 
> Int(f,x=a..b)=int(f,x=a..b,numeric = false);
# Ist die  obere oder untere  Grenze eines bestimmten Integrals variabel, so hängt der Integralwert noch von dieser (oberen oder unteren) Grenze ab und ist damit eine Funktion der (oberen oder unteren) Grenze. Um die (obere oder untere) Integrationsgrenze nicht mit der Integrationsvariablen zu verwechseln, wird die Integrationsvariable x formal in  t  umbenannt. 
> a:='a':f:=unapply(t^2,t);
> I1:=Int(f(t),t=a..x)=int(f(t),t=a..x):I2:=Int(f(t),t=x..a)=int(f(t),t=x..a):
> print(I1,I2);
> diff(I1,x);diff(I2,x);
# Maple kann auch uneigentliche Integrale lösen, das sind Integrale von folgender Form
# 1. Integrale mit unendlichem Integrationsintervall 
# 2. Integrale mit Unstetigkeitsstellen (Pole) im Integranden.
# Zu 1.: Man berechnet zunächst mit endlichen Grenzen und bildet dann den Grenzwert (das macht natürlich Maple für uns):
# int(f(x), x = a .. infinity) = limit(int(f(x), x = a .. b), b = infinity), int(f(x), x = -infinity .. b) = limit(int(f(x), x = -a .. b), a = infinity), int(f(x), x = -infinity .. infinity) = limit(int(f(x), x = -a .. a), a = infinity);
# Dazu ein Beispiel: f(x) = 1/(x^2+1); , (&int;)[0]^(infinity)f(x) &DifferentialD;x=(lim)(&int;)[0]^(b)f(x) &DifferentialD;x=(lim) arctan(b)=(Pi)/(2). ;
# Maple liefert korrekt:
> f:=1/(1+x^2);Int(f,x=0..infinity)=int(f,x=0..infinity);
# Zu 2.: Vorsicht ist geboten, wenn zur Berechnung des bestimmten Integrals über einen Pol hinwegintegriert werden soll. Wir zeigen das am Beispiel des Integranden  f(x) = 1/x^3; , der an der Stelle  x = 0  eine Polstelle besitzt.
> b:='b':f:=1/x^3;
> plot(f,x=-1..2,-10..10,title = "\nf(x) = 1/x^3\n",titlefont = ["ROMAN", 15],color="NavyBlue",thickness=2,filled=[color=blue,transparency =.5],axes=boxed,gridlines);
# Soll beispielsweise das uneigentliche  Integral   (int(f(x), x = -1 .. b))*berechnet*werden, dann*liefert*uns*Maple;
> Int(f,x=-1..b)=int(f,x=-1..b) assuming b::positive;
# In solchen Fällen ist wie folgt vorzugehen: Besitzt der Integrand an der Stelle x = u  eine Unstetigkeitsstelle (Pol), dann berechnet man zunächst die Integrale  links- und rechtsseitig davon und bildet dann, folgende Grenwerte, von denen man hofft, dass sie existieren.
# int(f(x), x = a .. b) = limit(int(f(x), x = a .. u-varepsilon), varepsilon = 0); + limit(int(f(x), x = u+varepsilon .. b), varepsilon = 0);
# Dazu ein Beispiel:   Wir integrieren f(x) zwischen x = -1 und x = b > 0:
#  (&int;)[-1]^(b)f(x) &DifferentialD;x= (lim)(&int;)[-1]^(-varepsilon)f(x) &DifferentialD;x +(lim) (&int;)[varepsilon]^(b)f(x) &DifferentialD;x = (lim) (1)/(2)(varepsilon^(2)-1)/(varepsilon^(2))+(lim) (1)/(2)(b^(2)-varepsilon^(2))/(varepsilon^(2)b^()) .;
> I1:=int(f,x=-1..-epsilon) assuming epsilon::positive;
> I2:=int(f,x=epsilon..b) assuming epsilon::positive, b::positive;
# Die Grenzwertbildung an den Einzelintegralen liefert: 
> limit(I1,epsilon=0); limit(I2,epsilon=0);
# Die Grenzwerte der einzelnen Integrale existieren leider nicht. Allerdings existiert der als Cauchy  Hauptwert  (CH) bezeichnete Grenzwert
# CH int(f(x), x = -1 .. b) = limit([int(f(x), x = -1 .. -varepsilon)+int(f(x), x = varepsilon .. b)], varepsilon = 0) and limit([int(f(x), x = -1 .. -varepsilon)+int(f(x), x = varepsilon .. b)], varepsilon = 0) = limit([(1/2)*(varepsilon^2-1)/varepsilon^2+(1/2)*(b^2-varepsilon^2)/(varepsilon^2*b^2)], varepsilon = 0) and limit([(1/2)*(varepsilon^2-1)/varepsilon^2+(1/2)*(b^2-varepsilon^2)/(varepsilon^2*b^2)], varepsilon = 0) = (1/2)*(b^2-1)/b^2;
# Der Cauchy Hauptwert existiert und berechnet sich zu
> limit(I1+I2,epsilon=0);
# Dieses Ergebnis liefert uns auch Maple mit der Zusatzinformation  CauchyPrincipalValue  
> int(f,x=-1..b,CauchyPrincipalValue) assuming b::positive;
# Wir betrachten einen weiteren Integranden, 
> f:=unapply(1/(x+a)^2,x,a);
# der an der Stelle x = -a eine Polstelle besitzt. Integrieren wir von x = 0  bis  x = 2, dann macht  Maple darauf aufmerksam, dass im Falle -2 < a < 0 eine Diskontinuität im Integranden vorliegt, und damit der Wert des Integrals  wird.  
> Int(f(x,a),x = 0..2)=int(f(x,a),x = 0..2);
# Mit der Zusatzoption continuous teilen wir Maple mit, dass innerhalb des gewählten Integrationsbereiches der Integrand stetig verläuft. 
> Int(1/(x+a)^2,x = 0 .. 2)=int(1/(x+a)^2,x = 0 .. 2,'continuous');
# Wir betrachten das uneigentliche Integral   Typesetting[delayDotProduct]((int(ln*x/(x^2+1), x = 0 .. infinity)) . Der, Graf, true)*des*Integranden, (der*an*der)*Stelle*x = verschwindet, zeigt*F*olgendes;
> f:=ln(x)/(1+x^2);
> plot(f,x=0..5,y=-2..0.5,color="NavyBlue",thickness=2,filled=[color="Blue",transparency =.5],gridlines);
> Int(f,x=0..1)=int(f,x=0..1); Int(f,x=1..infinity)=int(f,x=1..infinity);
# Die Catalansche Konstante  Catalan hat den nummerischen Wert
> evalf(Catalan);
# Damit folgt
> Int(f,x=0..infinity)=int(f,x=0..infinity);
# Mittels der Heaviside-Funktion lassen sich Integranden abschnittsweise definieren. Maple kann damit rechnen:
> Int(g(x)*Heaviside(x-1),x =-1..3)=int(g(x)*Heaviside(x-1),x =-1..3);
# Als weiteres Beispiel eines uneigentlichen Parameterintegrals betrachten wir die Laplace-Transformierte der Funktion f(t);. Die Laplace-Transformation ist in der Schwingungslehre von großer Bedeutung.
> f:='f': t:='t':
> p:=unapply(Int(f(t)*exp(-s*t),t=0..infinity),s);
# Wir berechnen die Laplace-Transformierte des Einheitssprungs an der Stelle  t = 2
> f(t):=Heaviside(t-2);
> plot(f(t),t=0..5,axes=boxed,color=blue,thickness=2,gridlines);
> value(p(s));
# Bequemer erfolgt die Transformation mit dem Befehl  laplace  aus der Bibliothek inttrans: 
> with(inttrans):
> LT:=laplace(f(t),t,s);
# Die Umkehrung der Laplace-Transformation in den Originalbereich erfolgt durch den Befehl  invlaplace  
> invlaplace(LT,s,t);
# Ein weiteres Beispiel
> f:=piecewise(t>0 and t<3,1,t>3,0);
> f:=convert(f,Heaviside);
> p(s)=value(%);
> laplace(f,t,s);
# Wie bereits erwähnt, lassen sich nicht zu allen Funktionen f(x) die Stammfunktionen in geschlossener Form angeben. Hierzu als weiteres Beispiel die von  Euler mittels eines uneigentlichen Parameterintegrals  definierte  Gamma-Funktion: 
> Int(t^(x-1)*exp(-t),t=0..infinity)=int(t^(x-1)*exp(-t),t=0..infinity);
# Die  Gamma-Funktion besitzt die Eigenschaft, dass die Werte von  GAMMA(x+1); für ganzzahlige  x  gleich  x!  (x-Fakultät) sind. Ein Zahlenbeispiel bestätigt das
> GAMMA(3)-2!;
> plot(GAMMA(x),x=-5..5,-5..5,axes=boxed,color=blue,gridlines);
# Ein weiteres Beispiel:
> f:=unapply(x^(-1/2)*cos(x),x);
# Wir stellen die Funktion  f(x)  grafisch dar.;
> plot(f(x),x=0..5*Pi,y=-1..2,color="NavyBlue",thickness=2,filled=[color="Blue",transparency =.5],axes=framed,gridlines);
# Maple liefert uns folgendes Integral
> Int(f(t),t=0..x)=int(f(t),t=0..x);
# Das FresnelCosinus-Integral ist in Maple wie folgt definiert:  FresnelC(x)=(&int;)[0]^(x)cos(1/2 Pi t^( 2)) &DifferentialD;t.   ;
> INTA:=unapply(int(f(t),t=0..x),x);
> plot(INTA(x),x=0..2*Pi,axes=boxed,gridlines,thickness=2,color=blue);
# Wir werten das Integral an der Stelle  x = 5/4À  nummerisch aus und erhalten
> evalf(INTA(5/4*Pi));
# Entwickeln wir den Integranden  f(t);  an der Stelle t = 0  in eine Potenzreihe 
> ser:=convert(series(f(t),t=0,10),polynom);
# und integrieren anschließend gliedweise, was im betrachteten Intervall  absolute  Konvergenz voraussetzt, dann liefert Maple
> INTR:=unapply(int(ser,t=0..x),x);
# Ein Vergleich mit der (analytisch) berechneten Stammfunktion ergibt  an der Stelle x = 5/4À
> evalf(INTR(5/4*Pi));
# Wir vergleichen grafisch die Näherung der Stammfunktion mit der nummerisch ermittelten Kurve
> pl1:=plot(INTA(x),x=0..2*Pi,y=0..2,thickness=2,color=blue,thickness=2):
> pl2:=plot(INTR(x),x=0..2*Pi,y=0..2,thickness=1,color=red,thickness =1):
> display([pl1,pl2],gridlines,axes=boxed);
# Wir lassen uns den wahren Fehler grafisch ausgeben, der für x < Pi;  immerhin kleiner als 0,42% ist
> plot(abs(INTA(x)-INTR(x)),x=0..Pi,thickness=2,color=blue,gridlines,axes=boxed);
# Auch Vektor- und Matrizenfunktionen können mit Maple integriert werden. Das funktioniert wieder mit dem Befehl map. Wir erzeugen uns zunächst einen Vektor und eine Matrix, deren Komponenten von x und y abhängen:
> j:='j':i:='i':
> f:= unapply(y^j*x^(j-1),j):  g:= unapply(y^(i)*x^(i-j+1),i,j):
> A:=Vector(3,f):B:=Matrix(3,3,g): print(A,B);
> AI:=map(Int,A,y) = map(int,A,y);
> BI:=map(Int,B,x) = map(int,B,x);
# Mit Maple lassen sich auch  Mehrfachintegrale  berechnen. Wir betrachten die folgenden Doppelintegrale
> Int(x*y^2,[x,y])=int(x*y^2,[x,y]);
> Int(x*y^2,[y,x])=int(x*y^2,[y,x]);
> Int(Int(x*y^2,y=0..sqrt(1-x)),x=0..1)= int(int(x*y^2,y=0..sqrt(1-x)),x=0..1);
> Int(x*y^2,[y=0..sqrt(1-x),x=0..1])= int(x*y^2,[y=0..sqrt(1-x),x=0..1]);
# Enthalten die Integrationsgrenzen Dezimalzahlen, dann wertet Maple das Integral nummerisch aus:
> Int(Int(x*y^2,y=0..sqrt(1-x)),x=0..1.)= int(int(x*y^2,y=0..sqrt(1-x)),x=0..1.);
# Die gemeine Schraublinie (Helix) hat die Parameterdarstellung x = a cos t,  y = a sin t, z = p t. 
> r:='r': t:='t': a:='a': p:='p': n:='n':
> r:=Vector([a*cos(t),a*sin(t),p*t]):rp:=map(diff,r,t):
> print(r,rp);
# Die Bogenlänge dieser Raumkurve berechnet sich zu:  L(t) = int(sqrt((diff(r(tau), tau))^2), tau = 0 .. t);     NULL;
> rpq:=simplify(DotProduct(rp,rp,conjugate=false),trig);
> L:=int(sqrt(rpq),tau=0..t);
# Ist h die Ganghöhe und n die Windungszahl, dann gilt p = h/(2À)  und  t = 2 n À 
> simplify(subs(p=h/2/Pi,t=2*n*Pi,L));
# Der Durchhang eines schweren  Seils wird beschrieben durch die Kettenlinie  y = a cosh(x/a)   a > 0
> y:=unapply(a*cosh(x/a),x,a);
> plot(y(x,1),x=-1..1,gridlines,axes=boxed,color="NavyBlue",thickness=2);
# Die Neigung der Tangente  an die  Kettenlinie ist 
> ystr:=unapply(diff(y(x,a),x),x,a);
# Die Seillänge L im Interval [0...x] ermitteln wir aus der Beziehung  L(x)=(&int;)[0]^(x)sqrt(1+y^('2)(t)) &DifferentialD;t; .
> L:=unapply(int(sqrt(1+ystr(t,a)^2),t=0..x),x,a) assuming x::positive;
> simplify(L(x,a)) assuming x::positive, a::positive;
# Mit Maple lassen sich auch Ableitung von Parameterintegralen vornehmen:
> restart:
# Das von einem Parameter (hier x) abhängige Integral 
> g:=int(f(x,u),u=a(x)..b(x));
# wird  als Parameterintegral bezeichnet. Der Parameter x  tritt  im Integranden und an der oberen und unteren  Grenze des Integrals auf. Die Ableitung eines solchen Integrals nach dem Parameter  x  erfolgt mittels der Leibnizregel. Maple liefert korrekt:
> Diff(g,x)=diff(g,x);
# Hinweis: Solche Integrale treten beispielsweise im Zusammenhang mit der Berechnung der Knicklasten von Stützen mit veränderlichem Querschnitt unter Eigengewichtsbelastung (spezifisches Gewicht ³) auf  (s.h. Technische Mechanik 2, Festigkeitslehre mit Maple-Anwendungen, Kap. 10).
> h:=hu+(ho-hu)/l*x;  #linear veränderliche Querschnittshöhe über die Stützenlänge l
;
> A:=unapply(b*h,x);  #Querschnittsfläche
;
> int(gamma*A(u)*(w(x)-w(u)),u=x..l); 
> simplify(diff(%,x));
> 
;
