Problem mit coth(x) Funktion

beamham92

Ich möchte ein Programm schreiben, das mithilfe einer Folge den Cotangens Hyperbolicus berechnet. Habe das Programm soweit fertig bekomme aber verrückte Werte als Ergebnis:

"Bitte geben sie einen Wert fuer x ein: 1

Bitte geben sie einen Wert fuer Epsilon ein: 0.0000001

f(x)=coth(x)= 3.95673e+306
Vergleichswert= 1.31304
Process returned 0 (0x0) execution time : 7.086 s
Press any key to continue.

#include <iostream>
#include <cmath>
using namespace std;"

// Das Programm berechnet den Kotangens Hyperbolicus einer Zahl x, für 0<|x|<pi;
// Zur Berechnung wird die Summenformel in der Form      (((-1)^(n+1)*2^(2*n))/((2*n)!))*B[n]*x^(2*n-1) benutzt;
// Die Berechnung wird solagne durchgeführt bis der Betrag der Differenz zweier;
// aufeinanderfolger Glieder kleiner als ein zudefinierender Abruchwert epsilon            ist;
// oder n den Wert 12 erreicht hat;

    // Funktion zur Potenzberechnung

    float potenz (float basis, float expo)
    {
        if (expo==0){
                    return 1;
                    }

    float pot=basis;

    for(int i=1; i<expo; i++){
            pot=pot*basis;
            }

    return pot;

    }

    // Funktion zur Fakultätsberechnung

    int fakultaet (int f)
{
    if (f==0)
    {
        return 1;
    }

    int fak = 1;

    for (int i=1;i<=f;i++)
    {
        fak = i*fak;
    }
        return fak;
}

int main()
{
    // Variabeldeklination

    double pote1, pote2, pote3, exp1, exp2, exp3, x, B[12], pi, epsilon,  ergebnis, zaehler, nenner, z1, z2;
    int faku, n;

    pi = 3.14159265;
    B[1] = 1.0/6;
    B[2] = 1.0/30;
    B[3] = 1.0/42;
    B[4] = 1.0/30;
    B[5] = 5.0/66;
    B[6] = 691.0/2730;
    B[7] = 7.0/6;
    B[8] = 3617.0/510;
    B[9] = 43867.0/798;
    B[10] = 174611.0/330;
    B[11] = 854513.0/138;
    B[12] = 236364091/2730;
    ergebnis = 1/x;
    z1 = 0;
    n = 1;
    zaehler = 0;
    nenner = 1;

    //Abfrage von x

    cout << "Bitte geben sie einen Wert fuer x ein: ";
    cin >> x;
    cout << '\n';

    if (fabs(x)<=0 || fabs(x)>=pi)
    {
        cout << "Ihr Wert liegt nicht im Definitionsbereich!" << '\n';
        cout << "Bitte geben sie einen neuen Wert fuer x ein: ";
        cin >> x;
        cout << '\n';
    }

    // Abfrage von Epsilon
    // Je kleiner Epsilon gewählt wird, desto genauer ist das Ergebnis

    cout << "Bitte geben sie einen Wert fuer Epsilon ein: ";
    cin >> epsilon;
    cout << '\n';

    if (epsilon<=0)
    {
        cout << "Ihr Wert liegt nicht im Definitionsbereich!" << '\n';
        cout << "Bitte geben sie einen neuen Wert fuer Epsilon ein: ";
        cin >> epsilon;
        cout << '\n';
    }

    do {

        //Berechnung des Zaehlers
        exp1 = n+1;
        exp2 = 2*n;
        exp3 = 2*n-1;
        pote1 = potenz(-1,exp1);
        pote2 = potenz(2,exp2);
        pote3 = potenz(x,exp3);
        zaehler = pote1 * pote2 * pote3 * B[n];

        //Berechnung des Nenners
        faku = 2*n;
        nenner = fakultaet(faku);

        //Summandenberechnung

        z2 = z1;
        z1 = zaehler/nenner;

        //Berechnung des Ergebnisses

        ergebnis = ergebnis + z1;

        //Zählvariable n erhöhen

        n++;
    }

    // Prüfung des Abbruchwertes
    while (fabs(fabs(z2)-fabs(z1)) > epsilon && n<12) ;

    //Ausgabe des Ergbenisses

    cout << "f(x)=coth(x)= " << ergebnis << '\n';

    //Ausgabe des Vergleichswertes aus Mathebibliothek mit coth(x)=cosh(x)/sinh(x)

    double coth = cosh(x)/sinh(x);
    cout << "Vergleichswert= " << coth ;

    }

Könnt ihr mir vielleicht weiterhelfen, ich finde den Fehler nicht.

Freundliche Grüße

SeppJ

Das riecht doch sehr nach Überläufen. Warum benutzt du nicht pow statt einer eigenen Potenzfunktion? Warum benutzt du float? Eine Fakultätsfunktion die int zurück gibt, wird in der Regel bei 13! überlaufen.

beamham92

Ich darf aus der cmath nur fabs und die trigonometrischen Funktionen verwenden. Ich dachte float, da ich Zahlen mit Nachkommastellen ausrechne und float weniger Speicher verbraucht. Wäre double besser? Und anstatt int, long int?

