2016-09-12 2 views
1

Ich schrieb ein kurzes Programm in C, um lineare Interpolation durchzuführen, die iteriert, um eine Wurzel einer Funktion zu einer gegebenen Anzahl von Dezimalpunkten zu geben. Bisher für eine Funktion f (x):Warum gibt C Python für diesen Algorithmus verschiedene Werte?

long double f(long double x) { 
     return (pow(x, 3) + 2 * x - 8); 
    } 

Das Programm konvergiert nicht zu einem 1dp-Wert. Das Programm aktualisiert Variablen a und b, zwischen denen die Wurzel von f (x) liegt, bis a und b beide auf die gleiche Zahl mit einer gegebenen Genauigkeit runden. Mit langen Doppelzimmer und die obige Funktion, zeigt ein Debugger für die ersten 2 Wiederholungen:

a = 1.5555555555555556 
    a = 1.6444444444444444 

obwohl es hätte sein sollen:

a = 1.5555555555555556 
    a = 1.653104925053533 

Das Programm schlägt fehl Werte danach zu aktualisieren. Die Gleichung für die lineare Interpolation, die ich verwende, ist eine umgeordnete Version des mathematischen here und der Code, den ich benutze, ist eine C-Version eines Python-Programms, das ich geschrieben habe. Warum erhält die C-Implementierung trotz des gleichen Algorithmus unterschiedliche Werte und wie kann ich sie beheben?

OK Ich bin immer noch den Dreh dieser bekommen, aber hoffentlich unten Ich habe einen Minimal, Complete, und prüfbare Beispiel:

#include <stdlib.h> 
#include <stdio.h> 
#include <math.h> 

long double a; long double b; long double c; // The values for interpolation 
long double fofa; long double fofb; long double fofc; // The values f(a), f(b) and f(c) 
const int dp = 1; // The number of decimal places to be accurate to 

long double f(long double x) { 
    return (pow(x, 3) + 2 * x - 8); 
} 

int main(void) { 
    a = 1; b = 2; 
    while(roundf(a * pow(10, dp))/pow(10, dp) != roundf(b * pow(10, dp))/pow(10, dp)) { // While a and b don't round to the same number... 

     fofa = f(a); fofb = f(b); // Resolve the functions 
     printf("So f(a) = %g, f(b) = %g\n", (double)fofa, (double)fofb); // Print the results 

     c = (b * abs(fofa) + a * abs(fofb))/(abs(fofb) + abs(fofa)); // Linear Interpolation 

     fofc = f(c); 

     if(fofc < 0) { 
      a = c; 
     } 
     else if(fofc == 0) { 
      a = c; 
      break; 
     } 
     else { 
      b = c; 
     } 

    } 

    printf("The root is %g, to %d decimal places, after %d iterations.\n", (double)a, dp, i); 
} 
+0

Sie können nicht einige Debug verwenden, um welche Werte zu sehen Sie bekommen für 'fofa' (was ist' f (a) ', nehme ich an),' fofb', und alle anderen Zwischenwerte? – Teepeemm

+0

