________________________________________________________________________________
Inhalt....: Numerisches Lösen von nicht-polynomialen Gleichungen
Kategorie.: Handwerkskasten
Mathematik: Numerik
MuPAD.....: 3.1.1
Datum.....: 2005-04-04
Autoren...: Kai Gehrs <acrowley@mupad.de>
Autoren...: Andreas Sorgatz <sorgatz@sciface.com>
Funktionen: solve, numeric::solve, numeric::realroots, plotfunc2d, YRange
Funktionen: lhs, rhs, AllRealRoots
________________________________________________________________________________
Numerisches Lösen von
nicht-polynomialen Gleichungen
Wir betrachten folgende nicht ganzrationale Gleichung und versuchen diese
mit Hilfe des üblichen in MuPAD Befehls solve zu lösen:
G:= x^3 + 1 = exp(x)
solve(G, x)
Dieser Ansatz führt aber nicht zum Ziel, denn eine algebraische Darstellung
der Lösung ist in diesem Fall ganz allgemein nicht möglich. Da hilft nur eine
näherungsweise (d.h. numerische) Berechnung der Nullstellen.
Versuchen wir es dabei zunächst mit einem direkten Ansatz und gehen später
zu einem effektiveren Verfahren über. Wir verwenden numeric::solve, um
Näherungen der Nullstellen zu berechnen:
L1:= numeric::solve(G, x)
Diese Funktion liefert uns höchstens eine Nullstelle der Gleichung, denn die
existierenden numerischen Verfahren sind Intervall-orientiert (Newton- oder
Bisektionsverfahren) und finden pro Intervall (hier die reellen [bzw. komplexen]
Zahlen) höchstens eine Nullstelle. Insbesondere bei stark oszillierenden
Funktionen ist es sogar so, dass man aus Effizienzgründen gar nicht alle
Nullstellen berechnet haben möchte (es gibt Funktionen, die in einem winzigen
Intervall viele tausend Nullstellen besitzen).
Beim direkten Ansatz ist es daher sinnvoll sich zunächst einen Eindruck von
den Graphen der beteiligten Funktionen zu verschaffen: Term der linken Seite
der Gleichung (blau), rechte Seite (rot) und die Differenz (grün).
plotfunc2d( lhs(G), rhs(G), lhs(G)-rhs(G),
x = -2..2, YRange = -2..6,
GridVisible )
Neben der bereits gefundenen Nullstelle der grünen Funktion weist die Zeichnung
auf die Existenz von weiteren Nullstellen in den Intervallen [-1..-0.5] und [1..2] hin.
Wir verwenden also erneut numeric::solve, schränken diesmal aber die Such-
intervalle entsprechend ein:
L2:= numeric::solve( G, x=-1..-0.5 );
L3:= numeric::solve( G, x=1..2 )
Sind das nun alle Lösungen? Oder liegen ggf. ausserhalb des gezeichneten
Bereiches noch Nullstellen? Wir betrachten daher noch die Randintervalle:
L4:= numeric::solve( G, x=-1000..-2 );
L5:= numeric::solve( G, x=1.55..1000 )
und
L6:= numeric::solve( G, x=4.57..1000 )
und vereinigen die einzelnen Lösungsmengen:
L = _union(L1, L2, L3, L4, L5, L6)
Die Lösungsmenge L scheint nun vollständig zu sein, was anhand der folgen-
den Zeichnung plausibel ist und mit Kenntnis über die Funktionen und deren
streng monotonem Verlauf ausserhalb des gezeichneten Intervalls begründet
werden kann.
plotfunc2d( lhs(G) - rhs(G),
x = -3..5, YRange = -20..20,
GridVisible )
Aus didaktischer Sicht ist an dieser Herangehensweise schön, dass man sich
aktiv mit den Funktionen und deren Graphen auseinandersetzt.
Aber es ist kein Verfahren, dass man als "Black-Box" einsetzen kann, um schnell
zu Ergebnissen zu kommen. Zudem muss man sich möglicher numerischer
Probleme bewusst sein und sollte die Teilintervalle nicht unsinnig groß wählen,
da die Nullstellen sonst ggf. nicht gefunden werden, sondern "durchrutschen":
numeric::solve( G, x=1..10 ) <>
numeric::solve( G, x=1..10^8 )
Ab MuPAD 3.1.1 steht für die Funktion numeric::solve eine neue Option
AllRealRoots ("Alle reellen Nullstellen") zur Verfügung. Diese ermittelt über
eine geschickte Heuristik Suchintervalle und berechnet möglichst viele Null-
stellen. Die Option arbeitet in vielen Praxisfällen einwandfrei und erspart eine
eigene Angabe von Suchintervallen. Aufgrund der im folgenden geschilderten
Problematiken kann das Verfahren aber nicht als ganz sicher angesehen werden.
L = numeric::solve( G, x, AllRealRoots )
Warning: problem in isolating search intervals. Some roots may be lost [numeric::allRealRoots]
Es ist oft günstiger sich für einen zu untersuchenden Bereich die Teilintervalle
zur Nullstellensuche gezielt von MuPAD berechnen zu lassen. Dazu dient die
Funktion numeric::realroots. Sie garantiert für den untersuchten Bereich,
dass ausserhalb der berechneten Intervalle keine Nullstelle existiert. Ihr über-
geben wir den zu untersuchenden Bereich x = -10^5..10^5 sowie die
Option Merge, die bewirkt, dass benachbarte Intervalle zusammengefasst
werden, sofern es sinnvoll erscheint:
LI:= numeric::realroots(G, x=-10^5..10^5, Merge)
Wir erzeugen aus den Teillisten [a,b] von LI nun Bereiche a..b
LB:= map(LI, l -> l[1]..l[2])
und, da wir aufgrund des Funktionsgraphen wissen, dass in jedem Intervall genau
eine Nullstelle liegt, übergeben diese an numeric::solve, um Näherungen der
Nullstellen zu berechnen:
numeric::solve(G, x = b) $ b in LB;
L = _union(%)
Genäherte Zahlwerte zu den einzelnen Intervallen erhalten wir auch indem wir die
Größe der Intervalle an die aktuelle Rechengenaugikeit 10^(-DIGITS) angleichen
und dann einfach die linke Grenze des Intervalls nehmen:
LI:= numeric::realroots(G, x = -10^5..10^5,
10^(-DIGITS), Merge);
L = {l[1] $ l in LI}
ABER VORSICHT:
Es können in LI auch Intervalle auftauchen, in denen keine oder mehrere Null-
stellen liegen. Die Funktion numeric::realroots liefert uns mit den Intervallen
also im wesentlichen ein Ausschlußkriterium für Nullstellen.
Betrachten wir daher die Funktion sin(1/x):
H:= sin(1/x) = 0;
plot(
plot::Scene2d( plot::Function2d(lhs(H), x = -15..15),
GridVisible ),
plot::Scene2d( plot::Function2d(lhs(H), x = -1..1),
GridVisible ),
Width = 180
)
Wir erkennen, dass der Funktionsgraph mit zu wenigen Stützpunkten gezeichnet
wurde (einer zu geringen Abtastrate). Mittels Mesh = 800 erhält man eine realere
Darstellung des Graphen. Das Attribut Mesh reguliert die Anzahl der Stützstellen,
die MuPAD intern verwendet, um einen Funktionsgraphen zu zeichnen (ganz
ähnlich zu der Vorgehensweise, wenn man per Hand einen Funktionsgraphen
zeichnet: Auch dann berechnet man zu einigen x-Koordinaten die zugehörigen
y-Koordinaten, trägt die Punkte in ein Koordinatensystem ein und zeichnet den
Graphen durch die vorgegebenen Punkte). Zusätzlich schränken wir den
Definitionsbereich auf das eigentlich interessante Intervall [-1;1] um den Null-
punkt ein:
plot(plot::Function2d(lhs(H), x = -1..1, Mesh = 800,
GridVisible), Width = 180)
Nahe des Nullpunktes hat die Funktion oszilliert die Funktion sin(1/x) wild und
besitzt unendlich viele Nullstellen.
Genauer: Jedes 1/(k*PI) mit k aus der Menge der ganzen Zahlen ohne die Null
ist Lösung.
solve(H, x)
Die Funktion numeric::solve findet ohne eine Bereichsangabe keine Nullstelle:
numeric::solve(H, x)
Wir berechnen uns Nullstellenintervalle und Näherungen als linke Grenzen der
Intervalle
IL:= numeric::realroots(H, x = -1..1, Merge);
Intervalle = nops(IL);
L = {l[1] $ l in IL}
und erhöhen die Genauigkeit auf eine Intervallbreite von 10^-4:
numeric::solve(H, x);
IL:= numeric::realroots(H, x = -1..1, 10^-4, Merge):
Intervalle = nops(IL);
L = {l[1] $ l in IL}
Fazit: Eine Nullstellensuche ist alles andere als trivial. Kann die Lösung nicht
geschlossen dargestellt oder nicht von den algebraischen Algorithmen ge-
funden werden, so bleibt nur eine numerische Intervalluntersuchung.
________________________________________________________________________________
Übungen:
1. Berechnen Sie alle reellen Lösungen der folgenden Gleichungen. Stellen Sie zuvor die ent-
__sprechenden linken Seiten der Gleichungen grafisch dar, um ein Gefühl für die Lokalität der
__Nullstellen zu erhalten. Findet solve eine symbolische (exakte) Lösung?
Anmerkungen:
1. Weitere Anregungen finden Sie in der Buchreihe Mathematik 1 x anders. In dieser Reihe
wird eine Vielzahl unterschiedlichster mathematischer Probleme mit MuPAD gelöst. Die
Bücher können unter www.mupad.de/schule+studium/literatur kostenfrei kopiert werden.
________________________________________________________________________________