SeppJ

beamham92 schrieb:

Ich darf aus der cmath nur fabs und die trigonometrischen Funktionen verwenden. Ich dachte float, da ich Zahlen mit Nachkommastellen ausrechne und float weniger Speicher verbraucht. Wäre double besser? Und anstatt int, long int?

Das wäre ein Anfang, aber du hast das eigentliche Problem glaube ich nicht verstanden. Die Datentypen zu ändern verschiebt das Problem nur. Wähle einen Algorithmus, bei dem diese Genauigkeitsprobleme gar nicht erst auftreten können. Die Zahlenwerte sollten nie zu extrem werden. Das kannst du erreichen, indem du Teilrechnungen sinnvoll zusammenfasst. Schließlich musst du nicht zwischendurch 10^1000 speichern und dann durch 30! teilen, um dann hinterher 0.15 zu erhalten. Wäre doch besser, wenn du gleich die 0.15 hättest.

Außerdem solltest du bei der Summierung mit den kleinsten Gliedern anfangen und dann zu den größeren gehen. Aus der gewünschten Genauigkeit kannst du schon im Voraus abschätzen, wo du anfangen musst.

Und falls es dir erlaubt ist, würde ich nicht die Taylorentwicklung benutzen. Die Taylorreihe ist was für die Theorie, um damit Beweise zu machen. Zum praktischen Ausrechnen von Zahlen konvergiert sie denkbar schlecht. Außerdem ist das Auftauchen der großen Zahlen zwischendurch numerisch problematisch, wie du hier schon gemerkt hast.

Werner Salomon

Hallo beamham92,

Willkommen im C++-Forum.

Ich unterstelle mal, dass das die Laurent-Reihe sein soll. Wie SeppJ schon gesagt hat, finden bei Deiner Berechnung Überläufe statt, die das Ergebnis ad absurdum führen. Man vermeidet das, indem man Zähler und Nenner in dem Bruch gemeinsam berechnet, dann werden die Zahlen nicht mehr so groß - siehe factor im Code unten.
Auch solltest Du die eigentliche Berechnung in eine Funktion fassen.

Ich glaube auch, dass die Reihe nur im Bereich von -PI/2 bis +PI/2 taugt. Die Gegend um 0 ist sowieso ausgeschlossen.

#include <iostream>
#include <cmath> // std::abs, std::tanh

namespace
{
    namespace bernoulli
    {
        // siehe <http://mathworld.wolfram.com/BernoulliNumber.html>
        const double numbers2[] = {
            1.,         // B_0
            1./6,       // B_2
            -1./30,     // B_4
            1./42,      // B_6
            -1./30,     // B_8
            5./66,      // B_10
            -691./2730, // B_12
            7./6,       // B_14
            -3617./510, // B_16
            43867./789, // B_18
            -174611./330, // B_20
            854513./138 // B_22
        };
    }
}

// --   Kotangens Hyperbolicus Berechnung nach der 'Laurent Reihe'
//      siehe <http://mathworld.wolfram.com/HyperbolicCotangent.html>
double coth( double x, double eps = 1.E-8 )
{
    double sum = 1/x;
    const double xx = x*x;
    double factor = 4.0*x/2; // 2^(2n) * x^(2*n-1) / (2n)! fuer n=1
    for( int n=1; n < sizeof(bernoulli::numbers2)/sizeof(*bernoulli::numbers2); ++n, factor *= 4*xx/((2*n-1)*2*n) )
    {
        const double delta = factor * bernoulli::numbers2[n];
        sum += delta;
        if( std::abs( delta ) <= eps )
            break;
    }
    return sum;
}

int main()
{
    using namespace std;
    for( double x; cout << "> ", cin >> x; )
        cout << "coth(" << x << ")= " << coth(x) << "  ref:" << 1/std::tanh(x) << endl;
    return 0;
}

~~{Ob das so stimmt, das weiß ich nicht. Die Ergebnisse liegen auch im 'harmlosen Bereich' zu weit von 1/tanh(x) entfernt. Entweder wird delta nicht stetig kleiner oder da steckt noch ein Fehler drin.~~

Gruß
Werner

@Edit: Fehler korrigiert - siehe nächstes Posting

Werner Salomon

Werner Salomon schrieb:

.. oder da steckt noch ein Fehler drin.

Ja - ich hatte mich gestern Abend vertan. Statt

for( int n=1; n < sizeof(bernoulli::numbers2)/sizeof(*bernoulli::numbers2); ++n, factor *= (4*xx/(n-1))/n )

muss es

for( int n=1; n < sizeof(bernoulli::numbers2)/sizeof(*bernoulli::numbers2); ++n, factor *= 4*xx/((2*n-1)*2*n) )

heißen.
Dann bekommt man auch vernünftige Werte - zumindest im Bereich bis ca. PI/2.

Gruß
Werner

beamham92

Ich danke euch beiden für eure Hilfe.
Habe verstanden was ihr meintet und das Programm nochmal mit einer for-Schleife geschrieben, jetzt funktioniert Alles wie gewollt.