1 /* interpolation/interp_poly.c
2 *
3 * Copyright (C) 2001 DAN, HO-JIN
4 * Copyright (C) 2013 Patrick Alken
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or (at
9 * your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 */
20
21 /* Modified for standalone use in polynomial directory, B.Gough 2001 */
22
23 #include <config.h>
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_poly.h>
26
27 int
gsl_poly_dd_init(double dd[],const double xa[],const double ya[],size_t size)28 gsl_poly_dd_init (double dd[], const double xa[], const double ya[],
29 size_t size)
30 {
31 size_t i, j;
32
33 /* Newton's divided differences */
34
35 dd[0] = ya[0];
36
37 for (j = size - 1; j >= 1; j--)
38 {
39 dd[j] = (ya[j] - ya[j - 1]) / (xa[j] - xa[j - 1]);
40 }
41
42 for (i = 2; i < size; i++)
43 {
44 for (j = size - 1; j >= i; j--)
45 {
46 dd[j] = (dd[j] - dd[j - 1]) / (xa[j] - xa[j - i]);
47 }
48 }
49
50 return GSL_SUCCESS;
51 }
52
53 int
gsl_poly_dd_taylor(double c[],double xp,const double dd[],const double xa[],size_t size,double w[])54 gsl_poly_dd_taylor (double c[], double xp,
55 const double dd[], const double xa[], size_t size,
56 double w[])
57 {
58 size_t i, j;
59
60 for (i = 0; i < size; i++)
61 {
62 c[i] = 0.0;
63 w[i] = 0.0;
64 }
65
66 w[size - 1] = 1.0;
67
68 c[0] = dd[0];
69
70 for (i = size - 1; i-- > 0;)
71 {
72 w[i] = -w[i + 1] * (xa[size - 2 - i] - xp);
73
74 for (j = i + 1; j < size - 1; j++)
75 {
76 w[j] = w[j] - w[j + 1] * (xa[size - 2 - i] - xp);
77 }
78
79 for (j = i; j < size; j++)
80 {
81 c[j - i] += w[j] * dd[size - i - 1];
82 }
83 }
84
85 return GSL_SUCCESS;
86 }
87
88 /*
89 gsl_poly_dd_hermite_init()
90 Compute divided difference representation of data
91 for Hermite polynomial interpolation
92
93 Inputs: dd - (output) array of size 2*size containing
94 divided differences, dd[k] = f[z_0,z_1,...,z_k]
95 za - (output) array of size 2*size containing
96 z values
97 xa - x data
98 ya - y data
99 dya - dy/dx data
100 size - size of xa,ya,dya arrays
101
102 Return: success
103 */
104
105 int
gsl_poly_dd_hermite_init(double dd[],double za[],const double xa[],const double ya[],const double dya[],const size_t size)106 gsl_poly_dd_hermite_init (double dd[], double za[], const double xa[], const double ya[],
107 const double dya[], const size_t size)
108 {
109 const size_t N = 2 * size;
110 size_t i, j;
111
112 /* Hermite divided differences */
113
114 dd[0] = ya[0];
115
116 /* compute: dd[j] = f[z_{j-1},z_j] for j \in [1,N-1] */
117 for (j = 0; j < size; ++j)
118 {
119 za[2*j] = xa[j];
120 za[2*j + 1] = xa[j];
121
122 if (j != 0)
123 {
124 dd[2*j] = (ya[j] - ya[j - 1]) / (xa[j] - xa[j - 1]);
125 dd[2*j - 1] = dya[j - 1];
126 }
127 }
128
129 dd[N - 1] = dya[size - 1];
130
131 for (i = 2; i < N; i++)
132 {
133 for (j = N - 1; j >= i; j--)
134 {
135 dd[j] = (dd[j] - dd[j - 1]) / (za[j] - za[j - i]);
136 }
137 }
138
139 return GSL_SUCCESS;
140 } /* gsl_poly_dd_hermite_init() */
141