Zurück zum Inhaltsverzeichnis

Kapitel 17

Das Integral

Die Integralrechnung ist vielleicht eine der wichtigsten Anwendungen für ein CAS: Statt in endlosen Tabellen nachzuschlagen und sich mit Integrationsregeln herumzuplagen, gibt man einen Befehl ein.

Was in der Forschung höchst willkommen ist hat aber auch Auswirkungen auf die Lehre. Inwieweit müssen Schüler und Studenten von morgen noch die Integration beherrschen? Neben der Fähigkeit eines CAS, Integrale symbolisch zu berechnen, spielt aber - einmal mehr - die Grafik eine wichtige Rolle und man könnte sich die Einführung des Integrals in Physik etwa so vorstellen:

Im Wechselstromkreis sind Spannung und Strom phasenverschoben. Welche Arbeit wird verrichtet?

> restart: with(plottools):

> P:='sin(t)*cos(t-Pi/5)'; # unausgewertet, damit Maple nicht vereinfacht

P := sin(t)*cos(t-1/5*Pi)

Das Produkt aus Spannung und Strom stellt den Momentanwert der Leistung P dar. Für genügend kurze Zeitintervalle Delta*t kann man die Arbeit durch Delta*W = P*Delta*t annähern. Das würde dann so aussehen:

> student[leftbox](P, t=0..5, 30, shading=yellow);

[Maple Plot]

Und die Arbeit während einer größeren Zeitspanne müsste dann ungefähr die Summe der einzelnen Beiträge sein:

> W:=student[leftsum](P, t=0..5, 30);

