# Lektion 8
# Differenzialrechnung und Taylor-Reihen
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart:with(LinearAlgebra): with(plots):
# Die Ableitung von Ausdrücken oder Operatoren erfolgt in Maple mit dem  diff-Kommando  oder dem  D-Operator. Das Diff-Kommando ist untätig und schreibt nur den Befehl, die Ausführung selbst erfolgt mit diff. Beispiele hierzu sind:
> Diff(sin(x),x)=diff(sin(x),x);
> Diff(x*sin(cos(x)),x)=diff(x*sin(cos(x)),x);
# Höhere Ableitungen
> Diff(x*sin(cos(x)),x,x,x)=diff(x*sin(cos(x)),x,x,x);
# Das geht auch etwas kürzer mit 
> Diff(x*sin(cos(x)),x$3)=diff(x*sin(cos(x)),x$3);
# Ebenfalls können partielle Ableitungen gebildet werden 
> Diff(x^2+x*y^2,x)=diff(x^2+x*y^2,x);
> Diff(2*x+y^2,y)=diff(2*x+y^2,y);
> Diff(x^2+x*y^2,y,x)=diff(x^2+x*y^2,y,x);
> diff(f(x,y),y,x);
# Ableitung transzendenter Funktionen:
> Diff(tan(x),x)=diff(tan(x),x);
> Diff(BesselJ(v,x),x)=diff(BesselJ(v,x),x);
# Die Ableitungen von Operatoren werden mit dem D-Operator durchgeführt. Zum Unterschied zwischen D und diff  wird auf  operators[D] verwiesen.
> D(sin);D(sqrt);D(ln);
> D(tan);%(x);D(tan)(x);
> g:=unapply(x^4,x);
# Bildung der 1. Ableitung mit dem Operator D
> D(g);
> D(g)(x);D(g)(2);
# Ohne auf die Begriffe Geschwindigkeit und Beschleunigung näher einzugehen, wird als Beispiel der schiefe Wurf ohne Luftwiderstand betrachtet.  In den folgenden Gleichungen bezeichnet t e 0 die Zeit, v0 ~ 0  die positive Anfangsgeschwindigkeit und 0 d ± d À/2 den Abwurfwinkel einer Punktmasse gegenüber der positiven x-Achse. Zur Zeit t  hat die Punktmasse die Koordinaten
> x:=unapply(v0*t*cos(alpha),t,v0,alpha); y:=unapply(-g*t^2/2+v0*t*sin(alpha),t,v0,alpha);
#  Eliminieren wir aus der zweiten Gleichung für y die aus der ersten Gleichung für x folgende Zeit t, also  t = x /(v0 cos±) , dann erhalten wir folgende explizite Darstellung der Wurfparabel
> ye:=unapply(subs(t=x/(v0*cos(alpha)),y(t,v0,alpha)),x,v0,alpha);
# Die Koordinate xS des Scheitelpunktes der Wurfparabel erhalten wir aus der Forderung
> xS:=solve(D[1](ye)(x,v0,alpha),x);
# Die maximale Wurfhöhe 
> yS:=ye(xS,v0,alpha);
# wird damit zu folgendem Zeitpunkt erreicht:
> tS:=xS/v0/cos(alpha);
# Die größte Wurfweite  ermitteln wir aus der Bedingung  ye = 0. Wir erhalten zwei Lösungen, wobei wir hier nur an xmax > 0 interessiert sind:
> solve(ye(x,v0,alpha),x);
> xmax:=%[2];
# Bei festem  Zielpunkt xmax  und fester Anfangsgeschwindigkeit  v0  existieren zwei Abwurfwinkel ±: 
> alpham:=solve(diff(xmax,alpha),alpha);
# Wir interessieren uns für den Wurf nach rechts (± < À/2) und wählen deshalb
> alphamax:=alpham[1];
# Die größte Wurfweite wird also bei einem Abwurfwinkel von ± = À/4 (± = 45°) erreicht, und die Wurfzeit beträgt 
> solve(subs(alpha=alphamax,y(t,v0,alpha)),t); tmax:=%[2];
# Wir beschaffen uns eine Animation der Wurfbahnen für unterschiedliche Abwurfwinkel ±
> g:=10: v0:=15: alpha1:=Pi/6: alpha2:=Pi/4: alpha3:=Pi/3: xm:=eval(subs(alpha=alphamax,xmax)): ym:=subs(alpha=Pi/3,yS):
> pl1:=animatecurve(ye(x,v0,alpha1),x=0..xm,frames=20,color=black):   # alpha=30°
> pl2:=animatecurve(ye(x,v0,alpha2),x=0..xm,frames=20,color= gold):   # alpha=45°
> pl3:=animatecurve(ye(x,v0,alpha3),x=0..xm,frames=20,color=  red):   # alpha=60°
;
> display([pl1,pl2,pl3],axes=boxed,gridlines,view=0..ym,title="\nWurfparabeln für unterschiedliche Abwurfwinkel\n",titlefont=["ROMAN", 15],labels=["Wurfweite","Wurfhöhe"],labeldirections=["horizontal", "vertical"], labelfont = ["HELVETICA", 10]);
# Die Abwurfwinkel ±1 = À/6 (Flachwurf)  und ±3 = À/3 (Steilwurf)  führen zu gleichen Wurfweiten xmax = v*`0`^2*sqrt(3)/(2*g); .
# Als weiteres Beispiel betrachten wir die Biegelinie eines Balkens der Länge Á, der links eingespannt und rechts drehbar gelagert ist. Im Einzelnen bedeuten (q0: Konstante Linienkraft, : Länge des Balkens, E: Elastizitätsmodul, Iyy: Axiales Flächenträgheitsmoment, B = EIyy: Biegesteifigkeit,  ¾ = x/Á: Normierte Stabachskoordinate) 
> K:=q0*l^4/(48*B);
> w:=unapply(K*(3*xi^2-5*xi^3+2*xi^4),xi);
# Aus der Biegelinie lassen sich durch Differenziationsoperationen folgen Größen berechnen:
# 1.)  w'                  Neigung der Tangente an die Biegelinie
# 2.)  M = - Bw''    Biegemoment 
# 3.)  Q  = - Bw'''   Querkraft 
# Achtung: In den folgenden Ableitungen ist die Kettenregel diff(w(x), x) = (diff(w(xi), xi))*(diff(xi(x), x)) and (diff(w(xi), xi))*(diff(xi(x), x)) = (diff(w(xi), xi))/`&ell;`; zu beachten.
> wst:=unapply(diff(w(xi),xi)/l,xi);
> M:=unapply(-B*diff(w(xi),xi$2)/l^2,xi);
> Q:=unapply(-B*diff(w(xi),xi$3)/l^3,xi);
# Die Durchbiegung wird dort extremal, wo deren 1. Ableitung verschwindet
> loe:=solve(wst(xi),xi);evalf(loe);
# Der physikalisch sinnvolle Wert im Bereich 0 < ¾ < 1 ist
> xiw:=loe[3];
# Dazu gehört die größte Durchbiegung im Feld
> maxw:=simplify(w(xiw));evalf(%);
# Das Schnittlasten (Biegemoment und Querkraft) am eingespannten Rand (x = 0) sind
> ME:=M(0);QE:=Q(0);
# Von Interesse ist noch das größte Feldmoment. Es wird dort extremal, wo die Querkraft verschwindet. 
> xst:=solve(Q(xi),xi);MFm:=simplify(M(xst));
# Wir konkretisieren unser Problem und plotten die normierten Zustandsgrößen für folgendes Zahlenbeispiel
> q0:=1;B:=1;l:=1;
> par:=xi=0..1,titlefont=["ROMAN",14],color="NavyBlue",thickness=3,filled=[color = "Blue",transparency = .5],gridlines,axes=boxed:
> p1:=plot(w(xi),par,title="\nNormierte Durchbiegung\n"):
> p2:=plot(wst(xi),par,title="\nNormierte Tangentenneigung\n"):
> p3:=plot(M(xi),par,title="\nNormiertes Biegemoment\n"):
> p4:=plot(Q(xi),par,title="\nNormierte Querkraft\n"):
> plotmatrix:=Matrix(1,4,[p1,p2,p3,p4]):
> plots[display](plotmatrix);
# Die größte (normierte) Durchbiegung im Feld ist
> maxw:=simplify(w(xiw));evalf(maxw);
# maxw ist eine irrationale Zahl, die durch eine Folge rationaler Zahlen dargestellt werden kann und gegen den Grenzwert maxw konvergiert. Ein solche konvergente Folge bildet die Kettenbruchentwicklung von maxw, die in Maple mit dem Befehl cfrac  aus der Bibliothek  numtheory  erzeugt wird:
> numtheory[cfrac](maxw,6);
# oder auch weniger umfangreich mit Verwendung einer Zusatzoption confrac
> convert(maxw,confrac,liste);
# Wir geben die Liste der Kettenbrüche aus
> liste;
> evalf(%);
# In Tabellenwerken der Ingenieurliteratur findet man den Näherungswert maxw = 2/369 q0 Á4/B. Wir lassen uns die obige Liste grafisch ausgeben und stellen fest, dass sich die Werte  für n > 1 nur wenig unterscheiden, dennoch ist der Vorfaktor 1/185 den Ingenieuren noch zu ungenau.
> p1:=listplot(liste,style=point,color=red,symbolsize=10): display(p1);

