File indexing completed on 2024-05-12 15:27:06
0001 /*************************************************************************** 0002 File : nsl_diff.c 0003 Project : LabPlot 0004 Description : NSL numerical differentiation functions 0005 -------------------------------------------------------------------- 0006 Copyright : (C) 2016 by Stefan Gerlach (stefan.gerlach@uni.kn) 0007 ***************************************************************************/ 0008 0009 /*************************************************************************** 0010 * * 0011 * This program is free software; you can redistribute it and/or modify * 0012 * it under the terms of the GNU General Public License as published by * 0013 * the Free Software Foundation; either version 2 of the License, or * 0014 * (at your option) any later version. * 0015 * * 0016 * This program is distributed in the hope that it will be useful, * 0017 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 0018 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 0019 * GNU General Public License for more details. * 0020 * * 0021 * You should have received a copy of the GNU General Public License * 0022 * along with this program; if not, write to the Free Software * 0023 * Foundation, Inc., 51 Franklin Street, Fifth Floor, * 0024 * Boston, MA 02110-1301 USA * 0025 * * 0026 ***************************************************************************/ 0027 0028 /* TODO: 0029 * add more orders 0030 */ 0031 0032 #include "nsl_diff.h" 0033 #include "nsl_common.h" 0034 #include "nsl_sf_poly.h" 0035 0036 const char* nsl_diff_deriv_order_name[] = {i18n("first"), i18n("second"), i18n("third"), i18n("fourth"), i18n("fifth"), i18n("sixth")}; 0037 0038 double nsl_diff_first_central(double xm, double fm, double xp, double fp) { 0039 return (fp - fm)/(xp - xm); 0040 } 0041 0042 int nsl_diff_deriv_first_equal(const double *x, double *y, const size_t n) { 0043 if (n < 3) 0044 return -1; 0045 0046 double dy=0, oldy=0, oldy2=0; 0047 size_t i; 0048 for (i=0; i < n; i++) { 0049 if (i == 0) /* forward */ 0050 dy = (-y[2] + 4.*y[1] - 3.*y[0])/(x[2]-x[0]); 0051 else if (i == n-1) { /* backward */ 0052 y[i] = (3.*y[i] - 4.*y[i-1] + y[i-2])/(x[i]-x[i-2]); 0053 y[i-1] = oldy; 0054 } 0055 else 0056 dy = (y[i+1]-y[i-1])/(x[i+1]-x[i-1]); 0057 0058 if (i > 1) 0059 y[i-2] = oldy2; 0060 if (i > 0 && i < n-1) 0061 oldy2 = oldy; 0062 0063 oldy = dy; 0064 } 0065 0066 return 0; 0067 } 0068 0069 int nsl_diff_first_deriv(const double *x, double *y, const size_t n, int order) { 0070 switch (order) { 0071 case 2: 0072 return nsl_diff_first_deriv_second_order(x, y, n); 0073 case 4: 0074 return nsl_diff_first_deriv_fourth_order(x, y, n); 0075 /*TODO: higher order */ 0076 default: 0077 printf("nsl_diff_first_deriv() unsupported order %d\n", order); 0078 return -1; 0079 } 0080 } 0081 0082 int nsl_diff_first_deriv_second_order(const double *x, double *y, const size_t n) { 0083 if (n < 3) 0084 return -1; 0085 0086 double dy=0, oldy=0, oldy2=0, xdata[3], ydata[3]; 0087 size_t i, j; 0088 for (i=0; i < n; i++) { 0089 if (i == 0) { 0090 /* 3-point forward */ 0091 for (j=0; j < 3; j++) 0092 xdata[j]=x[j], ydata[j]=y[j]; 0093 dy = nsl_sf_poly_interp_lagrange_2_deriv(x[0], xdata, ydata); 0094 } else if (i < n-1) { 0095 /* 3-point center */ 0096 for (j=0; j < 3; j++) 0097 xdata[j]=x[i-1+j], ydata[j]=y[i-1+j]; 0098 dy = nsl_sf_poly_interp_lagrange_2_deriv(x[i], xdata, ydata); 0099 } else if (i == n-1) { 0100 /* 3-point backward */ 0101 y[i] = nsl_sf_poly_interp_lagrange_2_deriv(x[i], xdata, ydata); 0102 y[i-1] = oldy; 0103 } 0104 0105 if (i > 1) 0106 y[i-2] = oldy2; 0107 if (i > 0 && i < n-1) 0108 oldy2 = oldy; 0109 0110 oldy = dy; 0111 } 0112 0113 return 0; 0114 } 0115 0116 int nsl_diff_first_deriv_fourth_order(const double *x, double *y, const size_t n) { 0117 if (n < 5) 0118 return -1; 0119 0120 double dy[5]={0}, xdata[5], ydata[5]; 0121 size_t i, j; 0122 for (i=0; i < n; i++) { 0123 if (i == 0) 0124 for (j=0; j < 5; j++) 0125 xdata[j]=x[j], ydata[j]=y[j]; 0126 else if (i > 1 && i < n-2) 0127 for (j=0; j < 5; j++) 0128 xdata[j]=x[i-2+j], ydata[j]=y[i-2+j]; 0129 0130 /* 5-point rule */ 0131 dy[0] = nsl_sf_poly_interp_lagrange_4_deriv(x[i], xdata, ydata); 0132 0133 if (i == n-1) 0134 for (j=0; j < 4; j++) 0135 y[i-j] = dy[j]; 0136 0137 if (i > 3) 0138 y[i-4] = dy[4]; 0139 for (j=4; j > 0; j--) 0140 if (i >= j-1) 0141 dy[j] = dy[j-1]; 0142 } 0143 0144 return 0; 0145 } 0146 0147 int nsl_diff_first_deriv_avg(const double *x, double *y, const size_t n) { 0148 if (n < 1) 0149 return -1; 0150 0151 size_t i; 0152 double dy=0, oldy=0; 0153 for (i=0; i<n; i++) { 0154 if(i == 0) 0155 dy = (y[1]-y[0])/(x[1]-x[0]); 0156 else if (i == n-1) 0157 y[i] = (y[i]-y[i-1])/(x[i]-x[i-1]); 0158 else 0159 dy = ( (y[i+1]-y[i])/(x[i+1]-x[i]) + (y[i]-y[i-1])/(x[i]-x[i-1]) )/2.; 0160 0161 if (i > 0) 0162 y[i-1] = oldy; 0163 0164 oldy = dy; 0165 } 0166 0167 return 0; 0168 } 0169 0170 int nsl_diff_second_deriv(const double *x, double *y, const size_t n, int order) { 0171 switch (order) { 0172 case 1: 0173 return nsl_diff_second_deriv_first_order(x, y, n); 0174 case 2: 0175 return nsl_diff_second_deriv_second_order(x, y, n); 0176 case 3: 0177 return nsl_diff_second_deriv_third_order(x, y, n); 0178 /*TODO: higher order */ 0179 default: 0180 printf("nsl_diff_second_deriv() unsupported order %d\n", order); 0181 return -1; 0182 } 0183 } 0184 0185 int nsl_diff_second_deriv_first_order(const double *x, double *y, const size_t n) { 0186 if (n < 3) 0187 return -1; 0188 0189 double dy[3]={0}, xdata[3], ydata[3]; 0190 size_t i, j; 0191 for (i=0; i<n; i++) { 0192 if (i == 0) 0193 for (j=0; j < 3; j++) 0194 xdata[j]=x[j], ydata[j]=y[j]; 0195 else if (i > 1 && i < n-1) 0196 for (j=0; j < 3; j++) 0197 xdata[j]=x[i-1+j], ydata[j]=y[i-1+j]; 0198 0199 /* 3-point rule */ 0200 dy[0] = nsl_sf_poly_interp_lagrange_2_deriv2(xdata, ydata); 0201 0202 if (i == n-1) { 0203 y[i] = dy[0]; 0204 y[i-1] = dy[1]; 0205 } 0206 0207 if (i > 1) 0208 y[i-2] = dy[2]; 0209 if (i > 0) 0210 dy[2] = dy[1]; 0211 0212 dy[1] = dy[0]; 0213 } 0214 0215 return 0; 0216 } 0217 0218 int nsl_diff_second_deriv_second_order(const double *x, double *y, const size_t n) { 0219 if (n < 4) 0220 return -1; 0221 0222 double dy[4]={0}, xdata[4], ydata[4]; 0223 size_t i, j; 0224 for (i=0; i<n; i++) { 0225 if (i == 0) { 0226 /* 4-point forward */ 0227 for (j=0; j < 4; j++) 0228 xdata[j]=x[j], ydata[j]=y[j]; 0229 dy[0] = nsl_sf_poly_interp_lagrange_3_deriv2(x[i], xdata, ydata); 0230 } 0231 else if (i == n-1) { 0232 /* 4-point backward */ 0233 for (j=0; j < 4; j++) 0234 xdata[j]=x[i-3+j], ydata[j]=y[i-3+j]; 0235 y[i] = nsl_sf_poly_interp_lagrange_3_deriv2(x[i], xdata, ydata); 0236 y[i-1] = dy[1]; 0237 y[i-2] = dy[2]; 0238 } 0239 else { 0240 /* 3-point center */ 0241 for (j=0; j < 3; j++) 0242 xdata[j]=x[i-1+j], ydata[j]=y[i-1+j]; 0243 dy[0] = nsl_sf_poly_interp_lagrange_2_deriv2(xdata, ydata); 0244 } 0245 0246 if (i > 2) 0247 y[i-3] = dy[3]; 0248 for (j=3; j > 0; j--) 0249 if (i >= j-1) 0250 dy[j] = dy[j-1]; 0251 } 0252 0253 return 0; 0254 } 0255 0256 int nsl_diff_second_deriv_third_order(const double *x, double *y, const size_t n) { 0257 if (n < 5) 0258 return -1; 0259 0260 double dy[5]={0}, xdata[5], ydata[5]; 0261 size_t i, j; 0262 for (i=0; i < n; i++) { 0263 if (i == 0) 0264 for (j=0; j < 5; j++) 0265 xdata[j]=x[j], ydata[j]=y[j]; 0266 else if (i > 2 && i < n-3) 0267 for (j=0; j < 5; j++) 0268 xdata[j]=x[i-2+j], ydata[j]=y[i-2+j]; 0269 0270 /* 5-point rule */ 0271 dy[0] = nsl_sf_poly_interp_lagrange_4_deriv2(x[i], xdata, ydata); 0272 0273 if (i == n-1) 0274 for (j=0; j < 4; j++) 0275 y[i-j] = dy[j]; 0276 0277 if (i > 3) 0278 y[i-4] = dy[4]; 0279 for (j=4; j > 0; j--) 0280 if (i >= j-1) 0281 dy[j] = dy[j-1]; 0282 } 0283 0284 return 0; 0285 } 0286 0287 int nsl_diff_third_deriv(const double *x, double *y, const size_t n, int order) { 0288 switch (order) { 0289 case 2: 0290 return nsl_diff_third_deriv_second_order(x, y, n); 0291 /*TODO: higher order */ 0292 default: 0293 printf("nsl_diff_third_deriv() unsupported order %d\n", order); 0294 return -1; 0295 } 0296 } 0297 0298 int nsl_diff_third_deriv_second_order(const double *x, double *y, const size_t n) { 0299 if (n < 5) 0300 return -1; 0301 0302 double dy[5]={0}, xdata[5], ydata[5]; 0303 size_t i, j; 0304 for (i=0; i < n; i++) { 0305 if (i == 0) 0306 for (j=0; j < 5; j++) 0307 xdata[j]=x[j], ydata[j]=y[j]; 0308 else if (i > 2 && i < n-3) 0309 for (j=0; j < 5; j++) 0310 xdata[j]=x[i-2+j], ydata[j]=y[i-2+j]; 0311 0312 /* 5-point rule */ 0313 dy[0] = nsl_sf_poly_interp_lagrange_4_deriv3(x[i], xdata, ydata); 0314 0315 if (i == n-1) 0316 for (j=0; j < 4; j++) 0317 y[i-j] = dy[j]; 0318 0319 if (i > 3) 0320 y[i-4] = dy[4]; 0321 for (j=4; j > 0; j--) 0322 if (i >= j-1) 0323 dy[j] = dy[j-1]; 0324 } 0325 0326 return 0; 0327 } 0328 0329 int nsl_diff_fourth_deriv(const double *x, double *y, const size_t n, int order) { 0330 switch (order) { 0331 case 1: 0332 return nsl_diff_fourth_deriv_first_order(x, y, n); 0333 case 3: 0334 return nsl_diff_fourth_deriv_third_order(x, y, n); 0335 /*TODO: higher order */ 0336 default: 0337 printf("nsl_diff_fourth_deriv() unsupported order %d\n", order); 0338 return -1; 0339 } 0340 } 0341 0342 int nsl_diff_fourth_deriv_first_order(const double *x, double *y, const size_t n) { 0343 if (n < 5) 0344 return -1; 0345 0346 double dy[5]={0}, xdata[5], ydata[5]; 0347 size_t i, j; 0348 for (i=0; i < n; i++) { 0349 if (i == 0) 0350 for (j=0; j < 5; j++) 0351 xdata[j]=x[j], ydata[j]=y[j]; 0352 else if (i > 2 && i < n-3) 0353 for (j=0; j < 5; j++) 0354 xdata[j]=x[i-2+j], ydata[j]=y[i-2+j]; 0355 0356 /* 5-point rule */ 0357 dy[0] = nsl_sf_poly_interp_lagrange_4_deriv4(xdata, ydata); 0358 0359 if (i == n-1) 0360 for (j=0; j < 4; j++) 0361 y[i-j] = dy[j]; 0362 0363 if (i > 3) 0364 y[i-4] = dy[4]; 0365 for (j=4; j > 0; j--) 0366 if (i >= j-1) 0367 dy[j] = dy[j-1]; 0368 } 0369 0370 return 0; 0371 } 0372 0373 int nsl_diff_fourth_deriv_third_order(const double *x, double *y, const size_t n) { 0374 if (n < 7) 0375 return -1; 0376 0377 double dy[7]={0}, xdata[7], ydata[7]; 0378 size_t i, j; 0379 for (i=0; i < n; i++) { 0380 if (i == 0) 0381 for (j=0; j < 7; j++) 0382 xdata[j]=x[j], ydata[j]=y[j]; 0383 else if (i > 3 && i < n-4) 0384 for (j=0; j < 7; j++) 0385 xdata[j]=x[i-3+j], ydata[j]=y[i-3+j]; 0386 0387 /* 7-point rule */ 0388 dy[0] = nsl_sf_poly_interp_lagrange_6_deriv4(x[i], xdata, ydata); 0389 0390 if (i == n-1) 0391 for (j=0; j < 6; j++) 0392 y[i-j] = dy[j]; 0393 0394 if (i > 5) 0395 y[i-6] = dy[6]; 0396 for (j=6; j > 0; j--) 0397 if (i >= j-1) 0398 dy[j] = dy[j-1]; 0399 } 0400 0401 return 0; 0402 } 0403 0404 int nsl_diff_fifth_deriv(const double *x, double *y, const size_t n, int order) { 0405 switch (order) { 0406 case 2: 0407 return nsl_diff_fifth_deriv_second_order(x, y, n); 0408 /*TODO: higher order */ 0409 default: 0410 printf("nsl_diff_fifth_deriv() unsupported order %d\n", order); 0411 return -1; 0412 } 0413 } 0414 0415 int nsl_diff_fifth_deriv_second_order(const double *x, double *y, const size_t n) { 0416 if (n < 7) 0417 return -1; 0418 0419 double dy[7]={0}, xdata[7], ydata[7]; 0420 size_t i, j; 0421 for (i=0; i < n; i++) { 0422 if (i == 0) 0423 for (j=0; j < 7; j++) 0424 xdata[j]=x[j], ydata[j]=y[j]; 0425 else if (i > 3 && i < n-4) 0426 for (j=0; j < 7; j++) 0427 xdata[j]=x[i-3+j], ydata[j]=y[i-3+j]; 0428 0429 /* 7-point rule */ 0430 dy[0] = nsl_sf_poly_interp_lagrange_6_deriv5(x[i], xdata, ydata); 0431 0432 if (i == n-1) 0433 for (j=0; j < 6; j++) 0434 y[i-j] = dy[j]; 0435 0436 if (i > 5) 0437 y[i-6] = dy[6]; 0438 for (j=6; j > 0; j--) 0439 if (i >= j-1) 0440 dy[j] = dy[j-1]; 0441 } 0442 0443 return 0; 0444 } 0445 0446 int nsl_diff_sixth_deriv(const double *x, double *y, const size_t n, int order) { 0447 switch (order) { 0448 case 1: 0449 return nsl_diff_sixth_deriv_first_order(x, y, n); 0450 /*TODO: higher order */ 0451 default: 0452 printf("nsl_diff_sixth_deriv() unsupported order %d\n", order); 0453 return -1; 0454 } 0455 } 0456 0457 int nsl_diff_sixth_deriv_first_order(const double *x, double *y, const size_t n) { 0458 if (n < 7) 0459 return -1; 0460 0461 double dy[7]={0}, xdata[7], ydata[7]; 0462 size_t i, j; 0463 for (i=0; i < n; i++) { 0464 if (i == 0) 0465 for (j=0; j < 7; j++) 0466 xdata[j]=x[j], ydata[j]=y[j]; 0467 else if (i > 3 && i < n-4) 0468 for (j=0; j < 7; j++) 0469 xdata[j]=x[i-3+j], ydata[j]=y[i-3+j]; 0470 0471 /* 7-point rule */ 0472 dy[0] = nsl_sf_poly_interp_lagrange_6_deriv6(xdata, ydata); 0473 0474 if (i == n-1) 0475 for (j=0; j < 6; j++) 0476 y[i-j] = dy[j]; 0477 0478 if (i > 5) 0479 y[i-6] = dy[6]; 0480 for (j=6; j > 0; j--) 0481 if (i >= j-1) 0482 dy[j] = dy[j-1]; 0483 } 0484 0485 return 0; 0486 } 0487