Nullstellenbestimmung

Wenn du dir nicht sicher bist, in welchem der anderen Foren du die Frage stellen sollst, dann bist du hier im Forum für allgemeine Fragen sicher richtig.
RauberRacing
User
Beiträge: 25
Registriert: Montag 19. November 2007, 10:03
Wohnort: Schwarzenbach
Kontaktdaten:

Dienstag 19. Februar 2008, 19:42

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
User
Beiträge: 4871
Registriert: Donnerstag 20. Juli 2006, 23:06
Wohnort: Berlin

Dienstag 19. Februar 2008, 20:12

Wenn niemals die Rekusrion betreten wird und du None zurückgeliefert bekommst, sind offensichtlich deine Bedingungen verkehrt.
RauberRacing
User
Beiträge: 25
Registriert: Montag 19. November 2007, 10:03
Wohnort: Schwarzenbach
Kontaktdaten:

Dienstag 19. Februar 2008, 20:16

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
User
Beiträge: 4871
Registriert: Donnerstag 20. Juli 2006, 23:06
Wohnort: Berlin

Dienstag 19. Februar 2008, 20:22

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 :oops:
Zuletzt geändert von EyDu am Dienstag 19. Februar 2008, 21:10, insgesamt 1-mal geändert.
BlackJack

Dienstag 19. Februar 2008, 20:32

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.
Benutzeravatar
numerix
User
Beiträge: 2696
Registriert: Montag 11. Juni 2007, 15:09

Dienstag 19. Februar 2008, 20:34

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
User
Beiträge: 25
Registriert: Montag 19. November 2007, 10:03
Wohnort: Schwarzenbach
Kontaktdaten:

Dienstag 19. Februar 2008, 20:39

EyDu hat geschrieben: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

Dienstag 19. Februar 2008, 20:47

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
Benutzeravatar
numerix
User
Beiträge: 2696
Registriert: Montag 11. Juni 2007, 15:09

Dienstag 19. Februar 2008, 20:54

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
User
Beiträge: 25
Registriert: Montag 19. November 2007, 10:03
Wohnort: Schwarzenbach
Kontaktdaten:

Mittwoch 20. Februar 2008, 12:06

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
Administrator
Beiträge: 16024
Registriert: Freitag 20. Juni 2003, 16:30
Kontaktdaten:

Mittwoch 20. Februar 2008, 12:24

RauberRacing hat geschrieben: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?
My god, it's full of CARs! | Leonidasvoice vs Modvoice
BlackJack

Mittwoch 20. Februar 2008, 12:32

@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
User
Beiträge: 2464
Registriert: Sonntag 29. August 2004, 19:47
Kontaktdaten:

Mittwoch 20. Februar 2008, 12:45

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
Benutzeravatar
numerix
User
Beiträge: 2696
Registriert: Montag 11. Juni 2007, 15:09

Mittwoch 20. Februar 2008, 17:30

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
User
Beiträge: 25
Registriert: Montag 19. November 2007, 10:03
Wohnort: Schwarzenbach
Kontaktdaten:

Mittwoch 20. Februar 2008, 18:35

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