labplot  2.8.2
About: LabPlot is an application for plotting and analysis of 2D and 3D functions and data. It is a complete rewrite of LabPlot1 and lacks in the first release a lot of features available in the predecessor. On the other hand, the GUI and the usability is more superior.
  Fossies Dox: labplot-2.8.2.tar.gz  ("unofficial" and yet experimental doxygen-generated source code documentation)  

nsl_diff.c
Go to the documentation of this file.
1 /***************************************************************************
2  File : nsl_diff.c
3  Project : LabPlot
4  Description : NSL numerical differentiation functions
5  --------------------------------------------------------------------
6  Copyright : (C) 2016 by Stefan Gerlach (stefan.gerlach@uni.kn)
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  * This program is distributed in the hope that it will be useful, *
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
19  * GNU General Public License for more details. *
20  * *
21  * You should have received a copy of the GNU General Public License *
22  * along with this program; if not, write to the Free Software *
23  * Foundation, Inc., 51 Franklin Street, Fifth Floor, *
24  * Boston, MA 02110-1301 USA *
25  * *
26  ***************************************************************************/
27 
28 /* TODO:
29  * add more orders
30 */
31 
32 #include "nsl_diff.h"
33 #include "nsl_common.h"
34 #include "nsl_sf_poly.h"
35 
36 const char* nsl_diff_deriv_order_name[] = {i18n("first"), i18n("second"), i18n("third"), i18n("fourth"), i18n("fifth"), i18n("sixth")};
37 
38 double nsl_diff_first_central(double xm, double fm, double xp, double fp) {
39  return (fp - fm)/(xp - xm);
40 }
41 
42 int nsl_diff_deriv_first_equal(const double *x, double *y, const size_t n) {
43  if (n < 3)
44  return -1;
45 
46  double dy=0, oldy=0, oldy2=0;
47  size_t i;
48  for (i=0; i < n; i++) {
49  if (i == 0) /* forward */
50  dy = (-y[2] + 4.*y[1] - 3.*y[0])/(x[2]-x[0]);
51  else if (i == n-1) { /* backward */
52  y[i] = (3.*y[i] - 4.*y[i-1] + y[i-2])/(x[i]-x[i-2]);
53  y[i-1] = oldy;
54  }
55  else
56  dy = (y[i+1]-y[i-1])/(x[i+1]-x[i-1]);
57 
58  if (i > 1)
59  y[i-2] = oldy2;
60  if (i > 0 && i < n-1)
61  oldy2 = oldy;
62 
63  oldy = dy;
64  }
65 
66  return 0;
67 }
68 
69 int nsl_diff_first_deriv(const double *x, double *y, const size_t n, int order) {
70  switch (order) {
71  case 2:
73  case 4:
75  /*TODO: higher order */
76  default:
77  printf("nsl_diff_first_deriv() unsupported order %d\n", order);
78  return -1;
79  }
80 }
81 
82 int nsl_diff_first_deriv_second_order(const double *x, double *y, const size_t n) {
83  if (n < 3)
84  return -1;
85 
86  double dy=0, oldy=0, oldy2=0, xdata[3], ydata[3];
87  size_t i, j;
88  for (i=0; i < n; i++) {
89  if (i == 0) {
90  /* 3-point forward */
91  for (j=0; j < 3; j++)
92  xdata[j]=x[j], ydata[j]=y[j];
93  dy = nsl_sf_poly_interp_lagrange_2_deriv(x[0], xdata, ydata);
94  } else if (i < n-1) {
95  /* 3-point center */
96  for (j=0; j < 3; j++)
97  xdata[j]=x[i-1+j], ydata[j]=y[i-1+j];
98  dy = nsl_sf_poly_interp_lagrange_2_deriv(x[i], xdata, ydata);
99  } else if (i == n-1) {
100  /* 3-point backward */
101  y[i] = nsl_sf_poly_interp_lagrange_2_deriv(x[i], xdata, ydata);
102  y[i-1] = oldy;
103  }
104 
105  if (i > 1)
106  y[i-2] = oldy2;
107  if (i > 0 && i < n-1)
108  oldy2 = oldy;
109 
110  oldy = dy;
111  }
112 
113  return 0;
114 }
115 
116 int nsl_diff_first_deriv_fourth_order(const double *x, double *y, const size_t n) {
117  if (n < 5)
118  return -1;
119 
120  double dy[5]={0}, xdata[5], ydata[5];
121  size_t i, j;
122  for (i=0; i < n; i++) {
123  if (i == 0)
124  for (j=0; j < 5; j++)
125  xdata[j]=x[j], ydata[j]=y[j];
126  else if (i > 1 && i < n-2)
127  for (j=0; j < 5; j++)
128  xdata[j]=x[i-2+j], ydata[j]=y[i-2+j];
129 
130  /* 5-point rule */
131  dy[0] = nsl_sf_poly_interp_lagrange_4_deriv(x[i], xdata, ydata);
132 
133  if (i == n-1)
134  for (j=0; j < 4; j++)
135  y[i-j] = dy[j];
136 
137  if (i > 3)
138  y[i-4] = dy[4];
139  for (j=4; j > 0; j--)
140  if (i >= j-1)
141  dy[j] = dy[j-1];
142  }
143 
144  return 0;
145 }
146 
147 int nsl_diff_first_deriv_avg(const double *x, double *y, const size_t n) {
148  if (n < 1)
149  return -1;
150 
151  size_t i;
152  double dy=0, oldy=0;
153  for (i=0; i<n; i++) {
154  if(i == 0)
155  dy = (y[1]-y[0])/(x[1]-x[0]);
156  else if (i == n-1)
157  y[i] = (y[i]-y[i-1])/(x[i]-x[i-1]);
158  else
159  dy = ( (y[i+1]-y[i])/(x[i+1]-x[i]) + (y[i]-y[i-1])/(x[i]-x[i-1]) )/2.;
160 
161  if (i > 0)
162  y[i-1] = oldy;
163 
164  oldy = dy;
165  }
166 
167  return 0;
168 }
169 
170 int nsl_diff_second_deriv(const double *x, double *y, const size_t n, int order) {
171  switch (order) {
172  case 1:
173  return nsl_diff_second_deriv_first_order(x, y, n);
174  case 2:
176  case 3:
177  return nsl_diff_second_deriv_third_order(x, y, n);
178  /*TODO: higher order */
179  default:
180  printf("nsl_diff_second_deriv() unsupported order %d\n", order);
181  return -1;
182  }
183 }
184 
185 int nsl_diff_second_deriv_first_order(const double *x, double *y, const size_t n) {
186  if (n < 3)
187  return -1;
188 
189  double dy[3]={0}, xdata[3], ydata[3];
190  size_t i, j;
191  for (i=0; i<n; i++) {
192  if (i == 0)
193  for (j=0; j < 3; j++)
194  xdata[j]=x[j], ydata[j]=y[j];
195  else if (i > 1 && i < n-1)
196  for (j=0; j < 3; j++)
197  xdata[j]=x[i-1+j], ydata[j]=y[i-1+j];
198 
199  /* 3-point rule */
200  dy[0] = nsl_sf_poly_interp_lagrange_2_deriv2(xdata, ydata);
201 
202  if (i == n-1) {
203  y[i] = dy[0];
204  y[i-1] = dy[1];
205  }
206 
207  if (i > 1)
208  y[i-2] = dy[2];
209  if (i > 0)
210  dy[2] = dy[1];
211 
212  dy[1] = dy[0];
213  }
214 
215  return 0;
216 }
217 
218 int nsl_diff_second_deriv_second_order(const double *x, double *y, const size_t n) {
219  if (n < 4)
220  return -1;
221 
222  double dy[4]={0}, xdata[4], ydata[4];
223  size_t i, j;
224  for (i=0; i<n; i++) {
225  if (i == 0) {
226  /* 4-point forward */
227  for (j=0; j < 4; j++)
228  xdata[j]=x[j], ydata[j]=y[j];
229  dy[0] = nsl_sf_poly_interp_lagrange_3_deriv2(x[i], xdata, ydata);
230  }
231  else if (i == n-1) {
232  /* 4-point backward */
233  for (j=0; j < 4; j++)
234  xdata[j]=x[i-3+j], ydata[j]=y[i-3+j];
235  y[i] = nsl_sf_poly_interp_lagrange_3_deriv2(x[i], xdata, ydata);
236  y[i-1] = dy[1];
237  y[i-2] = dy[2];
238  }
239  else {
240  /* 3-point center */
241  for (j=0; j < 3; j++)
242  xdata[j]=x[i-1+j], ydata[j]=y[i-1+j];
243  dy[0] = nsl_sf_poly_interp_lagrange_2_deriv2(xdata, ydata);
244  }
245 
246  if (i > 2)
247  y[i-3] = dy[3];
248  for (j=3; j > 0; j--)
249  if (i >= j-1)
250  dy[j] = dy[j-1];
251  }
252 
253  return 0;
254 }
255 
256 int nsl_diff_second_deriv_third_order(const double *x, double *y, const size_t n) {
257  if (n < 5)
258  return -1;
259 
260  double dy[5]={0}, xdata[5], ydata[5];
261  size_t i, j;
262  for (i=0; i < n; i++) {
263  if (i == 0)
264  for (j=0; j < 5; j++)
265  xdata[j]=x[j], ydata[j]=y[j];
266  else if (i > 2 && i < n-3)
267  for (j=0; j < 5; j++)
268  xdata[j]=x[i-2+j], ydata[j]=y[i-2+j];
269 
270  /* 5-point rule */
271  dy[0] = nsl_sf_poly_interp_lagrange_4_deriv2(x[i], xdata, ydata);
272 
273  if (i == n-1)
274  for (j=0; j < 4; j++)
275  y[i-j] = dy[j];
276 
277  if (i > 3)
278  y[i-4] = dy[4];
279  for (j=4; j > 0; j--)
280  if (i >= j-1)
281  dy[j] = dy[j-1];
282  }
283 
284  return 0;
285 }
286 
287 int nsl_diff_third_deriv(const double *x, double *y, const size_t n, int order) {
288  switch (order) {
289  case 2:
290  return nsl_diff_third_deriv_second_order(x, y, n);
291  /*TODO: higher order */
292  default:
293  printf("nsl_diff_third_deriv() unsupported order %d\n", order);
294  return -1;
295  }
296 }
297 
298 int nsl_diff_third_deriv_second_order(const double *x, double *y, const size_t n) {
299  if (n < 5)
300  return -1;
301 
302  double dy[5]={0}, xdata[5], ydata[5];
303  size_t i, j;
304  for (i=0; i < n; i++) {
305  if (i == 0)
306  for (j=0; j < 5; j++)
307  xdata[j]=x[j], ydata[j]=y[j];
308  else if (i > 2 && i < n-3)
309  for (j=0; j < 5; j++)
310  xdata[j]=x[i-2+j], ydata[j]=y[i-2+j];
311 
312  /* 5-point rule */
313  dy[0] = nsl_sf_poly_interp_lagrange_4_deriv3(x[i], xdata, ydata);
314 
315  if (i == n-1)
316  for (j=0; j < 4; j++)
317  y[i-j] = dy[j];
318 
319  if (i > 3)
320  y[i-4] = dy[4];
321  for (j=4; j > 0; j--)
322  if (i >= j-1)
323  dy[j] = dy[j-1];
324  }
325 
326  return 0;
327 }
328 
329 int nsl_diff_fourth_deriv(const double *x, double *y, const size_t n, int order) {
330  switch (order) {
331  case 1:
332  return nsl_diff_fourth_deriv_first_order(x, y, n);
333  case 3:
334  return nsl_diff_fourth_deriv_third_order(x, y, n);
335  /*TODO: higher order */
336  default:
337  printf("nsl_diff_fourth_deriv() unsupported order %d\n", order);
338  return -1;
339  }
340 }
341 
342 int nsl_diff_fourth_deriv_first_order(const double *x, double *y, const size_t n) {
343  if (n < 5)
344  return -1;
345 
346  double dy[5]={0}, xdata[5], ydata[5];
347  size_t i, j;
348  for (i=0; i < n; i++) {
349  if (i == 0)
350  for (j=0; j < 5; j++)
351  xdata[j]=x[j], ydata[j]=y[j];
352  else if (i > 2 && i < n-3)
353  for (j=0; j < 5; j++)
354  xdata[j]=x[i-2+j], ydata[j]=y[i-2+j];
355 
356  /* 5-point rule */
357  dy[0] = nsl_sf_poly_interp_lagrange_4_deriv4(xdata, ydata);
358 
359  if (i == n-1)
360  for (j=0; j < 4; j++)
361  y[i-j] = dy[j];
362 
363  if (i > 3)
364  y[i-4] = dy[4];
365  for (j=4; j > 0; j--)
366  if (i >= j-1)
367  dy[j] = dy[j-1];
368  }
369 
370  return 0;
371 }
372 
373 int nsl_diff_fourth_deriv_third_order(const double *x, double *y, const size_t n) {
374  if (n < 7)
375  return -1;
376 
377  double dy[7]={0}, xdata[7], ydata[7];
378  size_t i, j;
379  for (i=0; i < n; i++) {
380  if (i == 0)
381  for (j=0; j < 7; j++)
382  xdata[j]=x[j], ydata[j]=y[j];
383  else if (i > 3 && i < n-4)
384  for (j=0; j < 7; j++)
385  xdata[j]=x[i-3+j], ydata[j]=y[i-3+j];
386 
387  /* 7-point rule */
388  dy[0] = nsl_sf_poly_interp_lagrange_6_deriv4(x[i], xdata, ydata);
389 
390  if (i == n-1)
391  for (j=0; j < 6; j++)
392  y[i-j] = dy[j];
393 
394  if (i > 5)
395  y[i-6] = dy[6];
396  for (j=6; j > 0; j--)
397  if (i >= j-1)
398  dy[j] = dy[j-1];
399  }
400 
401  return 0;
402 }
403 
404 int nsl_diff_fifth_deriv(const double *x, double *y, const size_t n, int order) {
405  switch (order) {
406  case 2:
407  return nsl_diff_fifth_deriv_second_order(x, y, n);
408  /*TODO: higher order */
409  default:
410  printf("nsl_diff_fifth_deriv() unsupported order %d\n", order);
411  return -1;
412  }
413 }
414 
415 int nsl_diff_fifth_deriv_second_order(const double *x, double *y, const size_t n) {
416  if (n < 7)
417  return -1;
418 
419  double dy[7]={0}, xdata[7], ydata[7];
420  size_t i, j;
421  for (i=0; i < n; i++) {
422  if (i == 0)
423  for (j=0; j < 7; j++)
424  xdata[j]=x[j], ydata[j]=y[j];
425  else if (i > 3 && i < n-4)
426  for (j=0; j < 7; j++)
427  xdata[j]=x[i-3+j], ydata[j]=y[i-3+j];
428 
429  /* 7-point rule */
430  dy[0] = nsl_sf_poly_interp_lagrange_6_deriv5(x[i], xdata, ydata);
431 
432  if (i == n-1)
433  for (j=0; j < 6; j++)
434  y[i-j] = dy[j];
435 
436  if (i > 5)
437  y[i-6] = dy[6];
438  for (j=6; j > 0; j--)
439  if (i >= j-1)
440  dy[j] = dy[j-1];
441  }
442 
443  return 0;
444 }
445 
446 int nsl_diff_sixth_deriv(const double *x, double *y, const size_t n, int order) {
447  switch (order) {
448  case 1:
449  return nsl_diff_sixth_deriv_first_order(x, y, n);
450  /*TODO: higher order */
451  default:
452  printf("nsl_diff_sixth_deriv() unsupported order %d\n", order);
453  return -1;
454  }
455 }
456 
457 int nsl_diff_sixth_deriv_first_order(const double *x, double *y, const size_t n) {
458  if (n < 7)
459  return -1;
460 
461  double dy[7]={0}, xdata[7], ydata[7];
462  size_t i, j;
463  for (i=0; i < n; i++) {
464  if (i == 0)
465  for (j=0; j < 7; j++)
466  xdata[j]=x[j], ydata[j]=y[j];
467  else if (i > 3 && i < n-4)
468  for (j=0; j < 7; j++)
469  xdata[j]=x[i-3+j], ydata[j]=y[i-3+j];
470 
471  /* 7-point rule */
472  dy[0] = nsl_sf_poly_interp_lagrange_6_deriv6(xdata, ydata);
473 
474  if (i == n-1)
475  for (j=0; j < 6; j++)
476  y[i-j] = dy[j];
477 
478  if (i > 5)
479  y[i-6] = dy[6];
480  for (j=6; j > 0; j--)
481  if (i >= j-1)
482  dy[j] = dy[j-1];
483  }
484 
485  return 0;
486 }
487 
#define i18n(m)
Definition: nsl_common.h:38
int nsl_diff_sixth_deriv(const double *x, double *y, const size_t n, int order)
Definition: nsl_diff.c:446
int nsl_diff_fourth_deriv_first_order(const double *x, double *y, const size_t n)
Definition: nsl_diff.c:342
int nsl_diff_fifth_deriv_second_order(const double *x, double *y, const size_t n)
Definition: nsl_diff.c:415
int nsl_diff_third_deriv(const double *x, double *y, const size_t n, int order)
Definition: nsl_diff.c:287
int nsl_diff_second_deriv_first_order(const double *x, double *y, const size_t n)
Definition: nsl_diff.c:185
int nsl_diff_second_deriv(const double *x, double *y, const size_t n, int order)
Definition: nsl_diff.c:170
int nsl_diff_fourth_deriv_third_order(const double *x, double *y, const size_t n)
Definition: nsl_diff.c:373
int nsl_diff_second_deriv_third_order(const double *x, double *y, const size_t n)
Definition: nsl_diff.c:256
int nsl_diff_first_deriv_fourth_order(const double *x, double *y, const size_t n)
Definition: nsl_diff.c:116
int nsl_diff_fourth_deriv(const double *x, double *y, const size_t n, int order)
Definition: nsl_diff.c:329
int nsl_diff_deriv_first_equal(const double *x, double *y, const size_t n)
Definition: nsl_diff.c:42
int nsl_diff_first_deriv(const double *x, double *y, const size_t n, int order)
Definition: nsl_diff.c:69
double nsl_diff_first_central(double xm, double fm, double xp, double fp)
Definition: nsl_diff.c:38
int nsl_diff_first_deriv_second_order(const double *x, double *y, const size_t n)
Definition: nsl_diff.c:82
int nsl_diff_third_deriv_second_order(const double *x, double *y, const size_t n)
Definition: nsl_diff.c:298
const char * nsl_diff_deriv_order_name[]
Definition: nsl_diff.c:36
int nsl_diff_sixth_deriv_first_order(const double *x, double *y, const size_t n)
Definition: nsl_diff.c:457
int nsl_diff_second_deriv_second_order(const double *x, double *y, const size_t n)
Definition: nsl_diff.c:218
int nsl_diff_first_deriv_avg(const double *x, double *y, const size_t n)
Definition: nsl_diff.c:147
int nsl_diff_fifth_deriv(const double *x, double *y, const size_t n, int order)
Definition: nsl_diff.c:404
double nsl_sf_poly_interp_lagrange_6_deriv5(double v, double *x, double *y)
Definition: nsl_sf_poly.c:278
double nsl_sf_poly_interp_lagrange_4_deriv(double v, double *x, double *y)
Definition: nsl_sf_poly.c:218
double nsl_sf_poly_interp_lagrange_3_deriv2(double v, double *x, double *y)
Definition: nsl_sf_poly.c:176
double nsl_sf_poly_interp_lagrange_4_deriv4(double *x, double *y)
Definition: nsl_sf_poly.c:251
double nsl_sf_poly_interp_lagrange_2_deriv2(double *x, double *y)
Definition: nsl_sf_poly.c:146
double nsl_sf_poly_interp_lagrange_2_deriv(double v, double *x, double *y)
Definition: nsl_sf_poly.c:140
double nsl_sf_poly_interp_lagrange_6_deriv4(double v, double *x, double *y)
Definition: nsl_sf_poly.c:263
double nsl_sf_poly_interp_lagrange_4_deriv2(double v, double *x, double *y)
Definition: nsl_sf_poly.c:233
double nsl_sf_poly_interp_lagrange_6_deriv6(double *x, double *y)
Definition: nsl_sf_poly.c:289
double nsl_sf_poly_interp_lagrange_4_deriv3(double v, double *x, double *y)
Definition: nsl_sf_poly.c:244