> restart:
# Beispiele zu Prozeduren
# 
# 1) Allgemeine numerische Integration
# Gesucht ist naeherungsweise die Flaeche unter dem Graphen der Funktion  f .
# Wir teilen die Flaeche in  n  Streifen von gleicher Breite  h .
# Wenn wir die Kurve ueber jedem Streifen durch ihre Tangente in der
# Mitte, also an der Stelle  xm, ersetzen, dann ist der Flaecheninhalt jedes Streifens  h*f(xm).
# Zu pruefende Typen:
# procedure: Funktion,  realcons:  reelle Konstante,   posint: positive ganze Zahl (=natuerliche Zahl)
> numint:=proc(f::procedure,xu::realcons,xo::realcons,n::posint)
> local  xudez,xodez,h,xm,flaeche;
> description "Numerische Integration der Funktion f von xu bis xo mit n
> Streifen";
> xudez:=evalf(xu);
> xodez:=evalf(xo);
> h:=(xodez-xudez)/n;
> xm:=xudez+h/2;
> flaeche:=0;
> to n do
>     flaeche:=flaeche +h*f(xm);
>     xm:=xm+h
> end do;
> flaeche
> end proc:
# Integration einer eingebauten Funktion 
> numint(sin,0,Pi/2,1);
> numint(sin,0,Pi/2,10);
> numint(sin,0,Pi/2,1000);
> int(sin(x),x=0..Pi/2); # Zum Vergleich exakt
# Integration einer selbst definierten Funktion
> g:=proc(x) 3*x^2 end;
> numint(g,0,1,10);
> numint(x->3*x^2,0,1,10); # Das ist einfacher
> int(3*x^2,x=0..1); # Zum Vergleich exakt
> numint(3*x^2,0,1,10);
> numint(x->3*x^2,0,u,10);
> numint(x->3*x^2,0,1,1.5);
> numint(x->3*x^2,0,1,-3);
> numint(x->1/x,1,5,10);
> int(1/x,x=1.0..5.0); # Zum Vergleich
> int(1/x,x=1..5); # Zum Vergleich exakt
# Man sieht, dass Logarithmen sich als Integrale berechnen lassen.
# Auf dieser Grundlage schreiben wir jetzt eine Prozedur, die
# Logarithmen naeherungsweise berechnet:

# 2) Logarithmus
> restart:
> logar:=proc(x::realcons,n::posint)
> local xdez,h,xm,flaeche;
> description "Berechnung des natuerlichen Logarithmus von x durch
> numerische Integration mit n Streifen";
> xdez:=evalf(x);
> if xdez<=0 then error "Wert nicht reell" end if;
> h:=(xdez-1.)/n;
> xm:=1.+h/2;
> flaeche:=0;
> to n do
>     flaeche:=flaeche + h/xm;
>     xm:=xm+h
> end do;
> flaeche
> end proc:
> 
> logar(5,1);
> logar(5,100);
> ln(5.); # Zum Vergleich
> logar(-2,1);
> logar(2,1.5);
> logar(Pi,10);
> ln(Pi);
> ln(evalf(Pi));  # Zum Vergleich

# 3) Regula falsi
# Die Nullstelle eines Graphen wird iterativ berechnet, indem die Kurve
# durch eine Sekante durch zwei Punkte angenaehert wird.
# Beim naechsten Durchlauf ersetzt die soeben gefundene Naeherung fuer die
# Nullstelle einen der beiden Punkte des vorigen Durchlaufs.
# Hinweis: Wenn keine Konvergenz erreicht wird, sollte man die Prozedur
# noch einmal mit zwei besseren Schaetzwerten starten.
> restart;
> reg:=proc(f::procedure,s1::realcons,s2::realcons)
> local x1,x2,xneu,i;
> description "Nullstellenbestimmung mit der regula falsi";
> x1:=evalf(s1);
> x2:=evalf(s2);
> i:=0;
> while abs(f(x2))>1.0E-8 do
>    if abs(f(x2)-f(x1))<1.0E-6 or i=1000
>    then error "Keine Konvergenz" end if;
> xneu:=x2-f(x2)*(x1-x2)/(f(x1)-f(x2));
> i:=i+1;
> x1:=x2;
> x2:=xneu
> end do
> end proc:
# 
> reg(sin,1.5,4.5);
> fsolve(sin,1.5..4.5); # Zum Vergleich
> reg(x->x^2-1,0,10);
> reg(x->x^2-1,10,10);
> reg(x->x^2-1,-10,10);
> reg(x->x^2-1,0,-10);
> reg(x->x^2-64,10000000,10000001);
> reg(x->x^2+64,10,20); # Es existiert keine Nullstelle
> reg(x->sin(x)-1/2,0,Pi/2);
> fsolve(sin(x)-1/2,x=0..Pi/2); # Zum Vergleich
> s:=solve(sin(x)-1/2,x); evalf(s); # Zum Vergleich exakt
> reg(x->x^3-1000,5,20);
> reg(x->x^3-1000,0,100);
> debug(reg): # Suche nach der Ursache des Versagens
> reg(x->x^3-1000,0,100);