# Lektion 4
# Funktionen und ihre Darstellungen
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(plots):
# Unter einer Funktion wird eine Vorschrift verstanden, die jedem Element x einer Menge  D  genau ein Element y aus einer Menge W zuordnet. Man schreibt y = y(x) oder  auch  y = f(x).
# x:   Unabhängige Veränderliche
# y:   Abhängige Veränderliche
# D:  Definitionsbereich der Funktion
# W: Wertebereich der Funktion
# Beispiele:
# 1.) Gerade: Definitionsbereich: D = (-,); Wertebereich: W = (-,);
> y:=2*x+1;
> plot(y,x=-2..2,-3..3,thickness=2,gridlines=true,scaling=constrained);
# 2.) Parabel: Definitionsbereich: D = (-,); Wertebereich: W = [0,)
> y:=x^2;
> plot(y,x=-2..2,thickness=2,gridlines=true,scaling=constrained);
# y = f(x)        heißt explizite  Darstellung (die Funktion ist nach einer Variablen aufgelöst)
# F(x,y) = 0    heißt implizite  Darstellung (die Funktion ist nicht nach einer Variablen aufgelöst)
# Beispiele für implizite  Darstellungen von Funktionen sind die Kurven zweiter Ordnung:  F(x,y) = a x2 + 2b xy + c y2 + 2d x + 2e y + f  =  0. 
> x:='x': y:='y':
> F:=unapply(x^2*a+2*x*y*b+y^2*c+2*d*x+2*e*y+f,x,y);
# Für die folgenden Rechnungen setzen wir
> f:=-1; d:=0; e:=0;
# 1. Ellipse:   a = 1;  b = -0.5;  c = 1;  (ac - b2 = 0.75 > 0);
> a:=1.; b:=-0.5; c:=1.; a*c-b^2; F(x,y);
> implicitplot(F(x,y),x =-2..2,y=-2..2,scaling=constrained,thickness=2,axes=boxed,gridlines=true);
# Das geht noch etwas glatter:
> implicitplot(F(x,y),x=-2..2,y=-2..2,scaling=constrained,thickness=2,axes=boxed,gridlines = true,gridrefine = 2,crossingrefine = 3);
# 2. Kreis:   b = 0;  (ac -b2  =  1 > 0)
> b:=0.; a*c-b^2; F(x,y);
> implicitplot(F(x,y),x =-1..1,y=-1..1,thickness=2,axes=boxed,gridlines=true,scaling=constrained);
# 3. Hyperbel:   a = 1; b = 0.5; c = -1;  (ac - b2  = -1.25 < 0);
> a:=1; b:=0.5; c:=-1.; a*c-b^2; F(x,y);
> implicitplot(F(x,y),x =-5..5,y=-5..5,thickness=2,axes=boxed,gridlines = true,gridrefine = 2,scaling=constrained);
# 4. Parabel:  d = 3; e = -1; f = -1;  (ac - b2 = 0);
> d:=3;e:=-1;a:=-1;b:=1;c:=-1; a*c-b^2; F(x,y);
> implicitplot(F(x,y),x =-2..2,y=-2..2,thickness=2,axes=boxed,gridlines=true,scaling=constrained);
# Als Beispiel für die Parameterdarstellung einer ebenen Kurve betrachten wir die Zykloide. Ihre Parameterdarstellung lautet mit t als Parameter: x = r (t - sin t);  y = r (1 - cos t), und für r = 1 ergibt sich folgendes Bild
> plot([t-sin(t),1-cos(t),t=0..4*Pi],thickness=2,scaling=constrained,gridlines=true);
# Für die Darstellung von Raumkurven existiert in Maple der Befehl  spacecurve  aus der Bibliothek plots.  Wir wählen die Parameterdarstellung einer Schraubenlinie (Helix)
> spacecurve([cos(t),sin(t),0.1*t,t=0..8*Pi],axes=normal);
# und für eine größere Anzahl von Stützpunkten (der Standardwert ist n = 50): 
> spacecurve([cos(t),sin(t),0.1*t,t=0..8*Pi],axes=normal,numpoints=100,thickness=2);
# Eine weitere Darstellungsmöglichkeit erlaubt der Befehl  tubeplot , ebenfalls aus der Bibliothek plots.   
> tubeplot([cos(t),sin(t),0.1*t,t=0..8*Pi],axes=NORMAL,style=patchnogrid,numpoints=100,radius=0.03);
# Wir beschäftigen uns im Folgenden mit der Lösung von Gleichungen und beginnen mit der Nullstellenbestimmung von Polynomen 2. Grades. Dazu benutzen wir den Befehl solve 
> pq:=x^2+p*x+q;
> loepq:=solve(pq,x);
# Hier liefert Maple die analytische Lösung (pq-Formel). Als weiteres Beispiel betrachten wir das Polynom 
> p1:=x^2+9;
> plot(p1,thickness=2,axes=boxed,gridlines=true);
# und ermitteln die beiden rein imaginären Wurzeln:
> solve(p1,x);
# Der Befehl lässt sich auch abgekürzt schreiben:
> solve(p1);
# Das folgende Polynom hat zwei reelle Nullstellen.
> p2:=x^2-6*x+6;
> plot(p2,x=0..6,y=-5..5,thickness=2,axes=boxed,gridlines=true);
> solve(p2);
# Und hier ein Polynom vom Grade 5:
> p3:=x^5-9*x^4+24*x^3-14*x^2-16*x+32;
> plot(p3,x=-1.5..5,y=-10..50,thickness=2,axes=boxed,gridlines=true,title="\nGraf des Polynoms p3\n",titlefont=["ROMAN",15]);
#  Wir faktorisieren p3 
> fp3:=factor(p3);
# Wir entnehmen dieser Darstellung sofort die einfache Nullstelle x = -1 sowie die doppelte Nullstelle  x = 4. 
# Sollen die Nullstellen der quadratischen Gleichung f(x) = x^2-2*x+2*(zweiter*Faktor in fp3)*ber;echnet werden, dann geht das wie folgt
> op2:=op(2,fp3);
> solve(op2);
# Die vollständige Lösungsmenge des Polynoms p3 erhalten wir mit dem Befehl
> solve(p3);
# Zur weiteren Bearbeitung und Auswertung von Polynomem siehe Polynomials in der Hilfe-Funktion. Wir ändern jetzt das letzte Beispiel nur geringfügig ab
> p4:=p3-1;
# und faktorisieren diesen Ausdruck (was hier nichts Neues ergibt) 
> factor(p4);
> plot(p4,x=-1.5..5,y=-10..50,title="\nGraf des Polynoms p4\n",titlefont=["ROMAN",15],thickness=2,axes=boxed,gridlines=true);
# Wir bestimmen die Nullstellen mit  
> N4:=solve([p4],[x]);
# Der Ausdruck  RootOf   in der obigen Lösung ist ein Platzhalter für alle Lösungen der Gleichung p4. Maple signalisiert mit dieser Ausgabe, dass allgemeine Wurzeln dieses Polynoms 5. Grades nicht existieren. Die nummerischen Näherungen aller Wurzeln erhalten wir mit dem Befehl evalf
> evalf(N4);
> evalf(N4[1]);evalf(N4[1][1]);evalf(N4[1,1]);
# Sind wir nur an den reellen Lösungen interessiert, dann hilft der Befehl Polynomial  aus der Bibliothek SolveTools  mit dem Zusatz domain = real: 
> SolveTools[Polynomial](p4, x, domain = real);
> evalf(%);
# Die nummerischen Werte der reellen Lösungen unserer  Gleichungen können wir auch unmittelbar mit dem Befehl  fsolve  erhalten. Als Lösung von p1 = 0, die rein imaginäre Lösungen besitzt, wird dann natürlich die NULL-Sequenz ausgegeben. 
> fsolve(p1);
# Will man auch die nummerischen Werte der komplexen Wurzeln haben, so errreicht man das mit einem Zusatz im Aufruf.
> fsolve(p1,complex);
> fsolve(p3); fsolve(p3,complex);
> fsolve(p4,complex);
# Die Quadratwurzel als Umkehrfunktion der zweiten Potenz
# Die Normalparabel y = x^(2) (x>=0)  ;ist im 1. Quadranten streng monoton wachsend und besitzt dort eine Umkehrung  y=sqrt(x). ;Damit wird nur der Ast der Normalparabel im 1. Quadranten gespiegelt. Für x > 0 bedeutet y =  sqrt(x);  diejenige positive Zahl, deren Quadrat gleich x ist. Für die negative Zahl y, deren Quadrat auch gleich  x  ist, schreiben wir y = -sqrt(x); . Wir schauen uns folgende Grafen an
> plot([x^2,sqrt(x),x],x=-1.5..1.5,linestyle=[solid,solid,dashdot],thickness=[2,2,1],color=[black,blue,red],scaling=constrained,axes=boxed,gridlines=true);
# Sollen die Wurzeln aus einer Quadratzahl berechnet werden, also  y = sqrt(x^2); , dann lösen wir die Gleichung y^2 = x^2, also;
> solve(y^2=x^2,y);
# Wenn also das Vorzeichen der reellen Variablen x  nicht bekannt ist, dann kann es falsch sein,  sqrt(x^2);  durch  x  zu ersetzen. Maple trifft daher diese Vereinfachung nicht. Der Ausdruck 
> g:=sqrt(x^2);
# wird bei Anwendung des Befehls simplify   korrekt vereinfacht
> NULL;simplify(g);
# Teilen wir  Maple mit, dass es sich bei x um eine reelle Zahl handelt, dann folgt 
> NULL;simplify(g, assume = real);
# und nach weiterer Spezifizierung, wenn wir  x  als positiv voraussetzen
> NULL;simplify(g, assume = positive);
# Der Befehl symbolic  teilt Maple mit, dass eine formale symbolische Manipulation erlaubt ist. Dabei wird keine Rücksicht auf die möglichen Werte von  x  genommen. 
> NULL;simplify(g, symbolic);
# Die Dritte Wurzel als Umkehrung der dritten Potenz.
# Da auch die Funktion y = x^3; über den gesamten Wertebereich der reellen Achse x  nicht streng monoton verläuft, das gilt nur für x >= 0;, wird auch nur der Ast im ersten Quadranten gespiegelt. 
> plot([x^3,x^(1/3),x],x=-1.2..1.2,color=[black,blue,red],linestyle=[solid,solid,dashdot],thickness=[2,2,1],scaling=constrained,axes=boxed,gridlines=true);
# Für negative Werte von x  findet Maple keine reelle dritte Wurzel. Wir lassen uns das an einem Beispiel bestätigen:
> (-8)^(1/3);
# Wir nähern uns dem Problem von einer anderen Seite und lösen die  Polynomgleichung  y^3 = -8; :
> y:='y': solve(y^3=-8);
# Maple liefert uns alle drei Lösungen, von denen eine reell und zwei konjugiert komplex sind. 
# Will man als dritte Wurzel einer negativen Zahl  x  stets eine negative reelle Zahl erhalten, so kann man den Befehl  surd  einsetzen. Für  x < 0;  liefert  surd(x, n) = -(-x)^(1/n);.
> x:=-8; x^(1/3); simplify(%); surd(x,3);
# Lösungen Transzendenter Gleichungen. 
# Wir betrachten dazu zwei Beispiele
> x:='x': 
> solve(sin(x)=1/2); solve(cosh(x)=3);
# Nummerische Lösungen erhält man wieder mit dem Befehl  fsolve.
> fsolve(sin(x)=1/2,x); fsolve(cosh(x)=3,x);
> plot([sin(x),1/2],x=0..3*Pi,thickness=[2,1],color=[black,blue],axes=boxed,gridlines=true);
# Wir können vorgeben, in welchem abgeschlossenen Intervall Maple die nummerische Lösung suchen soll. 
> fsolve(sin(x)=1/2,x=Pi/2..Pi);
# Als Beispiel für das Auftreten einer transzendenten Gleichung in der Technischen Mechanik betrachten wir die Berechnung der kleinsten Knicklast für den 3. Euler-Fall : 
> z:='z': transz:=sin(z)-z*cos(z); solve(transz);
# Wir benötigen eine nichttriviale Lösung z <> 0;. Aufschluss gibt uns eine graphische Darstellung der Eigenwertgleichung
> plot(transz,z=-6..20,thickness=2,axes=boxed,gridlines=true);
# Die kleinste positive Lösung liegt zwischen 4 und 5. Ihren nummerischen Wert erhalten wir mit dem Aufruf
> z1:=fsolve(transz,z=4..5);
# Einen besseren Überblick über die Lage sämtlicher positiver Wurzeln unserer transzendenten Gleichung bekommen wir, wenn wir die Eigenwertgleichung durch cos z ` 0 dividieren. Das Ergebnis ist   tan z = z.  Wir stellen beide Seiten dieser Gleichung graphisch dar und können aus den Schnittpunkten der Kurven Näherungen für die Wurzeln entnehmen.
> plot([tan(z),z],z=0..3*Pi,-5..10,axes=boxed,gridlines=true);
# Die Lösung von Gleichungssystemen erfolgt ebenfalls mit dem Befehl solve. Wir wählen als Beispiel eine Gleichung zweiten und eine ersten Grades.
> x:='x': y:='y':      
> gl1:=2*y-x+1=0; gl2:=x^2+y^2-x*y-1=0; 
> solve([gl1,gl2]);
> fsolve([gl1,gl2]);
# Klarheit in der Frage, wieviele Lösungen zu erwarten sind und wo sie sich etwa befinden, gibt uns eine graphische Darstellung.
> implicitplot([gl1,gl2],x=-1.5..1.5,y=-1.5..1.5,axes=boxed,gridlines=true,scaling=constrained);
# Wir verschieben jetzt die Gerade nach unten:
> gl1:= 2*y - x + 7/4:
> plots[implicitplot]([gl1,gl2],x=-1.5..1.5,y=-1.5..1.5,axes=boxed,gridlines=true,scaling=constrained);
> solve([gl1,gl2]); assign(%):
# Die mit dem Befehl  assign  zugewiesenen (x,y)-Werte machen wir mit dem Befehl  allvalues  explizit. Die Umwandlung in nummerische Werte leistet wieder der Befehl  evalf.
> allvalues([x,y]); evalf(%);
# Mit dem Befehl   fsolve  erhalten wir lediglich die nummerischen Werte eines (x,y)-Paares 
> x:='x':  y:='y':  fsolve([gl1,gl2]);
# Der Zeichnung entnehmen wir, dass der zweite Schnittpunkt im abgeschlossenen Rechteck [(x, y)-1 <= x and x <= 0, -2 <= y and y <= 0];zu suchen ist. Wir erhalten seine Koordinaten, indem wir Maple diesen Hinweis geben.
> x:='x': y:='y':  fsolve({gl1,gl2},{x,y},{x=-1..0,y=-2..0});
# Funktionen von zwei Variablen
# Wir betrachten als Beispiel die Durchbiegung w(¾,·) einer dünnen elastischen Platte (ebenes Flächentragwerk) unter konstanter Flächenlast p0 mit den Abmessungen 2a und 2b, die an ihren vier Ecken punktgelagert ist. Die normierte Durchbiegung wird durch die Koordinatentransformationen  xi = x/a, eta = y/b;  auf ein Quadrat der Kantenlänge 2 transformiert und näherungsweise durch folgende Funktion beschrieben:
> w:=1/48*(10*(alpha^2+beta^2)+2-12*((alpha*xi)^2+(beta*eta)^2)+2*(xi^2+eta^2)+2*(alpha^2*xi^4+beta^2*eta^4)-6*xi^2*eta^2);
# Darin sind
> alpha:=1.2; beta:=1/alpha;
# Wir verschaffen uns ein dreidimensionales Bild dieser Funktion und verwenden für die grafische Darstellung den Befehl plot3d .
> plot3d(w,xi=-1..1,eta=-1..1,axes=boxed,orientation=[-115,-65,175]);
# Von Interesse ist die Durchbiegung mit der Darstellung der Konturlinien, das sind die Linien gleicher Durchbiegung w:
> plot3d(w,xi=-1..1,eta=-1..1,axes=boxed,style=patchcontour,contours=15,shading=zgrayscale,orientation=[-90,0,0],scaling=constrained);
# Zusätzlich eingezeichnet wird die nicht ausgebogene Platte (w = 0).
> plot3d([w,0],xi=-1..1,eta=-1..1,axes=boxed,orientation=[130,-60,0],lightmodel=light2,shading=zhue,transparency=.6);
# Die bezogene Durchbiegung w(ñ) einer dünnen kreisförmigen Platte unter Einzelkraft, die am Rand ñ = r/a = 1 eingespannt ist (a: positive Bezugslänge), wird beschrieben durch die Funktion  w = K(2*varrho^2*`ln&varrho;`-varrho^2+1);. An der Stelle ñ = 0 ersetzen wir die Funktion w durch ihren Grenzwert. 
# Zur Auswertung der Durchbiegung an diskreten Punkten (x,y) wird die  Prozedur Durchbiegung geschrieben.
> Durchbiegung:=proc(x,y,a,K)
> local rq,rhoq;
> global w;
> rq:=x^2+y^2;
> rhoq:= rq/a^2;
> if a<0 then print("Geometrie überprüfen");
> return end if;
> if rhoq =0 then w:=K;
> elif rhoq >1 then w:=0;
> else
> w:=K*(1-rhoq+2*rhoq*ln(sqrt(rhoq)));
> end if;
> end proc:
# Wir werten die Funktion am Punkt (x,y)  = ( 0.3, 0.3)  mit K = 1 und a = 2 aus und erhalten
> Durchbiegung(0.3,0.3,2,1);
# Wir machen uns ein Bild  dieser Funktion.
> plot3d((x,y)->Durchbiegung(x,y,1,1),-1..1,-1..1,axes=boxed,orientation=[-65,-55,175],scaling=constrained);
> plot3d((x,y)->Durchbiegung(x,y,1,1),-1..1,-1..1,axes=boxed,style=patchcontour,contours=10,shading=z,scaling=constrained,orientation=[-65,-65,175]);
# Im Zusammenhang mit der Berechnung von Grenzlasten für gerade Stäbe unter Eigengewichtsbelastung treten die Funktionen AiryAi(z) und AiryBi(z) auf, wobei z = x + Iy ist.  
> st:=time():  #Wir messen die dazu benötigte Zeit, das kann etwas dauern
;
> p1:=plot3d(Re(AiryAi(x+I*y)),x=-4..4,y=-4..4,view=-4..3,axes=boxed):
> p2:=plot3d(Im(AiryAi(x+I*y)),x=-4..4,y=-4..4,view=-4..3,axes=boxed):
> p3:=plot3d(Re(AiryBi(x+I*y)),x=-4..4,y=-4..4,view=-4..3,axes=boxed):
> p4:=plot3d(Im(AiryBi(x+I*y)),x=-4..4,y=-4..4,view=-4..3,axes=boxed):
> time()-st; #Ende der Zeitmessung
;
> plots[display](Matrix(1,4,[p1,p2,p3,p4]));
# Wir zeichnen noch die Höhenlinien (insgesamt 10) der Funktion Re(AiryAi)
> plot3d(Re(AiryAi(x+I*y)),x=-4..4, y=-4..4,view=-1..3,axes=boxed,style=patchcontour,contours=10,color=x,scaling=constrained,orientation=[-90,0,0]);
> restart:
# Maple gestattet das Arbeiten mit eingebetteten Objekten. Diese grafischen Komponenten ermöglichen einen interaktiven Zugriff auf den Maple-Code. Das Hinzufügen von grafischen Komponenten geschieht mit der  Palette Components. Wir zeigen das an einem einfachen Beispiel. Die Funktion y = tanh(ax); kann mit wachsendem Parameter a als stetige Approximation der unstetigen Signum-Funktion angesehen werden.  Wir wählen zur grafischen Darstellung des Problems aus der Palette Components 
# 1. einen Schieberegler (slider)
# 2. eine Komponente Plot und 
# 3. eine Box (Mathematical expression) , die uns den aktuellen Wert der Funktion  y = tanh(ax); anzeigen soll. Die Komponenten werden mittels drag and drop in einen Tabellenrahmen gezogen und dort wie gewünscht angeordnet. Ein rechter Mausklick auf den Schieberegler zeigt unter Component Properties dessen Werte. Sein Name ist Regler1. Der Schieberegler arbeitet im Wertebereich [0<=a<=10]. ;Unter Action When Value Changes wird das aus dem Regler kommende Signal der Variablen para1 zugeordnet, und die Funktionen signum(x)  sowie  tanh(ax) werden im Grafikfenster dargestellt. Durch Ziehen am Regler mit der Maus werden die Kurve  tanh(ax)  kontinuierlich dargestellt.
# Weitere Informationen und Beispiele können der Hilfefunktion  Embedded Components and Maplets  entnommen werden. 
> 
;
#              
#                                      
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
