[prev in list] [next in list] [prev in thread] [next in thread] 

List:       kde-commits
Subject:    [labplot/analysis_playground] src/backend/nsl: implemented Simpson rule in nsl int
From:       Stefan Gerlach <stefan.gerlach () uni-konstanz ! de>
Date:       2016-09-29 19:35:55
Message-ID: E1bph7b-0004hH-Tp () code ! kde ! org
[Download RAW message or body]

Git commit 4105a2105db73fc28cd79a95dff51dc181be84c0 by Stefan Gerlach.
Committed on 29/09/2016 at 19:35.
Pushed by sgerlach into branch 'analysis_playground'.

implemented Simpson rule in nsl int

M  +3    -3    src/backend/nsl/nsl_int.c
M  +2    -1    src/backend/nsl/nsl_int_test.c
M  +7    -2    src/backend/nsl/nsl_sf_poly.c

http://commits.kde.org/labplot/4105a2105db73fc28cd79a95dff51dc181be84c0

diff --git a/src/backend/nsl/nsl_int.c b/src/backend/nsl/nsl_int.c
index b8edde3..a8e7abe 100644
--- a/src/backend/nsl/nsl_int.c
+++ b/src/backend/nsl/nsl_int.c
@@ -93,11 +93,11 @@ size_t nsl_int_simpson(double *x, double *y, const size_t n, int abs) {
 	for (i = 0; i < n-2; i+=2) {
 		for (j=0; j < 3; j++)
 			xdata[j] = x[i+j], ydata[j] = y[i+j];
-		y[np] = sum;
-		x[np++] = x[i+1];
-		printf("i/sum: %zu-%zu %g\n", i, i+2, sum);
 
 		sum += nsl_sf_poly_interp_lagrange_2_int(xdata, ydata);
+		y[np] = sum;
+		x[np++] = x[i+1];
+		/*printf("i/sum: %zu-%zu %g\n", i, i+2, sum);*/
 	}
 	y[n-1] = sum;
 
diff --git a/src/backend/nsl/nsl_int_test.c b/src/backend/nsl/nsl_int_test.c
index bbeb691..be0c4ae 100644
--- a/src/backend/nsl/nsl_int_test.c
+++ b/src/backend/nsl/nsl_int_test.c
@@ -76,12 +76,13 @@ int main() {
 	puts("");
 
 	printf("integral (Simpson's 1/3, 3-point):\n");
+	/*double ydata5[]={2,2,2,2,2};*/
 	double ydata5[]={1,2,3,-1,-3};
 	size_t np = nsl_int_simpson(xdata, ydata5, n, 0);
 
 	for (i=0; i < np; i++)
 		printf("%g %g\n", xdata[i], ydata5[i]);
-	printf("sum = %g (n = %zu)\n", ydata5[n-1], np);
+	printf("sum = %g (n = %zu)\n", ydata5[np-1], np);
 	puts("");
 
 	return 0;
diff --git a/src/backend/nsl/nsl_sf_poly.c b/src/backend/nsl/nsl_sf_poly.c
index b362145..a7828f9 100644
--- a/src/backend/nsl/nsl_sf_poly.c
+++ b/src/backend/nsl/nsl_sf_poly.c
@@ -106,6 +106,7 @@ double complex nsl_sf_poly_reversed_bessel_theta(int n, double complex x) {
 /***************** interpolating polynomials *************/
 
 double nsl_sf_poly_interp_lagrange_0_int(double *x, double y) {
+	/* rectangle rule */
 	return (x[1]-x[0])*y;
 }
 
@@ -116,6 +117,7 @@ double nsl_sf_poly_interp_lagrange_1_deriv(double *x, double *y) {
 	return (y[0]-y[1])/(x[1]-x[0]);
 }
 double nsl_sf_poly_interp_lagrange_1_int(double *x, double *y) {
+	/* trapezoid rule (2-point) */
 	return (x[1]-x[0])*(y[0]+y[1])/2.;
 }
 double nsl_sf_poly_interp_lagrange_1_absint(double *x, double *y) {
@@ -147,8 +149,11 @@ double nsl_sf_poly_interp_lagrange_2_deriv2(double *x, double *y) {
 	return 2.*( y[0]/(h1*h12) - y[1]/(h1*h2) + y[2]/(h12*h2) );
 }
 double nsl_sf_poly_interp_lagrange_2_int(double *x, double *y) {
-	/* TODO */
-	return 0;
+	/* Simpson's 1/3 rule for non-uniform grid (3-point) */
+	double h1 = x[1]-x[0], h2 = x[2]-x[1];
+	double h12 = h1+h2, h1_2 = h1/h2;
+
+	return h12/6.*( (2.-1./h1_2)*y[0] + (2.+h1_2+1./h1_2)*y[1] + (2.-h1_2)*y[2] );
 }
 
 double nsl_sf_poly_interp_lagrange_3(double v, double *x, double *y) {
[prev in list] [next in list] [prev in thread] [next in thread]