LEHRSTUHL A FÜR MATHEMATIK                                           [23.mws]
    Prof. Dr. E. Görlich
                                                                         Sommersemester 2000
                                                                         08.06.2000

                                             MAPLE  für ANALYSIS II
                                             =====================
    22) Darstellung von Mengen mit "geometry[draw]" oder durch Hilfsprogramm, parametrische
          Plots, Banach-Fixpunktsatz, MMMittelwertsatz, totales Differential.
          (interface(plotdevice=gl), shading = zhue, parametric plot, scaling=constrained, solve,
           matrix, vector, D, norm( ,frobenius), fsolve, jacobian, eval, evalm, map)
    ########################################################################
    Aus Uebung 9, Aufgabe 1:
    a) Menge {(x,y); abs(x)<=1, abs(y)<=1,abs(x)+abs(y)<1}
    > restart:
    > ?plot
    > ?geometry[draw]
    > ?geometry[square]
    > with(geometry):
    > point(A,-1,-1):point(B,-1,1):point(C,1,1):point(E,1,-1):
    > square(sq1,[A,B,C,E]):
    > geometry[draw]([square(sq1,[A,B,C,E])]);
    > point(F,0,-2):point(G,-2,0),point(H,0,2):point(J,2,0):
    > geometry[draw]([square(sq1,[A,B,C,E]),square(sq2,[F,G,H,J])]);
    Mengendarstellung mittels Hilfsprogramm.
    ===============================
    Eine andere Moeglichkeit zur Visualisierung von Mengen: Definiere eine Funktion (procedure)
    auf der Ebene, die den Wert 1 auf der Menge hat und den Wert 0 sonst. Dann plotte sie mit
    der Option "shading = zhue".
    Vorteil: Man braucht nicht die Eckpunkte auszurechnen und sie in das "geometry[draw]"-
    Programm einzugeben. Nachteil: Aufwendiger als obige Methode: Man braucht eine relativ
    kleine Maschenweite um ein einigermassen scharfes Bild zu erhalten.
    a)
    > MengeA:=proc(x,y) if abs(x)<=1 and abs(y)<=1 and abs(x)+abs(y)<1 then 1 else 0 fi; end;
    > plot3d(MengeA, -2..2,-2..2);
    > ?plot3d,options
    > plot3d(MengeA, -2..2,-2..2,shading=zhue,orientation=[0,0],grid=[80,80]);
    Wiederholung aus Paragraph 16: plotdevice = gl unter Linux:
    > interface(plotdevice=gl):
    > plot3d(MengeA, -2..2,-2..2,shading=zhue,grid=[80,80]);
    > interface(plotdevice=default):
    b) logarithmische Spirale (Forster 2, S. 32)
    > #?implicitplot
    > #?curveplot
    > #?parametric
    > #?plots
    > #?plot,options
    > ?plot,parametric
    > restart:
    > plot([exp(t)*cos(t),exp(t)*sin(t),t=-2*Pi..2*Pi]);
    > plot([exp(t/2/Pi)*cos(t),exp(t/2/Pi)*sin(t),t=-2*Pi..2*Pi],scaling=constrained);
    > plot([exp(t/2/Pi)*cos(t),exp(t/2/Pi)*sin(t),t=-6*Pi..6*Pi],scaling=constrained);
    > plot([exp(t/2/Pi)*cos(t),exp(t/2/Pi)*sin(t),t=-infinity..0],scaling=constrained);
    #####################################################################################
    Zu Uebung 9, Aufgabe 4:
    Fixpunktsatz von Banach. (Satz IX.3.7 der Vorlesung mit V=(IR^n;,|| . ||2))
    > equ:={arctan(x+y)=3*x,exp(-(x^2+y^2))=5*y};
    > ?solve
    > solve(equ,{x,y});
    > F:=(x,y)->vector([arctan(x+y)/3,exp(-(x^2+y^2))/5]);
    > Iterate:=proc(a,n::posint) local j,result; result:=F(a[1],a[2]); for j from 1 to n do
    result:=F(result[1],result[2]) od; end;
    > Iterate(vector([1,1]),1);
    > evalf(%);
    > Iterate2:=proc(a,n::posint) local j,result; result:=F(a[1],a[2]); for j from 1 to n do
    result:=evalf(F(result[1],result[2])) od; end;
    > Iterate2(vector([1,1]),1);
    > Iterate2(vector([1,1]),2);
    > Iterate2(vector([1,1]),10);
    > Iterate2(vector([1,1]),20);
    > Iterate2(vector([-1,1]),20);
    > Iterate2(vector([1,-1]),20);
    > Iterate2(vector([-1,-1]),20);
    Lipschitzbedingung:
    Abschaetzung der Lipschitzkonstante mittes Satz (X.3.1) der Vorlesung
    (Mittelwertsatz fuer Funktionen f: IR^n->IR^n) (Elegantere Methode siehe unten!)
    > ?D
    > with(linalg):
    Zuerst muss F neu definiert werden:
    > TD:=matrix([D[1](F[1]),D[2](F[1]),D[1](F[2]),D[2](F[2])]);
    > D[1](F[1](x,y));
    > F:=vector([(x,y)->arctan(x+y)/3,(x,y)->exp(-(x^2+y^2))/5]);
    > D[1](F[1]);
    > D[2](F[1]);
    > D[1](F[2]);D[2](F[2]);
    > D[1](F[2])(x,y);
    Totales Differential, "von Hand":
    > TD:=matrix([[D[1](F[1])(x,y),D[2](F[1])(x,y)],[D[1](F[2])(x,y),D[2](F[2])(x,y)]]);
    > ?linalg,norm
    > norm(%,frobenius);
    Eine offensichtliche Majorante dieser Funktion im Einheitsquadrat ergibt sich so:
    > #1/15*sqrt(50*1/(abs(1+(x+y)^2)^2)+36*exp(-Re(x^2+y^2))^2*abs(x)^2+36*exp(-Re(x^2+y^2))^2*
    abs(y)^2)
    Ersetze abs(x) und abs(y) durch 1 :
    > #1/15*sqrt(50*1/(abs(1+(x+y)^2)^2)+36*exp(-Re(x^2+y^2))^2+36*exp(-Re(x^2+y^2))^2)
    Setze nun x=y=0:
    > subs({x=0,y=0},1/15*sqrt(50*1/(abs(1+(x+y)^2)^2)+36*exp(-Re(x^2+y^2))^2+36*exp(-Re(x^2+y^2))^2));
    > evalf(%);
    -->> Lipschitzkonstante ist kleiner als 1, also Banach-Fixpunktsatz anwendbar.
    Gegenprobe mit "fsolve":
    > ?fsolve
    > fsolve(equ,{x,y});
    Das ist die gleiche Loesung. Dass nur eine Loesung existiert, folgt aus dem Satz von Banach.
    ######################################################################
    Eleganter: Maple hat auch einen Befehl zur Berechnung einer Jacobi-Matrix (Def. X.2.4):
    > ?jacobian
    > F:=[1/3*arctan(x+y) , 1/5*exp(-y^2-x^2)];
    > A:=jacobian(F,[x,y]);
    Wir moechten x=0,y=0 einsetzen. das geht nicht mit "eval" oder "evalm":
    > x:=0:y:=0:
    > evalm(A);
    > eval(A);
    Verwende dazu "map":
    > map(eval,A);
    > x:='x':y:='y':
    -->> Ergebnis wie zuvor