Implementierung Brent-Decker-Verfahren



  • Hi,

    ich habe gerade die Implementierung des Brent-Decker-Verfahrens
    von matlab https://de.wikipedia.org/wiki/Brent-Verfahren#Algorithmus_von_Brent_f.C3.BCr_Matlab
    nach C++ übertagen.

    double a = funktionswerte.at(i).get_x();
    			double fa = fct->at(a);
    			double b = funktionswerte.at(i + 2).get_x();
    			double fb = fct->at(b);
    
    			if ((fb*fa) > 0)
    				throw std::exception("f(a) and f(b) should have different key signature");
    			double c = a;
    			double fc = fa;
    			double d = b - a; 
    			double e = d; 
    			int iter = 0;
    
    			while (iter < 1000) {
    			iter++;
    				if ((fb*fc) > 0)
    					c = a; fc = fa; d = b - a; e = d;
    				if (std::abs(fc) < std::abs(fb)) {
    					a = b;
    					b = c;
    					c = a;
    					fa = fb;
    					fb = fc;
    					fc = fa;
    				}
    				double t = std::pow(10,-20);
    				double tol = 2 * std::numeric_limits<double>::epsilon() * std::abs(b) + t;
    				double m = (c - b) / 2;
    				if ((std::abs(m) > tol) && (std::abs(fb) > 0)) {
    					if ((std::abs(e) < tol) || (std::abs(fa) <= std::abs(fb))) {
    						d = m;
    						e = m;
    					}
    					else {
    						double s = fb / fa;
    						double p;
    						double q;
    						double r;
    						if (a == c) {
    							p = 2 * m*s;
    							q = 1 - s;
    						}
    						else {
    							q = fa / fc;
    							r = fb / fc;
    							p = s *(2 * m*q*(q - r) - (b - a)*(r - 1));
    							q = (q - 1)*(r - 1)*(s - 1);
    						}
    						if (p > 0)
    							q *= -1;
    						else
    							p *= -1;
    						s = e;
    						e = d;
    						if ((2 * p < 3 * m*q - std::abs(tol*q)) && (p < std::abs(s*q / 2))) {
    							d = p / q;
    						}
    						else {
    							d = m;
    							e = m;
    						}
    					}
    					a = b;
    					fa = fb;
    					if (std::abs(d) > tol)
    						b += d;
    					else {
    						if (m > 0)
    							b += tol;
    						else
    							b -= tol;
    					}
    				} 
    				else{
    					break;
    				}
    				fb = fct->at(b);
    			}
    			double xs = b;
    			nst.push_back(xs);
    		}
    

    Das Problem dabei ist, dass das Verfahren mir, bei Polynomfunktionen, die Nullstellen nicht präzise liefert. Abweichungn von 0,5 bis 1 sind die Regel. Wende ich das Verfahren jetzt allerdings auf eine cos(x) Funktion an, werden mir Nullstellen bis auf viele Nachkommastellen exact geliefer.
    Was stimmt jetzt hier also nicht? Ist die Implementierung Fehlerhaft, oder meine Herangehensweise.

    PS: Als Eingangswerte verwende ich einen Wert vor dem Vorzeichenwechsel und einen nach dem Vorzeichenwechsel



  • if ((fb*fc) > 0)
    {
      c = a; fc = fa; d = b - a; e = d;
    }
    

    Und es dauerte keine 5 Minuten bis ich den Fehler fand. 🙂 Überprüfe nochmals ob dein Code stimmt.



  • Und es dauerte über fünf Minuten bis ich drauf gekommen bin, was du mir sagen willst. Natürlich hab ich die Klammern vergessen... Man sollte nicht nur Copy&Paste machen 🙄
    Danke! 🙂


  • Mod

    axels. schrieb:

    Und es dauerte über fünf Minuten bis ich drauf gekommen bin, was du mir sagen willst.

    Was wohl daran liegt, dass er den Code nicht zitiert hat, sondern seine korrigierte Variante präsentiert hat. Das ist ziemlich fies.


Log in to reply