# Lektion 3
# Grundlagen des Programmierens
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
# Zur Steuerung des Programmflusses stehen folgende Anweisungen und Operatoren zur Verfügung
# 1. Die if Anweisung
# 2. Der ` if ` Operator
# 3. Die for Schleife
# 4. Die while Schleife
> restart:
# Die Verschiebung  w(¾)  einer Balkenachse sei durch folgende Funktion beschrieben (0 d ¾ d 1):
> w:=unapply(4*xi^2-8*xi^3+5*xi^4-xi^5,xi);
> plot(w(xi),xi=0..1,color=blue,gridlines=true,axes=boxed,thickness=2);
# Um näherungsweise die maximale Durchbiegung zu ermitteln, tasten wir die Funktion w(¾), bei Null beginnend, mit einer Schrittweite ”¾  ab.  Die Abtastwerte werden in die Liste abtast  geschrieben. Die Aufgabe besteht im Folgenden darin, die größte Zahl aus dieser Liste zu ermitteln. Wir teilen den Definitionsbereich 0 d ¾ d 1 in  n  gleichlange Intervalle  der Länge  ”¾ = 1/n . Damit ergeben sich genau n + 1 Abtastpunkte. Wir wählen n = 5.
> n:=5; dxi:=1./n;
> abtast:=[seq(w(i),i=0..1,dxi)];
# Der erste Kandidat für das gesuchte Maximum ist 
> v:=abtast[1];
# Es werden nun nacheinander alle übrigen Elemente abgefragt. Wenn ein Element größer als der bisherige Kandidat ist, wird es zum neuen Kandidaten. Der letzte Kandidat repräsentiert schließlich das gesuchte Maximum. Die if-Anweisung ermöglicht es, einen logischen Ausdruck auszuwerten und in Abhängigkeit vom Ergebnis bestimmte Anweisungen auszuführen oder zu überspringen. Sie heißt darum auch bedingte Anweisung. Ist die in der if-Anweisung genannte Bedingung erfüllt, wird die Anweisung hinter dem Schlüsselwort  then  ausgeführt, anderenfalls wird der Befehl nach end if  ausgeführt.
> if abtast[2] > v then v:=abtast[2] end if;
> if abtast[3] > v then v:=abtast[3] end if; 
> if abtast[4] > v then v:=abtast[4] end if;
> if abtast[5] > v then v:=abtast[5] end if; 
> if abtast[6] > v then v:=abtast[6] end if; 
# Bei einer größeren Anzahl von Listenelementen kann diese Vorgehensweise sehr mühsam werden. Effektiver ist die Einführung einer bedingten Schleife, die auch do-Schleife  genannt wird. Die Anzahl der Listenelemente wird durch den Befehl nops ermittelt.
> v:=abtast[1]:
> for j from 2 to nops(abtast) do
>    if  abtast[j] > v then v:=abtast[j]  end if
> end do:
> v;
# In der folgenden Lösung durchläuft  die Variable ww alle Werte der Liste abtast, und es kann auf den Schleifenindex gänzlich verzichtet werden. 
> v:=abtast[1]:
> for ww in abtast do if ww > v then v:= ww  end if  end do:
> v;
# Neben der if-Anweisung stellt Maple den if-Operator  zur Verfügung. Er verlangt drei Argumente. Das erste Argument (im Folgebeispiel ist das der Wert 3) wird auf die Wahrheitsgehalte  true, false und FAIL geprüft. Ist der Wahrheitsgehalt true, dann wird das zweite Argument zurückgegeben. Ist der Wahrheitsgehalt dagegen  false oder FAIL, dann wird das dritte Argument zurückgegeben. In Boolschen Ausdrücken bedeutet FAIL einen unbekannten oder unbestimmten Wert.
> c:=4; `if`(c > 3,t,f); `if`(c < 3,t,f);
# Soll eine Schleife so lange durchlaufen werden, bis eine Bedingung erfüllt ist, so tritt  while .. do  an die Stelle von  for .. do. 
# Hierzu ein Beispiel zur Berechnung der Quadratwurzel einer positiven Zahl a, also x = sqrt(a);. Das sind die Nullstellen der Funktion f(x) = x^2-a; .  Nach Newton gilt die Iteratonsvorschrift:  x[i+1]=  x[i] - (f(x[i]))/(f '(x[i])); (i = 1...n). Die Ableitung der Funktion  f  nach x  ist f '(x) = 2*x;. Damit folgt die Iterationsvorschrift: x[i+1] = 1/2*(x[i]+a/x[i]);  mit x[i] <> 0;. Das Verfahren erwartet eine  (möglichst gute) Startnäherung xstart. Die Iteration wird abgebrochen, wenn der absolute relative  Fehler  frel = |(xi+1 - xi)/xi|  kleiner ist als die vorgegebene Genauigkeitsschranke genau. Gesucht werden der Näherungswert für sqrt(a);, die  Anzahl n der dazu benötigten Iterationen und der sich dabei einstellende absolute relative Fehler frel.
> f:=unapply(x^2-a,a,x);
# Wir wollen die Quadratwurzel aus a = 10 berechnen. Es werden nur positive Werte für a zugelassen. 
> a:=10.: if a <= 0 then print("Der Radikand muss positiv sein") end if: 
# Wir wählen mit xstart = 1 einen recht groben Startwert 
> xstart:=1.: xneu:=xstart: frel:=1.: n:=0: genau:=0.1e-6:
> while frel >= genau do
>   n   :=n+1;
>   xalt:=xneu;
>   xneu:=0.5*(xalt+a/xalt);
>   frel:=abs(xneu/xalt-1);
> end do:
# Wir geben das Ergebnis aus. Der letzte Wert von xneu entspricht der Näherung sqrt(a) , ;n bezeichnet die Anzahl der Iterationen und frel ist der absolute relative Fehler.
> xneu, n, frel;
# Hat der Ausdruck nach while den Wahrheitswert true, werden die Befehle in der Schleife ausgeführt. Das wird solange fortgesetzt, bis der Ausdruck den Wert false oder FAIL besitzt. Wird hinter end do ein Doppelpunkt gesetzt, dann erfolgt keine Ausgabe der innerhalb der Schleife berechneten Größen. Soll eine Ausgabe erfolgen, dann muss ein Semikolon gesetzt werden. 
#  Maple Prozeduren
# Maple besitzt eine Fülle von eingebauten Prozeduren, etwa min und max zum Aufsuchen des Minimums bzw. des Maximums einer Sequenz. Wir wenden diese Prozeduren auf die Liste abtast an, die zunächst mit dem Befehl op in eine Sequenz umgewandelt werden muss, da im Argument dieser Prozeduren Sequenzen erwartet werden.
# Hinweis: Die Berechnung des Maximums einer Sequenz ist als Beispiel unten auf der Hilfeseite zu  nargs  beschrieben.
> min(op(abtast)); max(op(abtast));
# Um die Wirkungsweise solcher Prozeduren zu verstehen, wird eine eigene Prozedur zur Berechnung der Quadratwurzel aus einer positiven Zahl  geschrieben. Eine Maple Prozedur-Definition hat die folgende allgemeine Syntax:
# proc ( P )
#    local L;
#    global G;
#    options O;
#    desription D;
#    B
# end proc:
# Es bedeuten
#   P:	 Formalparameter
#   L: 	Lokale Variable (optional)
#  G: 	Globale Variable (optional)
#  O: 	Optionen (optional)
#  D: 	Beschreibungen
#  B: 	Sequenz von Befehlen, die den Rumpf des Programms bilden
# Bei jeder Prozedur (procedure) ist zwischen Definition und Aufruf zu unterscheiden. Wir geben unserer Prozedur den Namen qwurz. Als Formalparameter übergeben wir die Werte a (Radikand), xstart (Startnäherung), nmax (größte Anzahl Iterationen) und eps (absoluter relativer Fehler).
> qwurz:=proc(a,xstart,nmax,eps)
> local fr,n,xalt,xneu;
> description "Berechnung der Quadratwurzel aus einer positiven Zahl";
> fr:=1.: n:=0: xneu:=xstart:;
> if a <= 0 then print(`Der Radikand muss positiv sein`);
>   return;
> end if;
> while fr >= eps do
>   n:=n+1;
>   if n > nmax then
>     print(`Die Anzahl der zulässigen Iterationen wurde überschritten`); 
>     return; 
>   end if;
>   xalt:=xneu;
>   xneu:=0.5*(xalt + a/xalt);
>   fr:=abs(xneu/xalt-1.);
> end do:
> print(xneu,n,fr);
> end proc:
# Aufruf der Prozedur mit folgenden Aktualparametern:
> r:=10.: start:=1.: nm:=20: genau:=0.1e-6:
> qwurz(r,start,nm,genau); 
# Die Werte  r (Radikand) , start (Startwert für die Lösung), nm (Maximalanzahl der Iterationen) und die vorgegebene Genauigkeit genau stellen beim Aufruf die  aktuellen Parameter der Prozedur dar. Als Platzhalter werden in der Prozedur die Namen a, xstart, nmax und eps benutzt, die erst beim Aufruf konkretisiert werden. Sie werden auch formale Parameter  der Prozedur (dummy variable) genannt. Die zweite Zeile definiert die lokalen Variablen. Ist eine Variable nicht deklariert, dann meldet sich Maple mit der Warnung: 
#   Warning , `name` is implicitly declared local to procedure `name`
# Die lokalen Variablen verlieren  ihren Wert nach Abarbeitung der Prozedur. Sollen diese jedoch später noch zugänglich sein, dann müssen sie  innerhalb der Prozedur als global deklariert sein. 
# Hinweis: Es empfiehlt sich, aus Gründen der Übersichtlichkeit, die Deklaration aller  Variablen einer Prozedur. Weitere Einzelheiten s. procedure und  parameter. 
# Die description Sequenz kann optional zur Beschreibung der Prozedur eingesetzt  und mit dem Befehl Desribe angezeigt werden.
> Describe(qwurz);
# Ein return Befehl verusacht den  Rücksprung zu dem Punkt, von dem die Prozedur aufgerufen wurde. Das erfolgt beispielsweise bei der Eingabe einer negativen Zahl und bei der Überschreitung der Anzahl der zulässigen Iterationen (hier 20). Wir rufen die Prozedur mit einem negativen Wert auf
> qwurz(-r,start,nm,genau); 
# Um die Definitionen einer Prozedur anzusehen, reicht es nicht, den Prozedurnamen einzugeben. Dazu benötigen wir den print-Befehl
> qwurz;
> print(qwurz);
# Mit dem Kommando showstat erhalten wir eine Darstellung des Programms, etwa im Zusammenhang mit einer Fehlersuche, in der jede neue Kommandozeile durchnummeriert ist. 
> showstat(qwurz);
# Eine Methode, die Ausführung einer Prozedur an einer bestimmten Stelle zu unterbrechen,  besteht u.a. darin, einen  breakpoint  zusetzen. Breakpoints werden durch Sterne nach der Zeilennummer gekennzeichnet. Das Setzen der breakpoints geschieht mit dem Kommando stopat . In unserem Beispiel werden breakpoints vor der 1. und 4. Zeile gesetzt. Maple stoppt dann vor Ausführung der 1. Zeile. Mit den  Befehlen  stopat, stopwhen und stoperror wird automatisch der Interactive Maple Debugger geöffnet.
> stopat(qwurz,1),stopat(qwurz,4);
> showstat(qwurz);
# Wir rufen die Prozedur qwurz auf und haben die Möglichkeit, das Programm im Debugg-Modus zu überprüfen. Das geschieht mit den Kommandos: next, step, into, list, outform, cont. 
> qwurz(10.,1.,20,0.1e-6):
# Um die mit stopat gesetzten breakpoints zu enfernen, wird das Kommando unstopat gegeben.
> unstopat(qwurz): showstat(qwurz);
> restart:
# Wir schreiben eine weitere Prozedur, die nach dem einfachen Newton-Verfahren die Nullstellen des folgenden Polynoms  f  berechnen soll.
> f:=proc(x) x^2-2*x-1 end;
> newton:=proc(xstart,nmax,eps)
> local fr,n,xalt,xneu,N;
> global f;
> description "Berechnung der Quadratwurzel aus einer positiven Zahl";
> fr:=1.: n:=0: xneu:=xstart:
> while fr >= eps do
>   n:=n+1;
>   if n > nmax then
>     print(`Die Anzahl der zulässigen Iterationen wurde überschritten`); 
>     return; 
>   end if;
>   xalt:=xneu;
>   N   := D(f)(xalt);                                
>   xneu:=xalt - f(xalt)/N;
>   fr  :=abs(xneu/xalt-1.);
> end do:
> print(evalf(xneu),n,fr);
> end proc:
# Wir rufen die Prozedur mit dem Startwert xstart = 1 auf. (An dieser Stelle verschwindet die 1. Ableitung von f, was innerhalb der Prozedur zu einem Problem führt!)
> newton(1,20,0.1e-6);
# Die Prozedur läuft auf einen Fehler (Division durch Null). Um dem Fehler auf die Spur zu kommen,  nummerieren wir zunächst die Kommandozeilen
> showstat(newton);
# Divisionen, und nur dabei kann der Fehler auftreten, finden wir in den Zeilen 11 und 12. Dort wird durch N und  xalt  dividiert. Wir setzen mit dem Befehl stopwhen  Beobachtungspunkte vor  die betreffenden Zeilen. 
> stopwhen([newton,N]): stopwhen([newton,xalt]):
# Wir starten die Prozedur nochmal und das Debugger-Fenster öffnet sich:
> newton(1,20,0.1e-6);
# Im Debugger-Mode sehen wir, dass die Variable N  im Laufe der Berechnung den Wert Null annimt. Weitere Informationen erhalten wir mit dem Befehl tracelast , der uns die Werte der Variablen vor Abbruch des Programms anzeigt.
> tracelast;
# Wir entfernen die Beobachtungspunkte aus der Prozedur
> unstopwhen():
# und  rufen die Prozedur abermals auf, jetzt aber mit dem Startwert xstart = 0
> newton(0,20,0.1e-6);
# Auch in diesem Fall  erhalten wir eine Fehlermeldung. 
> tracelast;
# In diesem Fall ist xalt = 0. Das Programm ist also noch sehr verbesserungswürdig.  
# Um verbesserte Startwerte zu finden, geben wir die Funktion f und deren 1. Ableitung  grafisch aus:
> plot([f,D(f)],-1..3,gridlines=true,color=[black,blue]);
# und wählen jetzt die verbesserten Startwerte x1 = 2.5 und x2 = -0.5, die zu  konvergierenden Lösungen führen :
> newton( 2.5,20,0.1e-6);
> newton(-0.5,20,0.1e-6);
# In einigen Fällen können wir uns den Maple- Quellcode ansehen, so beispielsweise für den natürlichen Logarithmus ln und die Exponentialfunktion exp:
> interface(verboseproc=2);
> eval(ln);
> 
;
> 
;