# Maple gestattet uns mit dem Befehl  map  das Differenzieren von Matrix- und Vektorkomponenten. Wir konstruieren dazu eine Matrix A:
> f := unapply((x-y)^(i-j),i,j);
> A:=Matrix(3, f): Ax:=map(diff,A,x): Axy:=map(diff,A,x,y):
> print(A,Ax,Axy);
# Das funktioniert  auch mit Vektoren. Wir beschaffen  uns dazu die dritte Spalte von A:
> V3_A := Column(A,3): 
# und differenzieren partiell
> V3_Ax:=map(diff,V3_A,x): V3_Ay:=map(diff,V3_A,y): V3_Axy:=map(diff,V3_A,x,y):
> print(V3_A, V3_Ax, V3_Ay, V3_Axy);
# Taylor-Reihen
# Die Taylor-Reihe ist eine Potenzreihe, deren Glieder Potenzen und damit Funktionen einer unabhängigen Variablen (hier x) sind.
# P(x) = sum(f(x[0])^n*(x-x[0])^n/factorial(n), n = 0 .. infinity);
# Der Punkt x0 heißt Entwicklungspunkt, und  ((D@@n)(f))(x[0]); bezeichnet die n-te Ableitung der Funktion f(x) an der Stelle x = x[0];. Wird die Reihe nach dem n-ten Glied abgebrochen, dann sprechen wir von einem Taylor-Polynom der Ordnung n und erhalten
# P[n](x) = sum(f(x-x[0])^nu*(x-x[0])^nu/factorial(n), nu = 0 .. n);
# Die Konstruktion des Taylor-Polynoms  P[n](x); erfolgt so, dass an der Stelle x[0]; die Taylor-Reihe denselben Funktionswert und dieselben Ableitungen bis zur Ordnung n besitzt wie die Funktion f(x) . 
# Als Beispiel betrachten wir den Hyperbelsinus 
> f:=sinh(x);
# Die folgenden Befehle liefern uns die Taylor-Polynome ersten  bzw. vierten Grades, entwickelt an der  Stelle x0 = 1. 
> x0:=1:
> n:=1: p1:=convert(series(f,x=x0,n+1),polynom); 
> n:=4: p4:=convert(series(f,x=x0,n+1),polynom); 
> plot([f,p1,p4],x=-3..3,-4..5,color=[black,blue,red],legend=["sinh","p1","p4"],axes=normal,thickness=2,gridlines);
# Die blaue Linie zeigt die Tangente an der Stelle x = 1. Das Polynom höheren Grades (n = 4) führt im positiven Bereich offensichtlich zu einer besseren Anpassung. 
# Wir betrachten die Funktion
> f:=1/(1-x); plot(f,x=-2..4,-5..5,axes=boxed,gridlines);
# und entwickeln sie in eine Potenzreihe an der Stelle x = 0
> series(f,x=0,6);
# Die obige Reihe wird geometrische Reihe genannt. Sie konvergiert nur für Werte |x| < 1 und hat dann den Summenwert 1/(1-x);
> n:='n':
# Das folgende Beispiel zeigt die Potenzreihenentwicklung der binomischen Reihe um den Punkt x[0] = 0;. 
> f:=(1+x)^n;       series(f,x=0,5);
# Ist n eine nichtnegative ganze Zahl, dann bricht die Reihe mit der höchsten Potenz xn ab:
> n:=3: series(f,x=0,5);
> n:=1/2:series(f,x=0,5);
# Reihenentwicklungen von Funktionen werden mit großem Erfolg im gesamten Ingenieurwesen eingesetzt. Wollen wir beispielsweise Werte der Kosinusfunktion im Intervall [0, À/4] näherungsweise berechnen, dann beschaffen wir uns zunächst die entsprechende Reihenentwicklung, indem wir cos(x) um den Punkt x[0] = 0; in eine Reihe entwickeln. 
> series(cos(x),x=0,10);
# Wir verwenden folgende Näherungspolynome
> TP_6:=convert(series(cos(x),x=0,8),polynom);
> TP_4:=convert(series(cos(x),x=0,6),polynom);
# und stellen die Funktion cos(x) und ihre Näherungen im Intervall 0  d x d À grafisch dar
> plot([cos(x),TP_6,TP_4],x=0..Pi,-1..1,color=[black,blue,red],legend=["cos","TP_6","TP_4"],axes=boxed,thickness=2,gridlines);
# Ein Blick auf die Grafen zeigt, dass für Werte 0 < x < À/2 die Kosinusfunktion durch beide Näherungen recht gut wiedergegeben wird. Aufschlußreich ist auch eine Darstellung des wahren relativen Fehlers für die Approximation TP_6 im Intervall [0, À/4]:
> plot((TP_6-cos(x))/cos(x),x=0..Pi/4,axes=boxed,gridlines);
# Für die Wahl x = À/4 vergleichen wir den Näherungswert mit dem exakten Wert cos(À/4) =  1/2 sqrt(2); und finden Übereinstimmung auf fünf Stellen hinter dem Komma:
> evalf(subs(x=Pi/4,TP_6)),1/sqrt(2.);
# Schlechter wird das Ergebnis, wenn wir den Kosinus des Winkels À/2 mit der Näherung TP_6 berechnen, der exakte Wert ist Null.
> evalf(subs(x=Pi/2,TP_6));
# Partielle Ableitungen
# Wir betrachten die normierte Biegefläche einer durch eine konstante Flächenkraft beanspruchte dünne Platte im Intervall [-1 <= xi and xi <= 1, -1 <= eta and eta <= 1];, die an ihren Ecken punktförmig gelagert ist.
> w:=unapply(11/24-(5/24)*xi^2-(5/24)*eta^2+(1/24)*xi^4+(1/24)*eta^4-(1/8)*xi^2*eta^2,xi,eta);
> plot3d(w(xi,eta),xi=-1..1,eta=-1..1,axes=boxed,labels=[xi,eta,w],grid=[60, 60],orientation=[125,-60,180],color=xi,lightmodel=light1);
# Von dieser Funktion bilden wir die partiellen Ableitungen nach ¾ und ·:
> Diff(w(xi,eta),xi) = D[1](w)(xi,eta);
> Diff(w(xi,eta),eta) = D[2](w)(xi,eta);
# Nun bilden wir noch die drei zweiten partiellen Ableitungen. 
> Diff(w,xi  ,xi)=D[1,1](w)(xi,eta);
> Diff(w,eta,eta)=D[2,2](w)(xi,eta);
> Diff(w,xi ,eta)=D[1,2](w)(xi,eta);
> Diff(w,eta ,xi)=D[2,1](w)(xi,eta);
# Maple bestätigt mit den letzten zwei Zeilen die Schwarzsche Vertauschungsregel. 
# Die partiellen Ableitungen der Biegefläche hängen eng mit den Schnittlasten der Platte zusammen. Im Fall verschwindender Querdehnung ist beispielsweise -bis auf eine Konstante- das Biegemoment m[`&xi;&xi;`];  gegeben durch
# m[xixi]=-((&PartialD;)^(2) w(xi,eta))/(&PartialD; xi^(2)) ;
> mxixi:=-D[1,1](w)(xi,eta);  
> plot3d(mxixi,xi=-1..1,eta=-1..1,style=Patch,axes=boxed,color=xi,lightmodel=light1,scaling=constrained,orientation=[-135,55,0]);
# Wir lassen uns die Konturlinien anzeigen, das sind Linien gleicher Funktionswerte, die in der Statik zur nummerischen Auswertung von Einflussflächen benötigt werden.
> plot3d(mxixi,xi=-1..1,eta=-1..1,axes=boxed,style=patchcontour,contours=10,scaling=constrained,orientation=[-90,0,0]);
# Zur Bildung der Ableitung von Funktionsoperatoren kommt der D -Operator zum Einsatz. Damit stellen sich die ersten und zweiten Ableitungen wie folgt dar:
> D[1](w); D[2](w);
> D[1,1](w); D[1,2](w); D[2,2](w);
# Wir ermitteln die Gleichung der Tangentialebene an die Fläche w(¾, ·) im Punkt (¾0,·0) , die wir aus der zweidimensionalen Taylor-Reihe vom Grade 1 erhalten. Dazu benötigen wir das Paket MultivariateCalculus  aus der Bibliothek  Student
> with(Student[MultivariateCalculus]):
> TF1:=unapply(TaylorApproximation(w(xi,eta),[xi,eta]=[xi0,eta0],1),xi,eta);
# Der folgende Befehl liefert die zweidimensionale quadratische Approximation der Funktion w(¾, ·) im Punkt (¾0,·0) 
> TF2:=unapply(TaylorApproximation(w(xi,eta),[xi,eta]=[xi0,eta0],2),xi,eta);
# Wir stellen die Fläche  w(¾, ·)  und deren Tangentialebene für den Entwicklungspunkt P = (¾0 = 0.5, ·0= 0.)  in einer Grafik dar.
> xi0:=0.5: eta0:=0.: 
> plot3d([w(xi,eta),TF1(xi,eta)],xi=-1..1,eta=-1..1,axes=framed,color=[blue,green],style=[patchcontour,patchnogrid],lightmodel=light1,transparency=[0,.6],orientation=[-120,-60,-180]);
# Die Approximation der Funktion durch ein Taylor-Polynom zweiten Grades, also die Annäherung der Fläche durch eine Fläche zweiten Grades, erfordert Übereinstimmungen nicht nur in den nullten und ersten, sondern auch in den zweiten partiellen Ableitungen. Wir stellen die Fläche w zusammen mit ihrer Näherung zweiten Grades grafisch dar.
> plot3d([w(xi,eta),TF2(xi,eta)],xi=-1..1,eta=-1..1,style=[patchcontour,patchnogrid],axes=framed,color=[blue,red],lightmodel=light1,transparency=[0,0.6],orientation=[-50,-85]);
# Man sieht, dass sich die Flächennäherung zweiten Grades in der Umgebung des Entwicklungspunktes bereits gut an die Ursprungsfläche anschmiegt. Um den Approximationsfehler genauer beurteilen zu können, plotten wir den Betrag der Differenz beider Funktionen aus.
> plot3d(abs(w(xi,eta)-TF2(xi,eta)),xi=-1..1,eta=-1..1,style=patchcontour,axes=boxed,shading=z,lightmodel=light1,grid=[50,50],orientation=[-65,60,0]);
> restart:
# Der Zusammenhang zwischen zwei Variablen x und y  kann in folgenden Formen gegeben werden:
# 1. Explizite Funktion:  y = f (x)                      Beispiel: y = 1/x
# 2. Implizite Funktion: F(x,y) = 0                    Beispiel: x2 + y2 - 1 = 0
# 
# Wir interessieren uns für die Ableitung einer in impliziter Form gegebenen Funktion F(x,y) = 0.  Als Beispiel wählen wir die als Kartesisches Blatt bezeichnete Funktion
> F:=unapply(x^3+y^3=3*a*x*y,x,y,a);
#  Die Gerade 
> g:=unapply(x+y+a,x,y,a);
# ist Asymptote dieser  Kurve. Wir wählen a = 2 und geben den Grafen der Kurve aus:
> plots[implicitplot]([F(x,y,2),g(x,y,2)],x=-5..4,y=-5..4,linestyle=[solid,longdash],coloring=["SteelBlue","DarkGreen"],filledregions=true,gridrefine=2,crossingrefine=2,resolution=1000);
# Sei y = y(x), dann erhalten wir aus  (&PartialD;)/(&PartialD; x) F(x,y) = 0=(&PartialD;)/(&PartialD; x) (x^(3)+y^(3)-3 a x y)=3 x^(2)+3 y^(2)y'-3 ay-3 axy'.;
# Lösen wir diese Beziehung nach  (diff(y(x), x))*auf, dann;erhalten wir diff(y(x), x); =(-x^2+ay)/(y^2-ax); .
# Dieses Ergebnis liefert uns auch Maple mit dem Befehl  implicitdiff  
> ystr:=unapply(implicitdiff(F(x,y,a),y,x),x,y,a);
# Mit  a = 2 liegt der Punkt  P0  (x0 = 3, y0 = 3) auf der Kurve:
> F(3,3,2);
# Dort besitzt die Kurve die Tangentenneigung
> ystr(3,3,2);
> 
;
# 
# 
# 