Was verwenden Sie als Wert für "dp" (1, 10, 15, etwas anderes)? Warum benutzt du 'float'-precision mit' roundf() ', wenn du' long double' für 'a' und' fofa' usw. hast? (Sie sollten 'roundl()' und 'powl()' (anstelle von 'pow()') verwenden, es sei denn, Sie verwenden '', aber Sie sollten erwähnen, dass, wenn Sie es sind und Sie nicht würden benutze 'roundf()', aber nur 'round()' ... –

+1

Bitte stelle ein [mcve], kein Puzzlespiel zur Verfügung, siehe auch [ask] – Olaf

Antwort

7

Die Funktion abs() (von <stdlib.h>) die Unterschrift int abs(int); - Sie sind Ganzzahlen aus Ihren Berechnungen erhalten.

Sie sollten long double fabsl(long double); von <math.h> verwenden.

Sie sollten auch powl() statt pow() (long double vs double) verwenden, und roundl() statt roundf() (long double vs float).

Stellen Sie sicher, dass Sie die richtigen Typen verwenden, mit anderen Worten.


Wenn Sie die Typ Probleme behoben haben, haben Sie immer noch ein Problem mit der Konvergenz. Es wäre hilfreich gewesen, wenn Sie einen MCVE (Minimal, Complete, Verifiable Example) produziert hätten. Dies ist jedoch ein MCVE, die ich aus Ihrer Frage ableiten kann:

#include <math.h> 
#include <stdio.h> 

static inline long double f(long double x) 
{ 
    return(powl(x, 3) + 2 * x - 8); 
} 

int main(void) 
{ 
    long double a = 1.0L; 
    long double b = 2.0L; 
    int dp = 6; 

    while (roundl(a * powl(10, dp))/powl(10, dp) != roundl(b * powl(10, dp))/powl(10, dp)) 
    { 
     long double fofa = f(a); 
     long double fofb = f(b); 
     long double c = (b * fabsl(fofa) + a * fabsl(fofb))/(fabsl(fofb) + fabsl(fofa)); 
     long double fofc = f(c); 
     printf("a = %+.10Lf, f(a) = %+.10Lf\n", a, fofa); 
     printf("b = %+.10Lf, f(b) = %+.10Lf\n", b, fofb); 
     printf("c = %+.10Lf, f(c) = %+.10Lf\n", c, fofc); 
     putchar('\n'); 

     if (fofc < 0.0L) 
     { 
      a = c; 
     } 
     else if (fofc == 0.0L) 
     { 
      a = c; 
      break; 
     } 
     else 
     { 
      b = c; 
     } 
    } 
    printf("Result: a = %Lg\n", a); 
    return 0; 
} 

Der Ausgang ich von ihm erhalten ist:

a = +1.0000000000, f(a) = -5.0000000000 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.5555555556, f(c) = -1.1248285322 

a = +1.5555555556, f(a) = -1.1248285322 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6531049251, f(c) = -0.1762579238 

a = +1.6531049251, f(a) = -0.1762579238 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6677455452, f(c) = -0.0258828049 

a = +1.6677455452, f(a) = -0.0258828049 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6698816424, f(c) = -0.0037639074 

a = +1.6698816424, f(a) = -0.0037639074 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6701919841, f(c) = -0.0005465735 

a = +1.6701919841, f(a) = -0.0005465735 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702370440, f(c) = -0.0000793539 

a = +1.6702370440, f(a) = -0.0000793539 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702435859, f(c) = -0.0000115206 

a = +1.6702435859, f(a) = -0.0000115206 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702445357, f(c) = -0.0000016726 

a = +1.6702445357, f(a) = -0.0000016726 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446735, f(c) = -0.0000002428 

a = +1.6702446735, f(a) = -0.0000002428 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446936, f(c) = -0.0000000353 

a = +1.6702446936, f(a) = -0.0000000353 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446965, f(c) = -0.0000000051 

a = +1.6702446965, f(a) = -0.0000000051 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446969, f(c) = -0.0000000007 

a = +1.6702446969, f(a) = -0.0000000007 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446970, f(c) = -0.0000000001 

a = +1.6702446970, f(a) = -0.0000000001 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446970, f(c) = -0.0000000000 

a = +1.6702446970, f(a) = -0.0000000000 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446970, f(c) = -0.0000000000 

a = +1.6702446970, f(a) = -0.0000000000 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446970, f(c) = -0.0000000000 

a = +1.6702446970, f(a) = -0.0000000000 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446970, f(c) = -0.0000000000 

Der Grund für die Endlos-Schleife ist klar; Der Unterschied zwischen a und b ist kein kleiner Bruchteil. Sie müssen Ihr Konvergenzkriterium überprüfen. Es sollte wahrscheinlich fofc mit 0.0 innerhalb der angegebenen Anzahl von Dezimalstellen vergleichen - oder etwas in dieser Richtung.

+0

Vielen Dank, überraschend konvergiert es jetzt (zu viele dp, aber ich werde das herausfinden!) – AWSM

+0

Den Hinweisen von [Lutzl] folgen (http : //stackoverflow.com/users/3088138/lutzl) [Antwort] (http://stackoverflow.com/a/39459981), finden Sie möglicherweise die Wikipedia articl e auf der [Falsche Positionsmethode] (https://en.wikipedia.org/wiki/False_position_method) erleuchtend. Ich habe nicht versucht, die Mathematik weiter zu verbessern. –

3

Was Sie implementieren Regula falsi genannt wird oder falsche Position Methode.

Es ist eigentlich nicht notwendig, die absoluten Werte zu verwenden, solange die Bedingung des umgekehrten Vorzeichens beibehalten wird.

Es gibt ein sehr bekanntes Problem mit Plain Vanilla Regula Falsi Problem, in dem Moment, wenn die Funktion über das verbleibende Intervall konvex ist, wird ein Endpunkt nicht weiter in Richtung der Wurzel bewegt. Es gibt einfache Modifikationen, um das zu vermeiden, zum Beispiel das Einfügen von Bisektionsschritten. Noch einfacher zu implementieren, aber schwerer zu verstehen ist die Illinois-Modifikation. Siehe Wikipedia Artikel für Regula Falsi für Details.

Oder diese Frage und Antwort: Regula-Falsi Algorithm?


Angepasst von Antwort in Verbindung:

#include<math.h> 
#include<stdio.h> 

long double f(long double x) { 
    return powl(x, 3) + 2 * x - 8; 
} 

int main(void) { 
    const int dp = 18; 
    long double eps=0.5*powl(10,-dp); 
    int i=0; 
    long double a=1, fofa = f(a); 
    long double b=2, fofb = f(b); 

    printf("\na=%.21Lf b=%.21Lf fofa=%.21Lf fofb=%.21Lf\n------\n",a,b, fofa,fofb); 

    if(signbit(fofb)==signbit(fofa)) { 
     printf("Warning, initial values do not have opposite sign!\n"); 
    } 
    do { 
     long double c=(a*fofb-b*fofa)/(fofb-fofa), fofc = f(c); 

     if(signbit(fofc)!=signbit(fofa)) { 
      b=a; fofb=fofa; 
      a=c; fofa=fofc; 
     } else { 
      a=c; fofa=fofc; 
      fofb *= 0.5; 
     } 
     i++; 
     printf("\na=c=%.21Lf b=%.21Lf fofa=fofc=%.21Lf fofb=%.21Lf",c,b, fofc,fofb); 

    } while(fabsl(b-a)>eps); 

    printf("\ngoal reached after %d iterations\n",i); 

    return 0; 
} 

mit

a=1.000000000000000000000 b=2.000000000000000000000 fofa=-5.000000000000000000000 fofb=4.000000000000000000000 
------ 

a=c=1.555555555555555555507 b=2.000000000000000000000 fofa=fofc=-1.124828532235939643549 fofb=2.000000000000000000000 
a=c=1.715539947322212467064 b=1.555555555555555555507 fofa=fofc=0.480046589479470829469 fofb=-1.124828532235939643549 
a=c=1.667685780603345490963 b=1.715539947322212467064 fofa=fofc=-0.026500999000164700194 fofb=0.480046589479470829469 
a=c=1.670189362207942139265 b=1.715539947322212467064 fofa=fofc=-0.000573759143326624515 fofb=0.240023294739735414734 
a=c=1.670297511133477909220 b=1.670189362207942139265 fofa=fofc=0.000547652143260468627 fofb=-0.000573759143326624515 
a=c=1.670244695550498326532 b=1.670297511133477909220 fofa=fofc=-0.000000014643676336194 fofb=0.000547652143260468627 
a=c=1.670244696962696986627 b=1.670297511133477909220 fofa=fofc=-0.000000000000373724489 fofb=0.000273826071630234313 
a=c=1.670244696962769068724 b=1.670244696962696986627 fofa=fofc=0.000000000000373706274 fofb=-0.000000000000373724489 
a=c=1.670244696962733028543 b=1.670244696962769068724 fofa=fofc=-0.000000000000000000434 fofb=0.000000000000373706274 
a=c=1.670244696962733028651 b=1.670244696962733028543 fofa=fofc=0.000000000000000000867 fofb=-0.000000000000000000434 
goal reached after 10 iterations 
Verwandte Themen