2017-04-07 4 views
3

Ich habe eine Reihe von Funktionen f_t mit mehreren Wurzeln (zwei tatsächlich). Ich möchte die "erste" Wurzel finden und dies mit fsolve funktioniert die meiste Zeit gut. Das Problem ist, dass die beiden Wurzeln konvergieren, wenn t in die Unendlichkeit geht. (Ein einfaches Beispiel für meine Funktionen wäre f_t(x) = x^2 - 1/t). Je größer t wird, desto mehr Fehler macht fsolve. Gibt es eine vordefinierte Funktion, ähnlich wie fsolve, zu der ich sagen kann, sollte es nur in einem bestimmten Bereich (z. B. finden Sie immer die Wurzel in [0, inf)).Finden Sie eine Wurzel einer Funktion in einem bestimmten Bereich

Die Frage ist im Wesentlichen das gleiche wie https://mathematica.stackexchange.com/questions/91784/how-to-find-numerically-all-roots-of-a-function-in-a-given-range?noredirect=1&lq=1, aber die Antworten gibt es für Mathematica, ich will sie in Python.

PS: Ich weiß jetzt, wie ich meinen eigenen Algorithmus schreiben kann, aber da diese langsamer als Builtins sind, habe ich gehofft, einen eingebauten zu finden, der dasselbe macht. Vor allem habe ich diesen Beitrag gelesen Find root of a function in a given interval

+1

(https://docs.scipy.org/doc/scipy-0.18.1/reference/optimize.html# root-finding) – Michael

Antwort

1

Es ist allgemein akzeptiert, dass für glatte, gut ertragene Funktionen die Brent method die schnellste Methode ist, die garantiert eine Wurzel gibt. Wie bei den anderen beiden aufgelisteten Methoden müssen Sie ein Intervall [a, b] angeben, über das die Funktion kontinuierlich ist und das Vorzeichen ändert.

Die Scipy-Implementierung ist dokumentiert here. Ein Beispiel Anwendungsfall für die Funktion erwähnt Sie könnte wie folgt aussehen: [? Versuchte dieses]

from __future__ import division 
import scipy 

def func(x,t): 
    return(x**2 - 1/t) 

t0 = 1 
min = 0 
max = 100000 # set max to some sufficiently large value 

root = scipy.optimize.brentq(func, min, max, args = (t0)) # args just supplies any extra 
                 # argument for the function that isn't the varied parameter 
+0

Nur Kosmetik, aber Sie könnten 't' definieren, für Python2 könnte es auch besser sein,' 1./t' zu verwenden. +1, nette Lösung. – Cleb

+1

guten Ruf, danke. Ich diskutierte, ob ich mich in diesem kleinen Code über die Teilung Gedanken machen sollte oder nicht. Ich habe gerade das faule Ding mit einem 'from __future__' gemacht. – Alex

+0

Irgendeine Idee, ob man das benutzen kann, um alle Wurzeln in einem gegebenen Bereich wie in meiner Antwort unten zu finden? Ich war überrascht, das nicht in Standard-Scipy zu finden. – Cleb

2

Sie scipy.optimize.bisect verwenden können, die a und b zwei Parameter nimmt, die das Startintervall definieren. Es gibt ein paar Einschränkungen, aber:

  • Das Intervall muss endlich sein. Sie können nicht in [0, inf] suchen.
  • Die Funktion muss Zeichen an der Wurzel Flip (f(a) und f(b) müssen entgegengesetzte Vorzeichen haben), so zum Beispiel, können Sie die Wurzel f(x) = abs(x) nicht finden können (Wenn das auch ein „root“ im mathematischen Sinne ist). Es funktioniert auch nicht für f(x) = x**2 - 1 und ein Intervall von [a, b] mit < -1 und b> 1.
  • Die Methode basiert nicht auf Gradienten. Dies kann ein Vorteil sein, wenn die Funktion sehr zackig oder teuer zu bewerten ist, aber bei anderen Funktionen möglicherweise langsamer ist.

Eine Alternative ist scipy.optimize.minimize zu verwenden, um abs(f(x)) zu minimieren. Diese Funktion kann bounds mit Unendlichkeit annehmen. Die Minimierung kann jedoch in einem lokalen Nicht-Root-Minimum der Funktion enden.

+0

Der Umgang mit Unendlichkeit ist in der Tat ein Problem, auch in meiner Antwort unten basierend auf 'Chebpy'. Aber ich denke, für praktische Zwecke kann man einfach eine große Zahl einstecken. – Cleb

+1

@Cleb, das ist richtig, aber eine große Anzahl wird die Anzahl der Bisektionsschritte erhöhen. Ein Schritt pro Doppelbereich. – kazemakase

1

Classically, könnten Sie root verwenden:

import numpy as np 
from scipy.optimize import root 

def func(x, t): 
    return x ** 2 - 1./t 

t = 5000 

res = root(func, 0.5, args=(t,)).x[0] 
print res 

, dass die positiven drucken würde, in diesem Fall 0.0141421356237.

Wenn Sie den Bereich angeben möchten und bestimmen alle Wurzeln innerhalb dieses Intervalls, können Sie chebpy:

from chebpy import chebfun 

x = chebfun('x', [-100000, 100000]) 
t = 5000 
f = x ** 2 - 1./t 

rts = f.roots() 
print rts 

Dies sowohl die positive drucken würde und die negative Wurzel, in diesem Fall

[-0.01413648 0.01413648] 

Wenn Sie nur im positiven Bereich suchen möchten, können Sie

x = chebfun('x', [-100000, 100000]) 
ändern

zu

x = chebfun('x', [0, 100000]) 

Ich bin aber nicht sicher, wie die Unendlichkeit verwenden, aber Sie können nur eine sehr hohe Zahl für praktische Zwecke verwenden, glaube ich.

Verwandte Themen