Algorithmus um Nullstellen von Polynomen zu berechnen

[umask007]

Hallo,
ich schreibe gerade einen neuen Algorithmus, um die Nullstellen von Polynomen beliebigen Grades zu berechnen.
Der Algorithmus basiert auf dem Newtonverfahren (https://de.wikipedia.org/wiki/Newtonverfahren).
Mittels der Formel x=x-f(x)/f'(x) nähert sich das Newtonverfahren schrittweise einer Nullstelle.
Das Problem bei dem Verfahren ist, dass es scheitert, wenn f'(x) null ist, da man bekannterweise nicht durch null teilen kann. Meine Idee zur Verbesserung des Verfahrens ist simpel: Man teilt
die Funktion in Intervalle auf, in denen f'(x) != null ist, führt das Verfahren separat für jedes Intervall aus
und fügt die Ergebnisse zusammen. Man braucht um die Nullstellen zu berechnen also zuerst die Nullstellen
der Ableitung, weshalb zuerst die Nullstellen aller Ableitungen in aufsteigender Reihenfolge berechnet werden.
Hier der Beispielcode:

export function superacid(polynom:Decimal[],accuracy:Decimal=new Decimal(0.01), newtonIterations:number=32):Decimal[]
{
 /*
  Simplify polynomial by dividing through the first coefficient.
  Smaller coefficients mean less chance of overflowing.
 */
 polynom=simplyfyPolynom(polynom);
 var degree = polynom.length-1;
 if (degree<=0) return null;
 var polynoms=calculateDerivates(polynom); //calculate all derivatives
 /*
   Add polynomial to the list of polynomials, whose roots need to be calculated
 */
 polynoms.unshift(polynom);
 var c= degree;
 var firstDegreeDerivation=polynoms[polynoms.length-2];
 var derivateZeros:Decimal[]=[firstDegreeDerivation[1].neg().div(firstDegreeDerivation[0])];
 var roots = derivateZeros;
 /*
  Calculate roots of polynomial and all its derivates, starting with lowest degree.
  Roots are needed for calculating the roots of the polynom with the next degree. 
 */
 while (--c)
 {
    var fx= polynoms[c-1];
    var fxx= polynoms[c];
    roots=[];
    /*
      Intervals are located from -infinity to the first stationary point, in between neigbouring
      stationary points, and from the last stationary point to +infinity. Because of that, the zeros
      of the derivate are needed, calculated in the last loop iteration.
      In each interval there is at most one zero, because the sign of the derivative is constant.
    */ 
    var intervals= getIntervalsInPolynom(...derivateZeros);
    for( var interval of intervals)
    {  //Do Newton iterations seperately for each interval
       var root=newtonIterate(fx,fxx,interval.left,interval.right,interval.start,newtonIterations,accuracy);
       if (root==null) continue;   
       //skip next interval, because root will be the same as in the current interval 
       if (interval.right.minus(root).abs().lessThan(accuracy)) intervals.next();
       roots.push(root); //Add root to the list of found roots   
    }
    derivateZeros=roots; //Roots are needed for calculating the next intervals   
 }
 return roots; // return found roots
}

Mich würde mal Eure Meinung zum neuen Algorithmus interessieren. Am interessantesten ist für mich, ob der Algorithmus auch alle Nullstellen findet. Da ich kein Mathematiker bin, habe ich auch keine Ahnung, wie man einen formellen Beweis für einen solchen Algorithmus schreiben könnte. Habe den Algorithmus stattdessen intensiv mit Unit Tests geprüft. Der Algorithmus kann sogar das Wilkinson Polynom lösen (https://en.wikipedia.org/wiki/Wilkinson's_polynomial) von daher denke ich, dass der Algorithmus auf jeden Fall auf die ein oder andere Weise nützlich sein wird.
Hier ist der komplette Code:
https://github.com/P3N9U1N/superacid
Hier ist eine GUI zum Testen:
https://p3n9u1n.github.io/superacid/gui/index.html

Viele Grüße, Thomas

 

[Rumbah]

Es ist ein numerisches Näherungsverfahren. Daher kannst du allein über die Rechengenauigkeit in Probleme laufen. Auch musst du die Suchstartpunkte geschickt wählen, damit du auch alle Nullstellen erwischst. Und im Reelen weißt du auch nicht, wie viele Nullstellen überhaupt da sind.
Du löst das mit den Startpunkten und der Anzahl ja durch die erste Ableitung, aber das verschiebt das Problem ja nur auf ein Polynom mit einem Grad niedriger. Da hast du dann natürlich wieder die gleichen Probleme.

z.B. findet deine Gui für das Polynom x^2+0.0001 die zwei nicht vorhandenen Nullstellen +-0.003 .
Wahrscheinlich könnte man auch was konstruieren, wo du mit dem falschen Anfangswert in eine Periode läufst.
Aber für "normale" Polynome bist du mit deiner Idee wahrscheinlich gut dran.

 

[umask007]

Ja es ist ein numerisches Näherungsverfahren. Laut dem Abel-Ruffini Theorem
gibt es allerdings keine Formel für Polynome mit 5. Grad oder höher, Näherungsverfahren
sind hierfür also alternativlos. Jedes Näherungsverfahren hat einen Schwellenwert,
innerhalb dessen Zahlen als gleich gelten. Der Punkt ist, dass man bei der Nullstellenberechnung
diesen Schwellenwert braucht, da die Null selten genau getroffen wird. Der von dir genannte Problemfall
lässt sich lösen, indem der Schwellenwert verringert und die Präzision erhöht wird, decimal.js hat hier noch
sehr viel Luft nach oben, auf Kosten der Performance.

Wahrscheinlich könnte man auch was konstruieren, wo du mit dem falschen Anfangswert in eine Periode läufst.

Das würde ich ausschließen lasse mich aber gerne eines besseren belehren, dient es doch zur Verbesserung des Algorithmus.

 

[BeBur]

Jedes Näherungsverfahren hat einen Schwellenwert,
innerhalb dessen Zahlen als gleich gelten.

Man muss hier noch unterscheiden zwischen dem Algorithmus an sich und der Implementierung. Eine Implementierung ausschließlich mit floats hat die zugehörigen bekannten Probleme, ggf. auch innerhalb eines im Algo angegebenen Schwellwertes.

Rumbah hat davon abgesehen etwas anderes angesprochen, nämlich das mehr Nullstellen gefunden werden als existieren. Das ist qualitativ noch etwas anderes als leicht daneben zu liegen.

Ich bin mir nicht ganz sicher, ob dir das klar ist, aber natürlich existiert eine lange Reihe von Varianten und Alternativen zum Newton-Verfahren. Ich will deine Arbeit nicht abwerten, das ist schon ganz cool, sich selber was auszudenken und dann zu bauen. Allerdings lautete dein erster Satz "ich schreibe gerade einen neuen Algorithmus, um die Nullstellen von Polynomen beliebigen Grades zu berechnen." und neu ist so etwas simples natürlich nicht.

 

[umask007]

"ich schreibe gerade einen neuen Algorithmus, um die Nullstellen von Polynomen beliebigen Grades zu berechnen." und neu ist so etwas simples natürlich nicht.

Klar ist der neu:
https://en.wikipedia.org/wiki/Root-finding_algorithms
Wie heißt der Algorithmus denn, den ich geschrieben habe?

Ich bin mir nicht ganz sicher, ob dir das klar ist, aber natürlich existiert eine lange Reihe von Varianten und Alternativen zum Newton-Verfahren.

Das weiß ich sehr gut, denn ich habe gründlich recherchiert. Ganz so viele sind es nicht, und wenn man noch die rausnimmt, die überhaupt nicht alle Nullstellen finden, bleiben nicht viele übrig:

Finding the real roots of a polynomial with real coefficients is a problem that has received much attention since the beginning of 19th century, and is still an active domain of research. Most root-finding algorithms can find some real roots, but cannot certify having found all the roots. Methods for finding all complex roots, such as Aberth method can provide the real roots. However, because of the numerical instability of polynomials (see Wilkinson's polynomial), they may need arbitrary-precision arithmetic for deciding which roots are real. Moreover, they compute all complex roots when only few are real.

It follows that the standard way of computing real roots is to compute first disjoint intervals, called isolating intervals, such that each one contains exactly one real root, and together they contain all the roots. This computation is called real-root isolation. Having isolating interval, one may use fast numerical methods, such as Newton's method for improving the precision of the result.

The oldest complete algorithm for real-root isolation results from Sturm's theorem. However, it appears to be much less efficient than the methods based on Descartes' rule of signs and Vincent's theorem. These methods divide into two main classes, one using continued fractions and the other using bisection. Both method have been dramatically improved since the beginning of 21st century. With these improvements they reach a computational complexity that is similar to that of the best algorithms for computing all the roots (even when all roots are real).

These algorithms have been implemented and are available in Mathematica (continued fraction method) and Maple (bisection method). Both implementations can routinely find the real roots of polynomials of degree higher than 1,000.

von https://en.wikipedia.org/wiki/Root-finding_algorithms

etwas simples natürlich nicht.

Wenn es so einfach ist nenne mir einen konkreten Verbesserungsvorschlag oder einen Algorithmus der besser funktioniert. Mein Algorithmus kann das Wilkinson Polynom lösen, das ist ganz bestimmt nicht einfach und das schaffen auch nicht alle Algorithmen.

Rumbah hat davon abgesehen etwas anderes angesprochen, nämlich das mehr Nullstellen gefunden werden als existieren. Das ist qualitativ noch etwas anderes als leicht daneben zu liegen.

Hab ich doch schon erklärt. Innerhalb der Eingabeparameter ist das Ergebnis korrekt. "x^2+0.001" funktioniert ohne Probleme. Mit höheren Präzisionseinstellungen geht auch "x^2+0.0001". Oder gibt es etwa einen Algorithmus, der auch mit Integer Werten Nullstellen präzise berechnen kann?

 

[BeBur]

Kannst du dir natürlich so zurechtlegen, aber in der Wissenschaft (das ist es, was du behauptest zu tun) ist man in der Bring-Schuld, du müsstest z.B. das Konvergenzverhalten deines Algos angeben (vermutlich wie beim herkömmlichen Newtonverfahren) und dir auch den Stand der Forschung in dem Bereich anschauen und nach bereits existierenden Newton-Verfahren recherchieren.
Bestimmte Probleme scheinen dir auch nicht bekannt zu sein, evtl. ist dir nur der deutsche Wikipedia Artikel ausschnittsweise bekannt? Link.
Zwei Paper zu dem Thema mit denen du anfangen könntest wären dieses hier oder dieses hier.

Das was du da vorschlägst ist eine naheliegende Wochenaufgabe für einen angehenden Informatiker oder Mathematiker im Bachelorstudium, es kann zwar sein, dass das einen Namen oder informallen Namen hat, ich würde daher aber auch nicht darauf wetten. Wenn man daher unbedingt wissen wollte, ob das einen Namen hat würde ich daher evtl. eher Übungsaufgaben in einschlägigen Veranstaltungen sichten als wissenschaftliche Veröffentlichungen.

Du hast ja sogar auch selber zitiert:

It follows that the standard way of computing real roots is to compute first disjoint intervals, called isolating intervals, such that each one contains exactly one real root, and together they contain all the roots. This computation is called real-root isolation. Having isolating interval, one may use fast numerical methods, such as Newton's method for improving the precision of the result.

Das ist ca. das was du machst und zwar mit der offensichtlichsten Methode Intervalle zu bilden, die überhaupt möglich ist. Mit der Einschränkung, dass du nicht weißt ob deine Intervalle eine Nullstelle enthalten. Scheint mir also kategeorisch schlechter zu sein als das was das von dir zitierte aussagt.

Du könntest dir auch ganz simpel mal anschauen, welche Algorithmen einschlägige Software verwenden, z.B. Matlab:

The roots function considers p to be a vector with n+1 elements representing the nth degree characteristic polynomial of an n-by-n matrix, A. The roots of the polynomial are calculated by computing the eigenvalues of the companion matrix, A.

Das scheint ein direktes Verfahren zu sein, kein iteratives. Lineare Algebra ist aber schon eine weile her bei mir. Oder Wolfram Alpha (da hab ich spontan nicht gefunden, was sie verwenden). Generell würde ich denken, dass das Konvergenzverhalten für diese simple Klasse von Funktionen das ist, was das eigentliche Interessante ist.

 

[umask007]

Das ist ca. das was du machst und zwar mit der offensichtlichsten Methode Intervalle zu bilden, die überhaupt möglich ist. Mit der Einschränkung, dass du nicht weißt ob deine Intervalle eine Nullstelle enthalten. Scheint mir also kategeorisch schlechter zu sein als das was das von dir zitierte aussagt.

Wenn es so offensichtlich ist, wieso erwähnt der Wikipedia Artikel nichts? Wenn es so offensichtlich ist, wieso ist Newton nicht auf die Idee gekommen?
https://en.wikipedia.org/wiki/Real-root_isolation
Du schreibst ja selbst dass ich etwas anderes mache. Meine Intervalle enthalten entweder 1 oder 0 Nullstellen. Ich kann dir auch erklären wieso. Zwischen den Nullstellen der Ableitung ist das Vorzeichen der Ableitung konstant. Wenn die Funktion im Intervall also die X-Achse kreuzt, kann es keine zweite Nullstelle geben, da hierfür sich das Vorzeichen der Ableitung umdrehen müsste. Das besondere an den Intervallen ist aber dass die Ableitung im Intervall niemals 0 ist, außerdem ist das Vorzeichen der Ableitung konstant. Das heißt das Newtonverfahren kann aufgrund Formel x=x-f(x)/f'(x) schonmal nicht wegen einer Nullstelle scheitern, außerdem springt es nicht hin und her, diese Fehlerquellen sind schonmal ausgeschlossen.

Bestimmte Probleme scheinen dir auch nicht bekannt zu sein, evtl. ist dir nur der deutsche Wikipedia Artikel ausschnittsweise bekannt? Link.

Natürlich ist mir das bekannt, deshalb habe ich den Algorithmus modifiziert.

Da ich kein Mathematiker bin, habe ich auch keine Ahnung, wie man einen formellen Beweis für einen solchen Algorithmus schreiben könnte

Kannst du dir natürlich so zurechtlegen, aber in der Wissenschaft (das ist es, was du behauptest zu tun) ist man in

Hab ich nie behauptet, ganz im Gegenteil, siehe Zitat oben. Deshalb frage ich um Hilfe, die du aber nicht bieten kannst, und das obwohl du meine Arbeit als "offensichtlich" "Wochenendarbeit für einen angehenden Informatiker" "einfach" bezeichnest.

 

[BeBur]

wieso erwähnt der Wikipedia Artikel nichts?

Weil es eine Enzyklopädie ist und kein Lehrbuch. Ein solches kannst zu zu rate ziehen. Auf Stack Exchange findet man sicherlich auch Fragen die deinen Algo behandeln.

wieso ist Newton nicht auf die Idee gekommen

Woher willst du das wissen? Seine Schriften sind verfügbar, die kannst du sichten.

Deshalb frage ich um Hilfe, die du aber nicht bieten kannst, und das obwohl du meine Arbeit als "offensichtlich" "Wochenendarbeit für einen angehenden Informatiker" "einfach" bezeichnest.

Du bist wie gesagt in der Bringschuld und ehrlich gesagt raubt es mir nicht den Schlaf, wenn du herum läufst und denkst, du hättest was ganz neues erfunden.

Du schreibst ja selbst dass ich etwas anderes mache. Meine Intervalle enthalten entweder 1 oder 0 Nullstellen.

Ja, du machst deutlich weniger. Die Erkenntnis die du da offenbar hattest ist Thema des ersten Semesters im Mathe (und Info-) Studium und ist im wesentlichen der Zwischenwertsatz.

Haben insbesondere f(a) und f(b) verschiedene Vorzeichen, so garantiert der Zwischenwertsatz die Existenz von mindestens einer Nullstelle von f im offenen Intervall (a,b).

 

[umask007]

Ja, du machst deutlich weniger. Die Erkenntnis die du da offenbar hattest ist Thema des ersten Semesters im Mathe (und Info-) Studium und ist im wesentlichen der Zwischenwertsatz.

Blödsinn, beim Zwischenwertsatz geht es um mindestens eine Nullstelle.

Du bist wie gesagt in der Bringschuld und ehrlich gesagt raubt es mir nicht den Schlaf, wenn du herum läufst und denkst, du hättest was ganz neues erfunden.

Weißt du, was mich vor Lachen um den Schlaf bringt? Wenn einer wie du so tut, als wäre alles so einfach aber sich dann als komplett ahnungslos herausstellt. Du hast nicht verstanden, dass der Algo keine root isolation betreibt, sondern die Intervalle so wählt, dass die Newton Iterationen ungestört laufen können. Auch hast du keine Ahnung von den anderen Algos, die Nullstellen finden. Nichtmal den Zwischenwertsatz kannst du korrekt einordnen.

Ich habe alle deine Fragen so gut es geht beantwortet. Wenn du mir unterstellst, zu plagieren bist du in der Bringschuld. Ich habe so viel recherchiert, wie man es für ein Hobbyprojekt erwarten kann.

Woher willst du das wissen? Seine Schriften sind verfügbar, die kannst du sichten.

Ganz einfach, der Newton Algorithmus funktioniert gar nicht immer, du hast selbst nen Link dazu gepostet. Ein solch wichtiger Zusatz wäre erwähnenswert.

Weil es eine Enzyklopädie ist und kein Lehrbuch. Ein solches kannst zu zu rate ziehen. Auf Stack Exchange findet man sicherlich auch Fragen die deinen Algo behandeln.

Also gleich auf Stack Exchange und nicht CB? Ich hoffe mal, dass die Leute dort hilfreicher sind als du.

 

[Web-Schecki]

Ich hoffe mal, dass die Leute dort hilfreicher sind als du.

Was genau erwartest du denn? Du wolltest Feedback.
Nun, zu allererst ist ein Algorithmus schwer zu analysieren, wenn er in Javascript notiert ist. Damit will sich niemand freiwillig herumschlagen.
Falls du denkst, dass dein Algorithmus wirklich sinnvoll ist, dann müsstest du das halt beweisen. Das heißt insbesondere, dass du ihn mit den derzeit besten bekannten Lösungsverfahren vergleichen müsstest. Das Newtonverfahren hat eher didaktischen und weniger praktisch-akademischen Wert. Die Messlatte ist also nicht besonders hoch, wenn du dich damit vergleichst. Diese Kritik musst du dir schon gefallen lassen.

Zumal immer noch nicht klar ist, was an deinem Verfahren besser als an der Vorgehensweise, die du selbst zitiert hast, nach der in jedem Intervall genau eine Nullstelle ist. Bei dir ist entweder keine oder genau eine Nullstelle im Intervall vorhanden. Das klingt erstmal schlechter. Die Idee, Intervalle zu verwenden, ist offenbar auch nicht brandneu, sondern - wie du selbst zitierst hast - eine Standardmethode für reellwertige Nullstellen.

Ich will deine Leistung nicht schmälern, aber du solltest dich einfach nicht wundern, wenn das niemanden vom Hocker haut. Wenn du da komplett selbst drauf gekommen bist, dann ist das auf jeden Fall trotzdem eine starke Leistung, nur womöglich nichts, was andere zu sehr beeindruckt.

 

[umask007]

Du hast natürlich Recht. Mein Ziel war es auch nicht, eine Lösung zu finden die besser ist als alles bisher dagewesene.
Alles hat angefangen, als eine quartische Gleichung lösen wollte. Also habe ich die Lösungsformel nachgeschlagen
und war einmal perplex. Die Lösungsformel ist kaum benutzbar. Ich wollte verstehen, wie eine solche Formel zustandekommt und
habe mich schlau gemacht, auch über kubische und quadratische Gleichungen.
Interessanterweise hatten Mathematiker früher ein geometrisches Verständnis von Polynomen und haben durch geometrische Umformungen versucht,
Lösungen für Polynome zu finden. Was hierbei allerdings immer noch benötigt wird ist ein Wurzel Algorhitmus, also entweder Heron oder Newtonverfahren.
Sogesehen handelt es sich bei den Lösungsformeln auch nur um Näherungsverfahren. Ein Wurzel Algorhitmus ist ein Algorhitmus um die Nullstellen des Polynoms
x^2-y=0 zu finden. Was ist also, wenn die Lösungsformeln gar nicht benötigt werden? Was ist, wenn man stattdessen einfach den Alghorhitmus so anpasst,
dass Polynome beliebigen Grades berechnet werden können?
Das ist noch gar nicht die verrückteste Idee, die ich hatte. Meine verrückteste Idee dreht sich um irrationale Zahlen. Wer sich mit der Entstehungsgeschichte irrationaler Zahlen
befasst, der weiß, dass diese eingeführt wurden, um kubische Gleichungen lösen zu können. Für negative Wurzeln gibt es keine Berechnungsmöglichkeit, braucht man allerdings
auch nicht, wenn man stattdessen ein Näherungsverfahren einsetzt. Jeder echte Mathematiker wird jetzt lachen, aber ich finde die Vorstellung faszinierend,
das irrationale Zahlen, eine der bedeuesten Entdeckungen der Mathematik, nur Schall und Rauch sind. Also wenn dich diese Idee nicht vom Hocker haut dann weiß ich auch nicht
was dich beeindrucken könnte.

Was die Kritik angeht, ich bin überzeugt. Ich kann nicht ausschließen, dass es den Algo nicht schon gibt. Werde die entsprechende Behauptung aus dem Repo entfernen,
außerdem füge ich einen Hinweis über die Beschränkungen des Algo hinzu. Es soll ja niemand den Algo in Flugzeugfirmware einbauen

 

[R3ddy]

So wie ich es verstanden habe, ermittelst du die Nullstellen von einem Polynom, indem du dir zunächst die lokalen Extremwerte suchst (was die Nullstellen der Ableitung sind).

Die gefundenen Extremwerte (und zusätzlich + und - unendlich) bilden dann die Grenzen von Intervallen, für welche du das Newtonverfahren anwendest. Als Startpunkt verwendest du für jedes Intervall den Mittelwert der jeweiligen Grenzen. (Die Sonderfälle mit den Intervallgrenzen von + und - unendlich ignoriere ich hier und im Folgenden.)

Du behauptest auf diese Weise alle Nullstellen finden zu können. Ich sehe hier aber keinen Mehrwert zum klassischen Newtonverfahren, da du bereits nach einem Iterationsschritt des Newtonverfahrens auf ein Extremwert landen kannst oder sich ein Zyklus ergeben kann. letzteren Fall möchte ich anhand eines Beispiels konkret zeigen:

Gegeben sei das folgende Polynom:

-1*x^8 + 8*x^6 - (73/4)*x^4 + 18*x^2 - (3/4)

Dieses hat lokale Extremwerte an den Stellen x=-2, x=0 und x=2.

Nehmen wir an, das Verfahren hat die Nullstellen exakt ermittelt. Als Startpunkte für das Newtonverfahren würden also x=-1 und x = 1 verwendet werden.

Nach einem Iterationsschritt würde auf x=-1 der Punkt x=1 folgenden. Ein Schritt weiter ergäbe wiederum x=-1 (also ein Zyklus).

Das Polynom hat die folgenden Nullstellen (auf 10 Stellen gerundet)
x = -2.261421293
x = -0.2086944943
x = +0.2086944943
x = +2.261421293
Das von dir implementierte Verfahren würde bei exakter / ungerundeter Rechnung daher nur 2 Nullstellen finden können.

Hinweis: Es kann sein, das aufgrund deiner beschränkten numerischen Genauigkeit du dieses Phänomen nicht siehst. Mit den richtigen (bzw. falschen) Koeffizienten, könnte man aber auch bei einer begrenzten numerischen Genauigkeit auf Zyklen stoßen.

Mein Fazit lautet daher, dass das Verfahren die Probleme des Newtonverfahrens nicht löst und daher keinen (oder nur einen sehr geringen) Mehrwert bietet.

 

[umask007]

Werde das mal testen, danke für dein Feedback.