Nullstellenbestimmung

[RauberRacing]

Servus.
Wollte mir mal ein Programm zur Näherungsweisen Nullstellenbestimmung einer stetigen Fkt. schreiben (nach Nullstellensatz) und habe jetzt folgenden Code.
Ich kriege jetzt aber mit jedem a, b wert und jeder Fkt nur "None" raus... wo ist mein Fehler? Ich find ihn ums verrecken net
Der geht auch net in die Rekursion, nie Oo (deswegen hab ich mal ausgabe von a, b, m reingeholt um das zu sehen

Code: Alles auswählen

def Nullstelle(a, b):
    m = (a + b)/2.0
    print "a = %f , b = %f , m = %f" %(a, b, m)
    fa = a**2 - 3
    fb = b**2 - 3
    fm = m**2 - 3
    if fa == 0:
        return a
    elif fb == 0:
        return b
    elif (fa < 0 and fm >= 0) or (fa > 0 and fm <= 0):
        return Nullstelle(a, m)
    elif (fb < 0 and fm >= 0) or (fb > 0 and fm <= 0):
        return Nullstelle(m, b)

[EyDu]

Wenn niemals die Rekusrion betreten wird und du None zurückgeliefert bekommst, sind offensichtlich deine Bedingungen verkehrt.

[RauberRacing]

hab jetzt mal ne neue version... jetzt komm ich aber nich mehr aus der rekursion raus obwohl der anker eindeutig is

Code: Alles auswählen

def Nullstelle(a, b):
    m = (a + b)/2.0
    fa = (a**2)
    fb = (b**2)
    fm = (m**2)

    if fa == 0:
        return a
    elif fb == 0:
        return b
    elif fm == 0:
        return m
    elif (fa < 0 and fm > 0) or (fa > 0 and fm < 0):
        return Nullstelle(a, m)
    elif (fb < 0 and fm > 0) or (fb > 0 and fm < 0):
        return Nullstelle(m, b)

[EyDu]

Du kannst nicht einfach eine Fließkommazahl mit "==" vergelichen, da diese typischerweise Rundungsfehlern unterliegt. Wenn du also testen willst:

Code: Alles auswählen

if a == x:

dann musst du so etwas verwenden:

Code: Alles auswählen

 if (x-epsilon) < a < (x+epsilon):

epsilon ist natürlich > 0 und geeignet zu wählen.

Edit: Es sollte in der Ungleichung natürlich "<" heißen

[BlackJack]

Ausserdem bekomme ich bei der neuen Version auch immer `None` und keine Endlosrekursion. Wobei "endlos" sind Rekursionen in CPython ja sowieso nicht - es kommt recht schnell eine Ausnahme.

[numerix]

Dazu kommt - neben anderen Dingen, die man eleganter lösen sollte -, dass es mit der im zweiten Beispiel gewählten Funktion nie funktionieren würde.
Sie hat zwar eine Nullstelle, aber keinen Vorzeichenwechsel. Darum wirst du sie mit dem gewählten Bisektionsverfahren nur zufällig finden, wenn dein Startintervall von der Form [-x;x] ist - und dann gleich beim ersten Mal treffen.
Aber auch das gilt natürlich nur, wenn du den Rat von EyDu beherzigst und mit den Fließkommazahlen "typangemessen" umgehst.

[RauberRacing]

Du kannst nicht einfach eine Fließkommazahl mit "==" vergelichen, da diese typischerweise Rundungsfehlern unterliegt. Wenn du also testen willst:

Code: Alles auswählen

if a == x:

dann musst du so etwas verwenden:

Code: Alles auswählen

 if (x-epsilon) > a > (x+epsilon):

epsilon ist natürlich > 0 und geeignet zu wählen.

das versteh ich grad irgendwie nicht. ich komm mit dem epsilon nich klar Oo was soll das sein?

@pütone: hab das eigtl mit f(x) = x²-1 gemacht habs anscheinend nur falsch kopiert

[BlackJack]

An `epsilon` musst Du eine kleine Zahl binden die den Bereich angibt, in dem Du zwei Zahlen noch als "gleich" betrachten möchtest. Fliesskommazahlen sind ungenau, darum kann man nicht immer auf exakte Gleichheit testen. Beispiel:

Code: Alles auswählen

In [437]: sum([0.1] * 10) == 1.0
Out[437]: False

In [438]: sum([0.1] * 10)
Out[438]: 0.99999999999999989

[numerix]

Es geht darum, dass Fließkommazahlen nur mit einer begrenzten Genauigkeit gespeichert werden können, und auch darum, dass bei einem Näherungsverfahren der (mathematisch) exakte Wert u.U. gar nicht getroffen werden KANN, etwa wenn es sich - wie bei deinem ersten Beispiel - um irrationale Nullstellen handelt.

Du darfst darum nicht prüfen, ob ein Funktionswert GLEICH Null ist, sondern musst abbrechen, wenn du eine von dir festgelegte (aber oberhalb der Maschinengenauigkeit liegende) Genauigkeit erreicht hast, also etwa wenn die Differenz zwischen dem Funktionswert und Null kleiner als 0,0000001 ist.

Das kannst du so machen wie EyDu es dir gezeigt hat oder z.B. durch

if abs(fa)<0.0000001: ...

wobei der von mir willkürlich gewählte Wert 0.0000001 dem Epsilon entspricht.

[RauberRacing]

so, hab jetzt verstanden was ihr meint... schonmal danke jetzt läufts... aber gibts irgendeine Möglichkeit jetzt die Funktion nur einmal eingeben zu müssen (im Idealfall als raw_input) und die variablen fa, fb und fm automatisch einen wert zugewiesen bekommen also irgendwie durch ne for schleife und string format oder irgendwas in der art?

Danke schonmal für die viele Hilfe

Code: Alles auswählen

def Nullstelle(a, b):
    m = (a + b)/2.0
    fa = (a**2)-1
    fb = (b**2)-1
    fm = (m**2)-1
    epsilon = 0.0000000001
    if -epsilon < fa < epsilon:
        return a
    elif -epsilon < fb < epsilon:
        return b
    elif -epsilon < fm < epsilon:
        return m
    elif (fa < 0 and fm > 0) or (fa > 0 and fm < 0):
        return Nullstelle(a, m)
    elif (fb < 0 and fm > 0) or (fb > 0 and fm < 0):
        return Nullstelle(m, b)

[Leonidas]

aber gibts irgendeine Möglichkeit jetzt die Funktion nur einmal eingeben zu müssen (im Idealfall als raw_input) und die variablen fa, fb und fm automatisch einen wert zugewiesen bekommen also irgendwie durch ne for schleife und string format oder irgendwas in der art?

Klar:

Code: Alles auswählen

def nullstelle(a, b):
    m = (a + b)/2.0
    fs = [(n**2)-1 for n in (a, b, m)]
    epsilon = 0.0000000001
    if -epsilon < f[0] < epsilon:
        return a
    elif -epsilon < f[1] < epsilon:
        return b
    elif -epsilon < f[2] < epsilon:
        return m
    elif (f[0] < 0 and f[2] > 0) or (f[0] > 0 and f[2] < 0):
        return nullstelle(a, m)
    elif (f[1] < 0 and f[2] > 0) or (f[1] > 0 and f[2] < 0):
        return nullstelle(m, b)

Aber ob's lesbarer ist?

[BlackJack]

@RauberRacing: Als erstes solltest Du die Nullstellen-Funktion so ändern, dass sie die Funktion als Argument bekommt.

Dann kann es, je nach Funktion, immer noch passieren, dass alle ``if``/``elif`` nicht greifen, Du solltest also noch ein ``else`` am Ende einfügen, das diesen Fall behandelt.

Letztlich würde ich die Funktion nicht rekursiv schreiben. Besser wäre eine Schleife mit einer maximalen Laufzahl, so dass man ggf. mit einer eigenen Ausnahme abbrechen kann, falls keine Lösung gefunden wird.

Wenn der Benutzer die Funktion eingeben können soll, kommt man wohl um `eval()` nicht herum.

[CM]

Hoi,

falls die zu evaluierende Funktion ein einfaches Polynom ist, geht auch:

Code: Alles auswählen

In [1]: from scipy import poly1d

In [2]: p = poly1d([3,4,5])

In [3]: print p
   2
3 x + 4 x + 5

In [4]: p.deriv()
Out[4]: poly1d([6, 4])

In [5]: p.integ()
Out[5]: poly1d([ 1.,  2.,  5.,  0.])

In [6]: p.variable
Out[6]: 'x'

In [7]: p.order
Out[7]: 2

Aber grundsätzlich sollte die Funktion auch erst einmal auf sicheren Füßen stehen - da hat BJ schon recht. (Oder Du fällst auf scipy.optimize.fsolve zurück .)

Gruß,
Christian

[numerix]

Ich würde es auch eher ohne Rekursion machen, obwohl das Bisektionsverfahren von der Idee her in der Tat ein rekursives Verfahren ist. Das kann man sich z.B. gut zu nutze machen, um auch schon mit einfachen Schultaschenrechnern der neuesten Generation (z.B. TI-30 X II S oder TI-30 Multiview oder die Casio-Modelle der ES-Reihe) solche Verfahren durch Wiederaufrufen der letzten Zeile(n) umzusetzen - ohne Programmierung oder Taschenrechner aus anderen "Gewichtsklassen".

Konkret könnte man es so machen:

Code: Alles auswählen

def nullstelle(f,a,b):
    while b-a>0.00000001:
        m = 0.5*(a+b)
        if f(a)*f(m)<0: b = m
        else: a = m
    return m

Darin ist die sinnvolle Empfehlung von BlackJack, die Funktion f mit an die Nullstellenfunktion zu übergeben schon eingebaut.

Den Tipp zur Verwendung von "eval" hast du ja schon bekommen.
Damit kannst du den Term einlesen und daraus dann die Funktion f basteln - entweder kurz mit einem Lambda-Ausdruck oder - wenn dir das nicht zusagt - über eine normale Funktionsdefinition.

P.S. Das ist natürlich jetzt alles ohne Netz und doppelten Boden. Es muss natürlich sichergestellt sein, dass f in [a,b] überall definiert ist, dass es eine Nullstelle mit Vorzeichenwechsel gibt, dass a<b ist usw.

[RauberRacing]

was ist das denn mit den In[1] Out[2] usw?
sry bin noob auf dem gebiet

[Leonidas]

was ist das denn mit den In[1] Out[2] usw?

Die Ausgaben von IPython, der das statt ``>>>`` anzeigt.

[numerix]

So könnte man es auch machen:

Code: Alles auswählen

def nullstelle(f,a,b):
    while b-a>0.0000001:
a,b=(not(f(a)*f(((a+b)*.5))<0))*((a+b)*.5)+(f(a)*f(((a+b)*.5))<0)*a,(f(a)*f(((a+b)*.5))<0)*((a+b)*.5)+(not(f(a)*f(((a+b)*.5))<0))*b
    return (a+b)*.5

P.S. Irgendwie klappt es mit der richtigen Darstellung der Einrückung nicht ..