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.


Anmelden zum Antworten