Numerische Integration Rechteckregel - Genauigkeitsproblem
-
Erstmal hallo an alle hier, ich bin neu aber ich hoffe ihr helft mir trotzdem.
Folgendes Problem:
Ich habe ein Programm zur numerischen Integration geschrieben. Es benutzt die Rechteckregel, teilt also die Kurve in Abschnitte der Größe h, berechnet den Funktionswert in der Mitte jedes Abschnitts und somit die approximierte Fläche unter der Kurve. Das ganz mit verschieden großer Anzahl an Abschnitten.Zu berechnen ist die Fläche zwischen den Kurven y² = 9x und y = x² - 4x +6. Die ist genau 5, Es geht jetzt darum wie genau die Berechnung der numerischen Integration mit der Rechteckregel ist. Erwartungsgemäß wird sie mit größerer Anzahl an Abschnitten besser, leider habe ich jedoch den Effekt dass sie wieder schlechter wird. Habt ihr eine Idee woran das liegt? Ist in meinem Code etwas falsch? Oder läuft da irgendwo eine Variable über?
Vielen Dank für eure Mühe!
Hier der etwas gekürzte Quellecode (es ist ein MPI Programm).
#include <stdio.h> #include <stdlib.h> #include <mpi.h> #include <math.h> #include <getopt.h> int main(int argc, char *argv[]) { // define functions double f(double x) { return (sqrtl(9 * x) - (x * x - 4 * x + 6)); } // init variables int myrank, size; // MPI variables double *recvbuf; double low_limit = 1; // global lower x value double up_limit = 4; // global upper x value double exact_integral = 5; // that is what the result should be long intervals; // number of intervals double h; // size of the interval // Init MPI MPI_Init(&argc, &argv); MPI_Comm_rank(MPI_COMM_WORLD, &myrank); MPI_Comm_size(MPI_COMM_WORLD, &size); // minimum intervals intervals = size; // Init private vars recvbuf = malloc(1 * sizeof(double)); double p_low_limit; // private lower limit long p_intervals; // private number of intervals double integral, h_half; // Calculate integral for different interval amounts while (integral != exact_integral) { // init local vars h = (up_limit - low_limit) / intervals; p_intervals = (intervals / size); p_low_limit = low_limit + ((up_limit - low_limit) / size) * myrank; integral = 0; h_half = h * 0.5f; for (long i = 0; i < p_intervals; ++i) { integral += h * f(p_low_limit + (i * h) + h_half); } // sum up everything MPI_Reduce(&integral, recvbuf, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); // print result if (myrank == 0) { printf("integral = %.51f error %.51f intervals %ld \n", recvbuf[0], fabs(recvbuf[0] - exact_integral), intervals); } intervals = intervals * 2; } // Finish MPI MPI_Finalize(); return 0; }
So starte ich das ganze:
mpicc -std=c99 -o KPP_2.out KPP_2.c mpirun -np 4 KPP_2.out -r 0
Und das ist die Ausgabe wo sich zeigt das der absolute Fehler wieder größer wird.
integral = 5.157827293307573057745685218833386898040771484375000 error 0.157827293307573057745685218833386898040771484375000 intervals 4 integral = 5.039525397367471626353108149487525224685668945312500 error 0.039525397367471626353108149487525224685668945312500 intervals 8 integral = 5.009886071709677146657213597791269421577453613281250 error 0.009886071709677146657213597791269421577453613281250 intervals 16 integral = 5.002471821729022494196215120609849691390991210937500 error 0.002471821729022494196215120609849691390991210937500 intervals 32 integral = 5.000617974565890300198134355014190077781677246093750 error 0.000617974565890300198134355014190077781677246093750 intervals 64 integral = 5.000154494839655683335877256467938423156738281250000 error 0.000154494839655683335877256467938423156738281250000 intervals 128 integral = 5.000038623784837099606193078216165304183959960937500 error 0.000038623784837099606193078216165304183959960937500 intervals 256 integral = 5.000009655950892195619417179841548204421997070312500 error 0.000009655950892195619417179841548204421997070312500 intervals 512 integral = 5.000002413988014815515725786099210381507873535156250 error 0.000002413988014815515725786099210381507873535156250 intervals 1024 integral = 5.000000603497022133581140224123373627662658691406250 error 0.000000603497022133581140224123373627662658691406250 intervals 2048 integral = 5.000000150874259752242778631625697016716003417968750 error 0.000000150874259752242778631625697016716003417968750 intervals 4096 integral = 5.000000037718559831034781382186338305473327636718750 error 0.000000037718559831034781382186338305473327636718750 intervals 8192 integral = 5.000000009429641956160139670828357338905334472656250 error 0.000000009429641956160139670828357338905334472656250 intervals 16384 integral = 5.000000002357408490638590592425316572189331054687500 error 0.000000002357408490638590592425316572189331054687500 intervals 32768 integral = 5.000000000589357007640955998795107007026672363281250 error 0.000000000589357007640955998795107007026672363281250 intervals 65536 integral = 5.000000000147362122504546277923509478569030761718750 error 0.000000000147362122504546277923509478569030761718750 intervals 131072 integral = 5.000000000036841640849161194637417793273925781250000 error 0.000000000036841640849161194637417793273925781250000 intervals 262144 integral = 5.000000000009232614672782801790162920951843261718750 error 0.000000000009232614672782801790162920951843261718750 intervals 524288 integral = 5.000000000002331468351712828734889626502990722656250 error 0.000000000002331468351712828734889626502990722656250 intervals 1048576 integral = 5.000000000000611954931173386285081505775451660156250 error 0.000000000000611954931173386285081505775451660156250 intervals 2097152 integral = 5.000000000000202504679691628552973270416259765625000 error 0.000000000000202504679691628552973270416259765625000 intervals 4194304 integral = 5.000000000000203392858111328678205609321594238281250 error 0.000000000000203392858111328678205609321594238281250 intervals 8388608 integral = 5.000000000000133226762955018784850835800170898437500 error 0.000000000000133226762955018784850835800170898437500 intervals 16777216 integral = 5.000000000000155431223447521915659308433532714843750 error 0.000000000000155431223447521915659308433532714843750 intervals 33554432 integral = 4.999999999999778843573494668817147612571716308593750 error 0.000000000000221156426505331182852387428283691406250 intervals 67108864 integral = 5.000000000000778932474077009828761219978332519531250 error 0.000000000000778932474077009828761219978332519531250 intervals 134217728 integral = 4.999999999998994582028899458236992359161376953125000 error 0.000000000001005417971100541763007640838623046875000 intervals 268435456 integral = 5.000000000001032063323691545519977807998657226562500 error 0.000000000001032063323691545519977807998657226562500 intervals 536870912 integral = 4.999999999994840571559961972525343298912048339843750 error 0.000000000005159428440038027474656701087951660156250 intervals 1073741824 integral = 5.000000000009165113112885592272505164146423339843750 error 0.000000000009165113112885592272505164146423339843750 intervals 2147483648 integral = 5.000000000031284308477097511058673262596130371093750 error 0.000000000031284308477097511058673262596130371093750 intervals 4294967296 integral = 4.999999999851998389033269631909206509590148925781250 error 0.000000000148001610966730368090793490409851074218750 intervals 8589934592 integral = 4.999999999930170524464756454108282923698425292968750 error 0.000000000069829475535243545891717076301574707031250 intervals 17179869184 integral = 5.000000000279551493065355316502973437309265136718750 error 0.000000000279551493065355316502973437309265136718750 intervals 34359738368 integral = 5.000000002061629977845313987927511334419250488281250 error 0.000000002061629977845313987927511334419250488281250 intervals 68719476736 integral = 5.000000005184316087536444683792069554328918457031250 error 0.000000005184316087536444683792069554328918457031250 intervals 137438953472 integral = 5.000000013307920099236980604473501443862915039062500 error 0.000000013307920099236980604473501443862915039062500 intervals 274877906944 integral = 5.000000033180974412516661686822772026062011718750000 error 0.000000033180974412516661686822772026062011718750000 intervals 549755813888
-
Ich würde sagen, dass die numerischen Fehler, die bei jeder einzelnen Rechteckauswertung passieren (z.B. beim Bestimmen der linken und rechten Grenze) irgendwann durch die schiere Anzahl der Intervalle zusammengenommen die Fehler, die durch die Rechteckapproximation entstehen, dominieren.
-
StephanP schrieb:
ich jedoch den Effekt dass sie wieder schlechter wird. Habt ihr eine Idee woran das liegt? Ist in meinem Code etwas falsch? Oder läuft da irgendwo eine Variable über?
Bei 549755813888 Rechtecken gibt es auch 549755813888 Rundungsfehler. Nur ganz winzige Fehlerchen in der letzten Stelle. Aber sauviele.
Und die Rechteckflächen werden im Vergleich zur bisherigen Flächensumme winzig klein, so daß nur wenige signifikate Stellen übrigbleiben, wenn zum Addieren die winzige Zahl ihren Exponenten der Summe anpaßt.
-
Also wird quasi der relative Fehler je Intervall mit zunehmender Anzahl an Intervallen größer und somit auch der gesamt Fehler. Das leuchtet ein. Vielen Dank!
Könnte ich die maximal erreichbare Genauigkeit durch einen größeren Datentyp (long double?) erhöhen? (Performance wird natürlich schlechter)
-
Wenn bei deinem Compiler long double genauer ist als double, warum nicht? Probier's doch einfach aus.
-
StephanP schrieb:
Könnte ich die maximal erreichbare Genauigkeit durch einen größeren Datentyp (long double?) erhöhen? (Performance wird natürlich schlechter)
Ja, aber kaum sinnvoll. Wenn Du numerisch was verbessern willst, addiere immer Zahlen, die ähnlich groß sind oder zuerst die kleinen und dann die großen.
Und wechsle das Verfahren. Nimm die Trapezregel. Dann wird's schon mit weniger Stützstellen recht genau. Oder noch besser die Simpson-Regel. Die ist echt ein Kracher, was die Genauigkeit angeht. Bei Deiner Parabel landet sie zufällig sofort einen Volltreffer. Aber zu anderen Funktionen ist sie auch sehr gutmütig. Und sie ist einfach zu programmieren.
-
volkard schrieb:
StephanP schrieb:
Könnte ich die maximal erreichbare Genauigkeit durch einen größeren Datentyp (long double?) erhöhen? (Performance wird natürlich schlechter)
Ja, aber kaum sinnvoll. Wenn Du numerisch was verbessern willst, addiere immer Zahlen, die ähnlich groß sind oder zuerst die kleinen und dann die großen.
Und wechsle das Verfahren. Nimm die Trapezregel. Dann wird's schon mit weniger Stützstellen recht genau. Oder noch besser die Simpson-Regel. Die ist echt ein Kracher, was die Genauigkeit angeht. Bei Deiner Parabel landet sie zufällig sofort einen Volltreffer. Aber zu anderen Funktionen ist sie auch sehr gutmütig. Und sie ist einfach zu programmieren.
Es geht eben um die Untersuchung der Rechteckregel (Relativer Fehler und Performance). Die Trapezregel ist als nächstes dran.
Vielen Dank an alle, ihr habt mir sehr geholfen.