W := 1/6*Sum(sin(1/6*i)*sin(1/6*i+3/10*Pi),i = 0 .....

> evalf(%);

1.890470645

Natürlich fangen dann alle Schüler sofort an, die Zerlegung zu verfeinern. Einige ändern auch die Grenzen und wundern sich über negative Ergebnisse und die Experimentierfreudigsten ändern die Phasenverschiebung von Spannung und Strom und strahlen, wenn die Arbeit 0 herauskommt.

Dann wird es Zeit, vom Einzelfall auf einen allgemein gültigen Zusammenhang zu schließen.

>

> Wn:=student[leftsum](P, t=0..5, n);

Wn := 5*Sum(sin(5*i/n)*sin(5*i/n+3/10*Pi),i = 0 .. ...

value(%); # mit diesem Befehl kann man Verwirrung stiften

Kann Maple einen Term für die Arbeit finden, wenn man die Zerlegung beliebig fein macht?

> limit(Wn,n=infinity);

1/4*sin(3/10*Pi)+5/2*cos(3/10*Pi)-1/4*sin(10+3/10*P...

Ja! Und dazu gehört auch eine Zahl:

> evalf(%);

1.921365057

>

Und weil man solche Berechnungen oft benötigt, schreibt man kurz:

> Int(P, t=0..5)=int(P, t=0..5);

Int(sin(t)*sin(t+3/10*Pi),t = 0 .. 5) = 5/2*cos(3/1...

> combine(%);

Int(1/2*cos(3/10*Pi)-1/2*cos(2*t+3/10*Pi),t = 0 .. ...

> evalf(%);

1.921365057 = 1.921365057

>

Anmerkung

Wenn man die Berechnung für volle Perioden durchführt, kommt man schon mit sehr groben Zerlegungen aus. Warum?

> Wp4:=student[leftsum](P, t=0..2*Pi, 4);

Wp4 := 1/2*Pi*Sum(sin(1/2*i*Pi)*sin(1/2*i*Pi+3/10*P...

> value(%);

Pi*sin(1/5*Pi)

> evalf(%);

1.846581831

> Wpn:=student[leftsum](P, t=0..2*Pi, n);

Wpn := 2*Pi*Sum(sin(2*i*Pi/n)*sin(2*i*Pi/n+3/10*Pi)...

> simplify(value(%));

cos(3/10*Pi)*Pi

> evalf(%);

1.846581830

> evalf(value(Wp4)-simplify(value(Wpn)));

.1e-8

> eval(int(P,t=0..tb),tb=2*Pi);

cos(3/10*Pi)*Pi

>

Ohne die Angabe von Grenzen erhält man eine Stammfunktion

> Int(P,t) = int(P, t);

Int(sin(t)*sin(t+3/10*Pi),t) = 1/2*cos(3/10*Pi)*t-1...

> W:=int(P, t);

W := 1/2*cos(3/10*Pi)*t-1/4*sin(2*t+3/10*Pi)

> plot({P,seq(W+n*eval(W,t=0),n=-3..0)},t=0..2*Pi);

[Maple Plot]

Integrationskonstante

> eval(W,t=0);

-1/4*sin(3/10*Pi)

Und das Schaubild legt die Vermutung nahe, dass für alle Stammfunktionen gilt:

> 'P'=Diff('W',t);

> P=diff(W,t);

P = Diff(W,t)

sin(t)*sin(t+3/10*Pi) = 1/2*cos(3/10*Pi)-1/2*cos(2*...

> combine(diff(W,t)-P,trig);

0

>

Zur Übung noch einmal das bestimmte Integral als Grenzwert:

> restart:

> dx:=(b-a)/n;

dx := (b-a)/n

> x:=a+i*dx;

x := a+i*(b-a)/n

> bestint:=limit(sum(f(x)*dx,i=1..n),n=infinity);

bestint := limit(sum(f(a+i*(b-a)/n)*(b-a)/n,i = 1 ....

> f:=t->t^2; bestint;

f := proc (t) options operator, arrow; t^2 end proc...

1/3*b^3-1/3*a^3

> f:=t->t^m; bestint;

f := proc (t) options operator, arrow; t^m end proc...

-(a^(m+1)-b^(m+1))/(m+1)

> f:=t->sin(omega*t); simplify(bestint);

f := proc (t) options operator, arrow; sin(omega*t)...

-(cos(omega*b)-cos(omega*a))/omega

> f:=t->exp(-t^2); bestint;

f := proc (t) options operator, arrow; exp(-t^2) en...

-1/2*sqrt(Pi)*erf(a)+1/2*sqrt(Pi)*erf(b)

> int(f(t),t);

1/2*sqrt(Pi)*erf(t)

>

Einfache Anwendungen

Krummlinig begrenzte Flächen

Ein Aufgabentyp, mit dem man viele Schülergenerationen geplagt hat, ist die Berechnung von Flächen zwischen Kurven. Nun trifft Maple die Fallunterscheidungen.

> restart:with(plottools):

> flaeche:=int(abs(f(x)-g(x)),x=a..b);

flaeche := int(abs(f(x)-g(x)),x = a .. b)

Zum Vergleich die orientierte Fläche:

> orflaeche:=int(f(x)-g(x),x=a..b);

orflaeche := int(f(x)-g(x),x = a .. b)

> f:=x->x: g:=x->x/2+1:

> flaeche;

PIECEWISE([b-1/4*b^2, b <= 2],[1/4*b^2-b+2, 2 < b])...

> orflaeche;

1/4*b^2-1/4*a^2-b+a

Ein etwas anspruchsvolleres Beispiel:

> f:=x->-x^3+4*x;

f := proc (x) options operator, arrow; -x^3+4*x end...

> g:=x->x^4-4*x^2+1;

g := proc (x) options operator, arrow; x^4-4*x^2+1 ...

>

Berechnung der Schnittpunkte

> sols:=evalf(solve(f(x)=g(x)));

sols := .2090569265, 1.956295201, -1.338261213, -1....

> eval(flaeche,{a=min(evalf(sols)),b=max(evalf(sols))});

4.630742281-2*RootOf(_Z^3-4*_Z+_Z^4-4*_Z^2+1,.20905...
4.630742281-2*RootOf(_Z^3-4*_Z+_Z^4-4*_Z^2+1,.20905...
4.630742281-2*RootOf(_Z^3-4*_Z+_Z^4-4*_Z^2+1,.20905...
4.630742281-2*RootOf(_Z^3-4*_Z+_Z^4-4*_Z^2+1,.20905...
4.630742281-2*RootOf(_Z^3-4*_Z+_Z^4-4*_Z^2+1,.20905...
4.630742281-2*RootOf(_Z^3-4*_Z+_Z^4-4*_Z^2+1,.20905...

> evalf(%);

8.412987595

Zur grafischen Darstellung benötigen wir noch ein paar Angaben

> a:=min(evalf(sols)); b:=max(evalf(sols));

a := -1.827090915

b := 1.956295201

> n:=round((b-a)/dx);

n := round(3.783386116*1/dx)

> dx:=1/100;n;

dx := 1/100

378

> zahl:=convert(flaeche,name);

zahl := `8.412987592`

> A := seq(rectangle([x,f(x)], [x+dx,g(x+dx)], color=yellow),x=seq(a+i*dx,i=0..n)):
plots[display](A,plot({f(x),g(x)},x=a-0.2..b+0.2), style=patchnogrid,title=`A = `||zahl);

[Maple Plot]

>

>

Rotationskörper

Auch Rotationskörper sind ein beliebtes Thema.

Folgende Fläche soll um die x-Achse rotieren

> restart: with(plots): with(plottools):

Warning, the name changecoords has been redefined

> kurven:=plot({sqrt(x),sqrt(2*(x-1))},x=0..3,thickness=2,color=black):

> schnitt2d:=polygonplot([[0,0],seq([x,sqrt(2*(x-1))],x=seq(1+i/20,i=0..20)),seq([x,sqrt(x)],x=seq(2-i/20,i=0..40))],axes=normal,tickmarks=[2,2],color=cyan):

> display(kurven,schnitt2d);

[Maple Plot]

>

Dreidimensional sieht das ungefähr so aus

> p1:=plot3d([t,sqrt(t)*sin(phi),sqrt(t)*cos(phi)],t=0..2,phi=0..3*Pi/2,axes=normal,scaling=constrained,color=yellow):

> p2:=plot3d([t,sqrt(2*(t-1))*sin(phi),sqrt(2*(t-1))*cos(phi)],t=1..2,phi=0..3*Pi/2,axes=normal,scaling=constrained,color=grey):

> #display(p1,p2);

> schnitt:=polygonplot3d([[0,0,0],seq([x,0,sqrt(2*(x-1))],x=seq(1+i/20,i=0..20)),seq([x,0,sqrt(x)],x=seq(2-i/20,i=0..40))],axes=normal,color=red,tickmarks=[2,2,2],color=cyan):

>

> display(p1,p2,schnitt,orientation=[-130,60],scaling=unconstrained);

[Maple Plot]

>

Und wie groß ist nun das Volumen?

>

> RotVol:=(radius,a,b)->Pi*int(radius(x)^2,x=a..b);

RotVol := proc (radius, a, b) options operator, arr...

> Aussen:=RotVol(x->sqrt(x),0,2);

Aussen := 2*Pi

> Innen:=RotVol(x->sqrt(2*(x-1)),1,2);

Innen := Pi

> Aussen-Innen;

Pi

>

Uneigentliche Integrale

Bei der Integration über ein unbeschränktes Intervall ist manchmal eine Fallunterscheidung gefragt:

> restart:

> int(exp(a*x),x=0..infinity);

Definite integration: Can't determine if the integral is convergent.

Need to know the sign of --> -a

Will now try indefinite integration and then take limits.

limit((exp(a*x)-1)/a,x = infinity)

> int(exp(a*x),x=0..xi);

(exp(a*xi)-1)/a

>

Integriert man über Singularitäten, so kann es sein, dass das Integral nicht definiert ist:

> int(1/sqrt(abs(x)),x=-1..1);

4

>

Das ging gut. Aber hier gibt es Probleme:

> int(1/x^3, x=-3..2);

int(1/(x^3),x = -3 .. 2)

> int(1/x^3, x=-3..2, CauchyPrincipalValue);

-5/72

>

Das heißt ausführlich:

> limit(int(1/x^3,x=-3..-x0)+int(1/x^3,x=x0..2),x0=0);

-5/72

Oder

> int(1/x^3,x=-3..-x0)+int(1/x^3,x=x0..2);

1/18*(-9+x0^2)/(x0^2)-1/8*(x0^2-4)/(x0^2)

> simplify(%);

-5/72

Wobei die Grenzwerte einzeln nicht existieren

> limit(int(1/x^3,x=-3..x0),x0=0,left)+limit(int(1/x^3,x=x0..2),x0=0,right);

undefined

>

Weitere Beispiele

> int(1/(x^2+x), x=-2..1, CauchyPrincipalValue);

-2*ln(2)

> evalf(%);

-1.386294361

> plot( 1/(x^2+x), x=-2..1, -20..20);

[Maple Plot]

>

> int(1/x, x=-2..1);

int(1/x,x = -2 .. 1)

> int(1/x,x=-2..1,CauchyPrincipalValue);

-ln(2)

> evalf(%);

-.6931471806

Aber

> int(1/x^2, x=-2..1);

infinity

>

Unstetigkeitsstellen mit endlichem Sprung

Über nicht singuläre Unstetigkeitsstellen kann ohne weitere Angaben integriert werden.

> restart:

> f:=cos(x)/sqrt(1-sin(x));

f := cos(x)/(sqrt(1-sin(x)))

> discont(f,x);

{1/2*Pi+2*Pi*_Z1}

> limit(f,x=Pi/2,left), limit(f,x=Pi/2,right);

sqrt(2), -sqrt(2)

> F:=int(f,x);

F := -2*sqrt(1-sin(x))

> int(f, x=0..Pi), int(f, x=0..Pi/2), int(f, x=Pi/2..Pi);

0, 2, -2

> plot({f,int(f, x)}, x=0..Pi,discont=true);

[Maple Plot]

>

>

Integraltabellen und Integrationsregeln

Integraltabellen

Nach dem Motto 'Jedem seine Tabelle' kann man sich eine Sammlung seiner Lieblingsintegrale zusammenstellen. Auch wenn man jedes Integral neu berechnen kann, ist das oft zu Vergleichszwecken sinnvoll, besonders wenn es sich nicht um elementare Funktionen handelt, für die man oft zusätzliche Annahmen machen muss. Ein Spreadsheet leistet hier gute Dienste, weil die Berechnung automatisiert werden kann. Einfach in der ersten Spalte eine neue oder leicht geänderte Funktion eintragen und das Spreadsheet neu auswerten. Zusätzlich können die Zellen des Spreadsheets von aussen angesprochen werden.

> restart:

f(x) Int(f(x),x) Int(Int(f(x),x),x)
x^n x^(n+1)/(n+1) x^(n+2)/((n+1)*(n+2))
exp(-x^2) 1/2*sqrt(Pi)*erf(x) 1/2*sqrt(Pi)*(x*erf(x)+exp(-x^2)/(sqrt(Pi)))
1/(sqrt(1-x^2)*sqrt(1-k^2*x^2)) EllipticF(x,csgn(k)*k) ein*etwas*umfangreicher*Term
Dirac(x)*f(x) f(0)*Heaviside(x) f(0)*Heaviside(x)*x

>

>

>

Maple macht das so ähnlich: Es führt eine Tabelle über die Integrale, die während einer Session berechnet wurden:

op(4,eval(`int/int`));

TABLE([([f(0)*Heaviside(x), x], 10, _EnvCauchyPrinc...
TABLE([([f(0)*Heaviside(x), x], 10, _EnvCauchyPrinc...
TABLE([([f(0)*Heaviside(x), x], 10, _EnvCauchyPrinc...
TABLE([([f(0)*Heaviside(x), x], 10, _EnvCauchyPrinc...
TABLE([([f(0)*Heaviside(x), x], 10, _EnvCauchyPrinc...
TABLE([([f(0)*Heaviside(x), x], 10, _EnvCauchyPrinc...
TABLE([([f(0)*Heaviside(x), x], 10, _EnvCauchyPrinc...

>

>

Integrationsregeln

Die Integrationsregeln sind seit Maple eigentlich kein Thema mehr, deshalb stehen sie auch im Package für Studenten. Dennoch soll es ja vorkommen, dass auch Maple eine Stammfunktion, die es eigentlich geben müsste nicht findet. Dann kann man sich immerhin bei einer Substitution oder partiellen Integration beraten lassen:

> restart: with(student):

> changevar(x=sin(u), Int(sqrt(1-x^2), x=a...b), u);

Int(sqrt(1-sin(u)^2)*cos(u),u = arcsin(a) .. arcsin...

> changevar(x=g(u), Int(f(x), x=a...b), u);

Int(f(g(u))*diff(g(u),u),u = RootOf(g(_Z)-a) .. Roo...

>

> intparts(Int(u(x)*v(x),x),u(x));

u(x)*int(v(x),x)-Int(diff(u(x),x)*int(v(x),x),x)

> int(diff(f(x),x)/f(x),x);

ln(f(x))

>

Kurvenintegrale

Kurvenintegrale können z.B. zur Berechnung krummlinig begrenzter Flächen oder zur Berechnung von Bogenlängen verwendet werden. Das Standardbeispiel ist die Ellipse.

>

Ellipsenfläche

> restart:

> kurvint:=-int(y(t)*diff(x(t),t),t=ta..tb);

kurvint := -int(y(t)*diff(x(t),t),t = ta .. tb)

> x:=t->a*cos(t)+3; y:=t->b*sin(t)+2;

x := proc (t) options operator, arrow; a*cos(t)+3 e...

y := proc (t) options operator, arrow; b*sin(t)+2 e...

> kurvint;

-1/2*a*b*cos(tb)*sin(tb)+1/2*a*b*tb-2*a*cos(tb)+1/2...

> eval(kurvint,{ta=0,tb=2*Pi});

a*b*Pi

>

Die Grafik darf nicht fehlen:

> with(plots): with(plottools):

Warning, the name changecoords has been redefined

> #implicitplot((x-3)^2/4+(y-2)^2=1,x=0..6,y=0..4,scaling=constrained,view=[0..6,0..4]);

> #solve((x-3)^2/4+(y-2)^2=1,y);

> a,b:=2,1;

a, b := 2, 1

> unten:=polygonplot([[1,0],seq([x(t),y(t)],t=seq(i*Pi/30,i=30..60)),[5,0]],color=grey):#unten;

> elli:=display(ellipse([3,2], a, b, filled=true, color=yellow),scaling=constrained):
# Fehler beim Einfärben konkaver Polygone wird durch display überdeckt (bei richtiger Reihenfolge)

>

> display(elli,unten,plot([[1,2],[5,2]],thickness=2),view=[0..6,0..4],scaling=constrained);

[Maple Plot]

> a,b:='a','b':

> repell:=eval(kurvint,{ta=0,tb=Pi});

repell := 1/2*a*b*Pi+4*a

> remell:=-eval(kurvint,{ta=Pi,tb=2*Pi});

remell := -1/2*a*b*Pi+4*a

> repell-remell;

a*b*Pi

>

Ellipsenumfang (oder Bogen für allgemeine Integrationsgrenzen)

> restart:

> bogen:=int(sqrt(diff(x(t),t)^2+diff(y(t),t)^2),t=0..2*Pi);

bogen := int(sqrt(diff(x(t),t)^2+diff(y(t),t)^2),t ...

> x:=t->a*cos(t)+3; y:=t->b*sin(t)+2;

x := proc (t) options operator, arrow; a*cos(t)+3 e...

y := proc (t) options operator, arrow; b*sin(t)+2 e...

>

> bogen;

4*EllipticE(sqrt(-(-a^2+b^2)/(a^2)))*sqrt(b^2/(a^2)...

>

Leider scheint das nicht so einfach zu sein wie die Berechnung der Fläche: Mit diesem elliptischen Integral haben sich aber schon andere Köpfe geplagt...

> assume(b>0,a>0):

> simplify(bogen);

4*EllipticE(sqrt(a^2-b^2)/a)*a

> eval(bogen,{a=2,b=1});

8*EllipticE(1/2*sqrt(3))

> evalf(%);

9.688448216

>

Aber seit es Maple gibt, kann es einem passieren, dass man wie der Reiter über den gefrorenen Bodensee kommt, jedenfalls mit festen Integrationsgrenzen. Mit einem allgemeinen Ansatz (unbestimmte Integrationsgrenzen) kann man dagegen ziemlich einbrechen:

>

> restart:

> bogen:=int(sqrt(diff(x(t),t)^2+diff(y(t),t)^2),t=ta..tb);

bogen := int(sqrt(diff(x(t),t)^2+diff(y(t),t)^2),t ...

> x:=t->a*cos(t)+3; y:=t->b*sin(t)+2;

x := proc (t) options operator, arrow; a*cos(t)+3 e...

y := proc (t) options operator, arrow; b*sin(t)+2 e...

> simplify(bogen);

(-EllipticE(cos(tb),sqrt(-(-a^2+b^2)/(a^2)))*sqrt((...
(-EllipticE(cos(tb),sqrt(-(-a^2+b^2)/(a^2)))*sqrt((...
(-EllipticE(cos(tb),sqrt(-(-a^2+b^2)/(a^2)))*sqrt((...
(-EllipticE(cos(tb),sqrt(-(-a^2+b^2)/(a^2)))*sqrt((...

> a,b:=2,1:

> simpb:=simplify(bogen);

simpb := 2*(-sqrt(1-cos(tb)^2)*EllipticE(cos(tb),1/...

'Eigentlich sollte Maple doch sehen, dass man unter geeigneten Voraussetzungen mit sin(tb)*sin(ta) kürzen kann'. Aber im ganzen Arsenal der Befehle und Annahmen gibt es keine einfache Möglichkeit, das Maple mitzuteilen. Wer den Bleistift gerne durch die Maus ersetzt, löscht in obiger Ausgabe einfach die entsprechenden Terme und kopiert das Ganze wieder in eine Eingabezeile.

Wer mit Befehlen weiterjonglieren will, kann es so versuchen:

>

> simpb:=subs(sqrt(-cos(tb)^2+1)=sin(tb),sqrt(-cos(ta)^2+1)=sin(ta),simpb);

simpb := 2*(-sin(tb)*EllipticE(cos(tb),1/2*sqrt(3))...

> simpb:=simplify(simpb);

simpb := -2*EllipticE(cos(tb),1/2*sqrt(3))+2*Ellipt...

> eval(simpb,{ta=0,tb=2*Pi});

0

>

Und nun?

Das liegt daran:

> plot(EllipticE(cos(tb),1/2*sqrt(3)),tb=0..2*Pi,numpoints=500);

[Maple Plot]

> #plot(EllipticE((tb),1/2*sqrt(3)),tb=0..1,numpoints=500);

Also

> 2*eval(simpb,{ta=0,tb=Pi});

8*EllipticE(1/2*sqrt(3))

> 4*eval(simpb,{ta=0,tb=Pi/2});

8*EllipticE(1/2*sqrt(3))

>

So geht es sicher nicht:

> eval(bogen,{ta=0,tb=2*Pi});

Error, division by zero

Aber mit der Fußgängermethode bekommt man eine Idee:

> eval(bogen,{ta=0.0001,tb=2*Pi-0.0001});

4.844024188

Das ist für eine Ellipse mit den Halbachsen 2 und 1 etwas zu kurz

> eval(bogen,{ta=0.0001,tb=Pi-0.0001});

4.844024188

> eval(bogen,{ta=0.0001,tb=Pi/2-0.0001});

2.421812094

>

Oder mit einer Näherungsformel zur Kontrolle

> evalf(Pi*(1.5*(a+b)-sqrt(a*b)));

9.694284005

Aber wozu haben wir Grenzwerte?

> tab:=limit(bogen,ta=0,right);

tab := 2*(-sqrt(1-cos(tb)^2)*EllipticE(cos(tb),1/2*...

> limit(tab,tb=Pi,left);

4*EllipticE(1/2*sqrt(3))

> 2*evalf(%);

9.688448216

>

So gesehen ist es also fast erstaunlich, dass Maple ohne allzu langes Nachdenken von einer festen Grenze zur anderen festen Grenze über den Bodensee kommt!

>

Integration komplexer Funktionen, Residuen

Analytische Funktionen können wie Funktionen einer reellen Veränderlichen integriert werden, da in diesem Fall der Wert des Integrals vom Weg unabhängig ist.

> restart:

> int(z^2,z=1..I);

-1/3-1/3*I

> int( z*cos(z^2), z=1..I);

-sin(1)

>

Im allgemeinen Fall kann der Integrationsweg z.B. parametrisiert angegeben werden:

> kint:=int(f(z(t))*diff(z(t),t),t=ta..tb);

kint := int(f(z(t))*diff(z(t),t),t = ta .. tb)

Erster Weg: Einheitskreis um den Ursprung

> z:=t->cos(t)+I*sin(t);

z := proc (t) options operator, arrow; cos(t)+I*sin...

> kint;

int(f(cos(t)+I*sin(t))*(-sin(t)+I*cos(t)),t = ta .....

> f:=z->z^2;

f := proc (z) options operator, arrow; z^2 end proc...

> kint;

cos(tb)^3+1/3*I*cos(tb)^2*sin(tb)+2/3*I*sin(tb)-I*s...
cos(tb)^3+1/3*I*cos(tb)^2*sin(tb)+2/3*I*sin(tb)-I*s...

> eval(kint,{ta=0,tb=Pi/2});

-1/3-1/3*I

Zweiter Weg: geradlinige Verbindung:

> z:=t->t+I*(-t+1);

z := proc (t) options operator, arrow; t+I*(-t+1) e...

> kint;

(-2/3-2/3*I)*(tb^3-ta^3)+2*tb^2-2*ta^2+(-1+I)*(tb-t...

> eval(kint,{ta=1,tb=0});

-1/3-1/3*I

Diese Funktion ist nicht analytisch:

> f:=z->abs(z);

f := abs

Geradliniger Weg:

> eval(kint,{ta=1,tb=0});

-1/4*sqrt(2)*ln(1+sqrt(2))-1/2+1/2*I+1/4*I*sqrt(2)*...

Auf dem Einheitskreis

> z:=t->cos(t)+I*sin(t);

z := proc (t) options operator, arrow; cos(t)+I*sin...

> eval(kint,{ta=0,tb=Pi/2});

-1+I

>

Nicht sehr hilfreich ist hier der Vorschlag von Maple, wenn man ohne die Angabe eines Wegs integriert:

> z:='z':

> int(abs(z),z=1..I);

limit(PIECEWISE([-1/2*z^2, z <= 0],[1/2*z^2, 0 < z]...

>

Man kann übrigens eine Funktion so auf Analytizität testen:

> with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> laplacian(evalc(abs(x+I*y)), [x,y]);

2*1/(sqrt(x^2+y^2))-x^2/((x^2+y^2)^(3/2))-y^2/((x^2...

> laplacian((x+I*y)^n, [x,y]);

0

>

Geschlossener Weg

Wenn das Integral vom Weg unabhängig ist, sollte das Umlaufintegral verschwinden:

> f:=z->z^n;

f := proc (z) options operator, arrow; z^n end proc...

> z:=t->cos(t)+I*sin(t);

> eval(kint,{ta=0,tb=2*Pi});

z := proc (t) options operator, arrow; cos(t)+I*sin...

0

Auch für nicht analytische Funktionen kann das Umlaufintegral verschwinden:

> f:=z->abs(z);

f := abs

> kint;

cos(tb)+I*sin(tb)-cos(ta)-I*sin(ta)

> eval(kint,{ta=0,tb=2*Pi});

0

>

Allerdings sollte das Umlaufintegral nicht verschwinden, wenn im Inneren eine Singularität liegt (und Maple sollte das bei z^n monieren):

> f:=z->1/z;

f := proc (z) options operator, arrow; 1/z end proc...

> eval(kint,{ta=0,tb=2*Pi});

0

So klappt es:

> int(f(z(t))*diff(z(t),t),t=0..2*Pi);

2*I*Pi

>

Warum klappt es mit eval nicht? Maple wendet bei der oben gewählten Formulierung von kint zuerst die Regel für die logarithmische Integration an

> kint;

ln(cos(tb)+I*sin(tb))-ln(cos(ta)+I*sin(ta))

> int(f(z(t))*diff(z(t),t),t);

ln(cos(t)+I*sin(t))

und dann verschwindet das Integral für die Grenzen 0 .. 2*Pi .

Mit der trägen Integration wird zuerst der Integrand ausgewertet und dann der richtige Wert berechnet

> eval(Int(f(z(t))*diff(z(t),t),t=ta..tb),{ta=0,tb=2*Pi});

Int((-sin(t)+I*cos(t))/(cos(t)+I*sin(t)),t = 0 .. 2...

> value(%);

2*I*Pi

>

Singularitäten

Wenn man die logarithmische Integration vor Maple versteckt, wird das Umlaufintegral um die Singularität von 1/z korrekt berechnet:

> restart:

> z:=cos(t) + I*sin(t); dz:=diff(z,t);

z := cos(t)+I*sin(t)

dz := -sin(t)+I*cos(t)

> f:=1/z;

f := 1/(cos(t)+I*sin(t))

> int( f*dz, t=0..2*Pi);

2*I*Pi

Maple kommt auch mit komplizierteren Funktionen zurecht

> f:=sin(z)/z^4;

f := sin(cos(t)+I*sin(t))/((cos(t)+I*sin(t))^4)

> int(f *dz, t=0..2*Pi);

-1/3*I*Pi

>

Residuen

Natürlich wird man bei Umlaufintegralen um Singularitäten den Residuensatz verwenden. Dennoch lohnt es sich, an dieser Stelle die Integration im Komplexen etwas näher zu beleuchten.

> restart: with(plots):

Warning, the name changecoords has been redefined

> f:=(2*z^3+3*z+50)/(z^3-3*z-52);

f := (2*z^3+3*z+50)/(z^3-3*z-52)

>

Mit discont bekommt man nur die reellen Singularitäten

> discont(f,z);

{4}

> solve(denom(f));

4, -2+3*I, -2-3*I

>

Die Residuen kann man an der Reihenentwicklung ablesen

> series(f,z=4,3);

series(38/9*(z-4)^(-1)+29/27+O((z-4)),z=-(-4),1)

>

> series(f,z=-2+3*I,3);

series((-19/9+49/18*I)*(z+2-3*I)^(-1)+(59/60+19/270...

>

Oder mit einem fertigen Befehl

> residue(f,z=4);

38/9

> residue(f,z=-2+3*I);

-19/9+49/18*I

> residue(f,z=-2-3*I);

-19/9-49/18*I

>

Nähere Untersuchung

> z:=x+I*y;

z := x+I*y

> discont(f,y);

{I*(x-4), I*(x+2+3*I), I*(x+2-3*I)}

> solve(denom(f),y);

I*x-3+2*I, -4*I+I*x, 3+2*I+I*x

> solve(denom(f),x);

-I*y+4, -2-3*I-I*y, -2+3*I-I*y

(Wenn man sich dazudenkt, dass x und y reell sein müssen, sind das die gleichen Ergebnisse wie oben).

>

Darstellung

> plr:=plot3d(Re(f),x=-8..8,y=-8..8,view=[-8..8,-8..8,-2..5],axes=framed,grid=[40,40],style=patchcontour,contours=15): plr;

[Maple Plot]

> pli:=plot3d(Im(f),x=-8..8,y=-8..8,view=[-8..8,-8..8,-2..5],axes=framed,grid=[40,40]): pli;

[Maple Plot]

>

Parameterdarstellung für die Integration

> r,x0,y0:='r','x0','y0':

> fp:=subs(x=r*cos(t)+x0, y=r*sin(t)+y0,f);

fp := (2*(r*cos(t)+x0+I*(r*sin(t)+y0))^3+3*r*cos(t)...

>

> dz:=diff(r*(cos(t)+x0+I*(sin(t)+y0)),t);

dz := r*(-sin(t)+I*cos(t))

>

> r:=5: x0:=1: y0:=0:

> spr:=spacecurve([r*cos(t)+x0, r*sin(t)+y0,Re(fp),t=0..2*Pi],color=black, thickness=2,numpoints=500, axes=normal): spr;

[Maple Plot]

>

> display(display(plr,color=yellow),spr,orientation=[100,30]);

[Maple Plot]

>

> spi:=spacecurve([r*cos(t)+x0, r*sin(t)+y0,Im(fp),t=0..2*Pi],color=black, thickness=2,numpoints=500,axes=normal): spi;

[Maple Plot]

>

> display(display(pli,shading=ZHUE),spi);

[Maple Plot]

>

Noch einmal die Residuen

> z:='z':

> r1:=residue(f,z=4); r2:=residue(f,z=-2+3*I); r3:=residue(f,z=-2-3*I);

r1 := 38/9

r2 := -19/9+49/18*I

r3 := -19/9-49/18*I

Umlaufintegrale für verschiedene Wege

a) Wie in obiger Darstellung (alle Singularitäten im Inneren)

> 2*Pi*I*(r1+r2+r3);

0

b)

> 2*Pi*I*(r2+r3);

-76/9*I*Pi

c)

> 2*Pi*I*r1;

76/9*I*Pi

>

Explizite Berechnung der entsprechenden Umlaufintegrale (das dauert etwas länger):

a)

> x0,y0,r;

> int(fp*dz,t=0..2*Pi);

1, 0, 5

0

>

Ggf. Kontrollplots

plot(Re(fp*dz),t=0..2*Pi);

plot(Re(fp),t=0..2*Pi);

plot(Im(fp*dz),t=0..2*Pi);

plot(Im(fp),t=0..2*Pi);

>

b)

> evalf(sqrt(13));

3.605551276

> st:=time():

> x0:=0: y0:=0: r:=37/10:

> int(fp*dz,t=0..2*Pi);
time()-st;

-76/9*I*Pi

18.256

c)

> evalf(sqrt(36+9));

6.708203934

> x0:=4:r:=67/10;

> int(fp*dz,t=0..2*Pi);

r := 67/10

76/9*I*Pi

> x0:=4:r:=68/10;

> int(fp*dz,t=0..2*Pi);

r := 34/5

0

>

>

Mehrfachintegrale

Weil die Welt nicht eindimensional ist, sind Mehrfachintegrale das tägliche Brot des Ingenieurs und Physikers: Die Gebiete und Räume, in denen Mehrfachintegrale eingesetzt werden, sind so umfassend, dass sie in dieser kurzen Einführung nur angedeutet werden können.

Wir beginnen mit einem einfachen Beispiel, der Berechnung des Volumens einer Pyramide:

> restart: with(plots):

Warning, the name changecoords has been redefined

> pyr:=polygonplot3d([[0,0,2],[1,0,0],[0,1,0] ,[0,0,2],[0,0,0],[1,0,0]],axes=normal,tickmarks=[2,2,2],labels=[x,y,z],color=yellow,thickness=2,orientation=[-20,70]):#pyr;

>

Mit einem 'Volumenelement'

> elem:=polygonplot3d([[1/2,0,1],[0,1/2,1],[0,0,1]],axes=normal,tickmarks=[2,2,2],labels=[x,y,z],color=red,thickness=5):

> display(elem,pyr,title=`Bitte drehen`);

[Maple Plot]

>

Für Studenten gibt es eine Kurzfassung von int(int(int...)))

> with(student):

> Tripleint(1,z=0..2-2*x-2*y,y=0..1-x,x=0..1);

Int(Int(Int(1,z = 0 .. 2-2*x-2*y),y = 0 .. 1-x),x =...

> value(%);

1/3

Das Doppelintegral kann man sich so ansehen:

> map(value,%%);

Int(2-2*x-2*x*(1-x)-(1-x)^2,x = 0 .. 1)

>

In jedem Fall kommt es auf die richtige Reihenfolge an:

> int(1,z = 0 .. 2-2*x-2*y);

2-2*x-2*y

> int(%,y=0..1-x);

2-2*x-2*x*(1-x)-(1-x)^2

> int(%,x=0..1);

1/3

>

Während man das Volumen einer Pyramide auch in der Formelsammlung nachschlagen kann, wird es beim Trägheitsmoment schon etwas komplizierter:

> Tripleint(y^2+z^2,z=0..2-2*x-2*y,y=0..1-x,x=0..1);

Int(Int(Int(y^2+z^2,z = 0 .. 2-2*x-2*y),y = 0 .. 1-...

> map(value,%);

Int(-7/6*(1-x)^4+1/3*(10-10*x)*(1-x)^3+1/2*(1/3*(2-...
Int(-7/6*(1-x)^4+1/3*(10-10*x)*(1-x)^3+1/2*(1/3*(2-...

>

> value(%);

1/6

>

Allgemein werden in Maple also Dreifachintegrale so formuliert:

> Int(Int(Int(g(x,y,z),z=z1(x,y)..z2(x,y)),y=y1(x)..y2(x)),x=x1..x2);

Int(Int(Int(g(x,y,z),z = z1(x,y) .. z2(x,y)),y = y1...

>

Oder etwas bequemer

> Tripleint(g(x,y,z),z=z1(x,y)..z2(x,y),y=y1(x)..y2(x),x=x1..x2);

Int(Int(Int(g(x,y,z),z = z1(x,y) .. z2(x,y)),y = y1...

>

>

Krummlinige Koordinaten

Kartesische Koordinaten sind in der Praxis eher die Ausnahme. Meistens müssen für die Integration die Koordinaten dem Problem angepasst werden. Dazu benötigt man die Funktionaldeterminante der Koordinatentransformation. In Maple kann mit 'jacobian' die erforderliche Matrix aufgestellt werden:

> restart: with(linalg):

Warning, the protected names norm and trace have been redefined and unprotected

> xyz:=[x(u,v,w),y(u,v,w),z(u,v,w)];

xyz := [x(u,v,w), y(u,v,w), z(u,v,w)]

> jacobian(xyz, [u,v,w]);

matrix([[diff(x(u,v,w),u), diff(x(u,v,w),v), diff(x...

> fundet:=det(jacobian(xyz, [u,v,w]));

fundet := diff(x(u,v,w),u)*diff(y(u,v,w),v)*diff(z(...
fundet := diff(x(u,v,w),u)*diff(y(u,v,w),v)*diff(z(...
fundet := diff(x(u,v,w),u)*diff(y(u,v,w),v)*diff(z(...

>

Ein bekanntes Beispiel (Kugelkoordinaten):

> x:=(u,v,w)-> u*cos(v)*sin(w);
y:=(u,v,w)-> u*sin(v)*sin(w);
z:=(u,v,w)-> u*cos(w);

x := proc (u, v, w) options operator, arrow; u*cos(...

y := proc (u, v, w) options operator, arrow; u*sin(...

z := proc (u, v, w) options operator, arrow; u*cos(...

> fundet;

-cos(v)^2*sin(w)^3*u^2-sin(v)^2*sin(w)^3*u^2-cos(w)...

> simplify(%);

-sin(w)*u^2

>

Weitere Koordinatensysteme finden Sie mit

> ?coords

>

Numerische Integration

Es gibt auch Integranden, die sich nicht symbolisch integrieren lassen

> restart:

> int( sin(1/x^3), x=1..2);

int(sin(1/(x^3)),x = 1 .. 2)

Dann ist immer noch eine numerische Berechnung möglich

> evalf(%);

.3548334332

Mit der trägen Variante Int sucht Maple erst gar nicht nach einer geschlossenen Ausdruck

> evalf(Int( sin(1/x^3), x=1..2));

.3548334332

Man kann die gewünschte Stellenzahl und verschiedene numerische Methoden angeben

> evalf(Int(sin(x)/x, x=1..100, 15, _NCrule));

.616142396521873

> evalf(Int(sin(x)/x, x=1..100, 15));

.616142396521873

>

Wenn die symbolische Integration nicht durchgeführt wird, kann das auch an RootOfs liegen

> f:=1/(1+x^7);

f := 1/(1+x^7)

> int( f, x=1..2);

1/7*ln(3)-1/7*sum(_R*ln(2-_R),_R = RootOf(_Z^6-_Z^5...
1/7*ln(3)-1/7*sum(_R*ln(2-_R),_R = RootOf(_Z^6-_Z^5...

> value(%);

1/7*ln(3)-1/7*sum(_R*ln(2-_R),_R = RootOf(_Z^6-_Z^5...
1/7*ln(3)-1/7*sum(_R*ln(2-_R),_R = RootOf(_Z^6-_Z^5...

>

> evalf(%);

.1163017046+0.*I

>

Oder direkt

> evalf( Int(f, x=1..2), 20);

.11630170437709733349

>

Uneigentliche Integrale numerisch berechnen?

> int( 1/x^2, x=0..infinity);

infinity

> evalf( Int( 1/x^2, x=0..infinity) );

Float(infinity)

>

Man kann auch ein falsche Ergebnis provozieren:

> evalf( Limit(Int( 1/x^2, x=n..infinity), n=0) );

2.000000000

> Limit(Int( 1/x^2, x=n..infinity), n=0) ;

Limit(Int(1/(x^2),x = n .. infinity),n = 0)

> value(%);

undefined

>

> Limit(Int( 1/x^2, x=n..infinity), n=0,right) ;

Limit(Int(1/(x^2),x = n .. infinity),n = 0,right)

> value(%);

infinity

> evalf(%);

Float(infinity)

>

> evalf( Limit(Int( 1/x^2, x=n..infinity), n=0,right) );

2.000000000

>

Kontrolle der Integration

Die Ableitung der Stammfunktion sollte wieder der Integrand sein. Aber manchmal ist dieser Nachweis mit Maple nicht ganz einfach, besonders wenn trigonometrische Funktionen im Spiel sind.

> restart:

> f:=sin(a+1/x);

f := sin(a+1/x)

> inte:=int(f,x);

inte := sin(a+1/x)*x-Si(-1/x)*sin(a)-Ci(1/x)*cos(a)...

> integrand:=diff(inte,x);

integrand := -cos(a+1/x)/x+sin(a+1/x)-sin(1/x)*sin(...

>

> simplify(integrand);

(-cos((a*x+1)/x)+sin((a*x+1)/x)*x-sin(1/x)*sin(a)+c...

> expand(integrand);

sin(a)*cos(1/x)+cos(a)*sin(1/x)

> combine(%,trig);

sin(a+1/x)

>

Falls die Termumformung mit Maple beim besten Willen nicht gelingt (man kann Stunden damit zubringen, die passenden Befehle in der richtigen Reihenfolge abzusetzen), und die Terme so umfangreich sind, dass auch eine Rechnung von Hand weitere Stunden kosten würde, gibt es weitere Möglichkeiten:

>

> testexpr:=integrand-f;

testexpr := -cos(a+1/x)/x-sin(1/x)*sin(a)/x+cos(1/x...

> solve(testexpr);

{a = a, x = x}

Demnach ist testexpr für alle x und a gleich Null.

Man kann aber auch mit Zufallszahlen eine an Sicherheit grenzende Wahrscheinlichkeit berechnen:

>

> test:='5000-rand()/1e8';

test := 5000-.1e-7*rand()

> seq( evalf(subs(x=test, a=test, testexpr)), i=1..5);

-.2893e-9, .24e-11, .6581e-9, .450e-10, .2134e-9

> map(abs, [%] );

[.2893e-9, .24e-11, .6581e-9, .450e-10, .2134e-9]

> max( op(%) );

.6581e-9

>

>

Maple beim Integrieren zusehen

> restart:

> infolevel[int]:=3:

> int( 1/(1+ln(x)), x);

int/indef1: first-stage indefinite integration

int/indef2: second-stage indefinite integration

int/ln: case of integrand containing ln

int/indef1: first-stage indefinite integration

int/indef2: second-stage indefinite integration

int/exp: case of integrand containing exp

-exp(-1)*Ei(1,-ln(x)-1)

> infolevel[`evalf/int`]:=2:

> evalf(Int(sin(x)/x, x=1..2));

evalf/int/control: integrating on 1 .. 2 the integrand

sin(x)/x

evalf/int/control: Trying easyproc

evalf/int/control: from ccquad: error = .9167111514330e-12
integrand evals = 19
result = .6593299064355

.6593299064

>

Zurück zum Inhaltsverzeichnis