Qucs-core  0.0.18
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
spline.cpp
Go to the documentation of this file.
1 /*
2  * spline.cpp - spline class implementation
3  *
4  * Copyright (C) 2005, 2006, 2009 Stefan Jahn <stefan@lkcc.org>
5  *
6  * This 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 2, or (at your option)
9  * any later version.
10  *
11  * This software is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this package; see the file COPYING. If not, write to
18  * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19  * Boston, MA 02110-1301, USA.
20  *
21  * $Id$
22  *
23  */
24 
25 #if HAVE_CONFIG_H
26 # include <config.h>
27 #endif
28 
29 #include <stdio.h>
30 #include <stdlib.h>
31 #include <string.h>
32 #include <assert.h>
33 
34 #include "logging.h"
35 #include "complex.h"
36 #include "object.h"
37 #include "vector.h"
38 #include "poly.h"
39 #include "tvector.h"
40 #include "tridiag.h"
41 #include "spline.h"
42 
43 namespace qucs {
44 
45 // Constructor creates an instance of the spline class.
47  x = f0 = f1 = f2 = f3 = NULL;
48  d0 = dn = 0;
49  n = 0;
51 }
52 
53 // Constructor creates an instance of the spline class with given boundary.
54 spline::spline (int b) {
55  x = f0 = f1 = f2 = f3 = NULL;
56  d0 = dn = 0;
57  n = 0;
58  boundary = b;
59 }
60 
61 // Constructor creates an instance of the spline class with vector data given.
63  x = f0 = f1 = f2 = f3 = NULL;
64  d0 = dn = 0;
65  n = 0;
67  vectors (y, t);
68  construct ();
69 }
70 
71 // Constructor creates an instance of the spline class with tvector data given.
73  x = f0 = f1 = f2 = f3 = NULL;
74  d0 = dn = 0;
75  n = 0;
77  vectors (y, t);
78  construct ();
79 }
80 
81 #define t_ (t)
82 #define y_ (y)
83 
84 // Pass interpolation datapoints as vectors.
86  int i = t.getSize ();
87  assert (y.getSize () == i && i >= 3);
88 
89  // create local copy of f(x)
90  realloc (i);
91  for (i = 0; i <= n; i++) {
92  f0[i] = real (y_(i)); x[i] = real (t_(i));
93  }
94 }
95 
96 // Pass interpolation datapoints as tvectors.
98  int i = t.getSize ();
99  assert (y.getSize () == i && i >= 3);
100 
101  // create local copy of f(x)
102  realloc (i);
103  for (i = 0; i <= n; i++) {
104  f0[i] = y_(i); x[i] = t_(i);
105  }
106 }
107 
108 // Pass interpolation datapoints as pointers.
109 void spline::vectors (nr_double_t * y, nr_double_t * t, int len) {
110  int i = len;
111  assert (i >= 3);
112 
113  // create local copy of f(x)
114  realloc (i);
115  for (i = 0; i <= n; i++) {
116  f0[i] = y[i]; x[i] = t[i];
117  }
118 }
119 
120 // Reallocate vector data if necessary.
121 void spline::realloc (int size) {
122  if (n != size - 1) {
123  n = size - 1;
124  if (f0) delete[] f0;
125  f0 = new nr_double_t[n+1];
126  if (x) delete[] x;
127  x = new nr_double_t[n+1];
128  }
129  if (f1) { delete[] f1; f1 = NULL; }
130  if (f2) { delete[] f2; f2 = NULL; }
131  if (f3) { delete[] f3; f3 = NULL; }
132 }
133 
134 // Construct cubic spline interpolation coefficients.
135 void spline::construct (void) {
136 
137  // calculate first derivative h = f'(x)
138  int i;
139  nr_double_t * h = new nr_double_t[n+1];
140  for (i = 0; i < n; i++) {
141  h[i] = x[i+1] - x[i];
142  if (h[i] == 0.0) {
143  logprint (LOG_ERROR, "ERROR: Duplicate points in spline: %g, %g\n",
144  x[i], x[i+1]);
145  }
146  }
147 
148  // first kind of cubic splines
150 
151  // setup right hand side
152  nr_double_t * b = new nr_double_t[n+1]; // b as in Ax = b
153  for (i = 1; i < n; i++) {
154  nr_double_t _n = f0[i+1] * h[i-1] - f0[i] * (h[i] + h[i-1]) +
155  f0[i-1] * h[i];
156  nr_double_t _d = h[i-1] * h[i];
157  b[i] = 3 * _n / _d;
158  }
159  if (boundary == SPLINE_BC_NATURAL) {
160  // natural boundary condition
161  b[0] = 0;
162  b[n] = 0;
163  } else if (boundary == SPLINE_BC_CLAMPED) {
164  // start and end derivatives given
165  b[0] = 3 * ((f0[1] - f0[0]) / h[0] - d0);
166  b[n] = 3 * (dn - (f0[n] - f0[n-1]) / h[n-1]);
167  }
168 
169  nr_double_t * u = new nr_double_t[n+1];
170  nr_double_t * z = b; // reuse storage
171  if (boundary == SPLINE_BC_NATURAL) {
172  u[0] = 0;
173  z[0] = 0;
174  } else if (boundary == SPLINE_BC_CLAMPED) {
175  u[0] = h[0] / (2 * h[0]);
176  z[0] = b[0] / (2 * h[0]);
177  }
178 
179  for (i = 1; i < n; i++) {
180  nr_double_t p = 2 * (h[i] + h[i-1]) - h[i-1] * u[i-1]; // pivot
181  u[i] = h[i] / p;
182  z[i] = (b[i] - z[i-1] * h[i-1]) / p;
183  }
184  if (boundary == SPLINE_BC_NATURAL) {
185  z[n] = 0;
186  } else if (boundary == SPLINE_BC_CLAMPED) {
187  nr_double_t p = h[n-1] * (2 - u[n-1]);
188  z[n] = (b[n] - z[n-1] * h[n-1]) / p;
189  }
190 
191  // back substitution
192  f1 = u; // reuse storage
193  f2 = z;
194  f3 = h;
195  f2[n] = z[n];
196  f3[n] = 0;
197  for (i = n - 1; i >= 0; i--) {
198  f2[i] = z[i] - u[i] * f2[i+1];
199  f1[i] = (f0[i+1] - f0[i]) / h[i] - h[i] * (f2[i+1] + 2 * f2[i]) / 3;
200  f3[i] = (f2[i+1] - f2[i]) / (3 * h[i]);
201  }
202 
203  // set up last slot for extrapolation above
204  if (boundary == SPLINE_BC_NATURAL) {
205  f1[n] = f1[n-1] + (x[n] - x[n-1]) * f2[n-1];
206  } else if (boundary == SPLINE_BC_CLAMPED) {
207  f1[n] = dn;
208  }
209  f2[n] = 0;
210  f3[n] = 0;
211  }
212 
213  // second kind of cubic splines
214  else if (boundary == SPLINE_BC_PERIODIC) {
215  // non-trigdiagonal equations - periodic boundary condition
216  nr_double_t * z = new nr_double_t[n+1];
217  if (n == 2) {
218  nr_double_t B = h[0] + h[1];
219  nr_double_t A = 2 * B;
220  nr_double_t b[2], det;
221  b[0] = 3 * ((f0[2] - f0[1]) / h[1] - (f0[1] - f0[0]) / h[0]);
222  b[1] = 3 * ((f0[1] - f0[2]) / h[0] - (f0[2] - f0[1]) / h[1]);
223  det = 3 * B * B;
224  z[1] = ( A * b[0] - B * b[1]) / det;
225  z[2] = (-B * b[0] + A * b[1]) / det;
226  z[0] = z[2];
227  }
228  else {
230  tvector<nr_double_t> o (n);
233  b.setData (&z[1], n);
234  for (i = 0; i < n - 1; i++) {
235  o(i) = h[i+1];
236  d(i) = 2 * (h[i+1] + h[i]);
237  b(i) = 3 * ((f0[i+2] - f0[i+1]) / h[i+1] - (f0[i+1] - f0[i]) / h[i]);
238  }
239  o(i) = h[0];
240  d(i) = 2 * (h[0] + h[i]);
241  b(i) = 3 * ((f0[1] - f0[i+1]) / h[0] - (f0[i+1] - f0[i]) / h[i]);
242  sys.setDiagonal (&d);
243  sys.setOffDiagonal (&o);
244  sys.setRHS (&b);
246  sys.solve ();
247  z[0] = z[n];
248  }
249 
250  f1 = new nr_double_t[n+1];
251  f2 = z; // reuse storage
252  f3 = h;
253  for (i = n - 1; i >= 0; i--) {
254  f1[i] = (f0[i+1] - f0[i]) / h[i] - h[i] * (z[i+1] + 2 * z[i]) / 3;
255  f3[i] = (z[i+1] - z[i]) / (3 * h[i]);
256  }
257  f1[n] = f1[0];
258  f2[n] = f2[0];
259  f3[n] = f3[0];
260  }
261 }
262 
263 // Returns pointer to the first value greater than the given one.
264 nr_double_t * spline::upper_bound (nr_double_t * first, nr_double_t * last,
265  nr_double_t value) {
266  int half, len = last - first;
267  nr_double_t * centre;
268 
269  while (len > 0) {
270  half = len >> 1;
271  centre = first;
272  centre += half;
273  if (value < *centre)
274  len = half;
275  else {
276  first = centre;
277  ++first;
278  len = len - half - 1;
279  }
280  }
281  return first;
282 }
283 
284 // Evaluates the spline at the given position.
285 poly spline::evaluate (nr_double_t t) {
286 
287 #ifndef PERIOD_DISABLED
288  if (boundary == SPLINE_BC_PERIODIC) {
289  // extrapolation easy: periodically
290  nr_double_t period = x[n] - x[0];
291  while (t > x[n]) t -= period;
292  while (t < x[0]) t += period;
293  }
294 #endif /* PERIOD_DISABLED */
295 
296  nr_double_t * here = upper_bound (x, x+n+1, t);
297  nr_double_t y0, y1, y2;
298  if (here == x) {
299  nr_double_t dx = t - x[0];
300  y0 = f0[0] + dx * f1[0];
301  y1 = f1[0];
302  return poly (t, y0, y1);
303  }
304  else {
305  int i = here - x - 1;
306  nr_double_t dx = t - x[i];
307  // value
308  y0 = f0[i] + dx * (f1[i] + dx * (f2[i] + dx * f3[i]));
309  // first derivative
310  y1 = f1[i] + dx * (2 * f2[i] + 3 * dx * f3[i]);
311  // second derivative
312  y2 = 2 * f2[i] + 6 * dx * f3[i];
313  return poly (t, y0, y1, y2);
314  }
315 }
316 
317 // Destructor deletes an instance of the spline class.
319  if (x) delete[] x;
320  if (f0) delete[] f0;
321  if (f1) delete[] f1;
322  if (f2) delete[] f2;
323  if (f3) delete[] f3;
324 }
325 
326 } // namespace qucs
nr_double_t d0
Definition: spline.h:70
poly evaluate(nr_double_t)
Definition: spline.cpp:285
void setDiagonal(tvector< nr_type_t > *)
Definition: tridiag.cpp:83
void setData(nr_type_t *, int)
Definition: tvector.cpp:185
matrix real(matrix a)
Real part matrix.
Definition: matrix.cpp:568
size
Definition: parse_vcd.y:203
int getSize(void) const
Definition: vector.cpp:192
void vectors(qucs::vector, qucs::vector)
Definition: spline.cpp:85
void construct(void)
Definition: spline.cpp:135
t
Definition: parse_vcd.y:290
h
Definition: parse_vcd.y:214
i
Definition: parse_mdl.y:516
nr_double_t * f0
Definition: spline.h:66
nr_double_t * f1
Definition: spline.h:67
void setType(int t)
Definition: tridiag.h:55
void solve(void)
Definition: tridiag.cpp:117
nr_complex_t det(matrix a)
Compute determinant of the given matrix.
Definition: matrix.cpp:762
#define t_
Definition: spline.cpp:81
value
Definition: parse_vcd.y:315
nr_double_t * upper_bound(nr_double_t *, nr_double_t *, nr_double_t)
Definition: spline.cpp:264
void realloc(int)
Definition: spline.cpp:121
#define y_
Definition: spline.cpp:82
void setOffDiagonal(tvector< nr_type_t > *)
Definition: tridiag.cpp:89
#define B(con)
Definition: evaluate.cpp:70
nr_double_t * f3
Definition: spline.h:69
#define A(a)
Definition: eqndefined.cpp:72
y
Definition: parse_mdl.y:499
int boundary
Definition: spline.h:72
nr_double_t * f2
Definition: spline.h:68
#define LOG_ERROR
Definition: logging.h:28
void logprint(int level, const char *format,...)
Definition: logging.c:37
void setRHS(tvector< nr_type_t > *)
Definition: tridiag.cpp:107
int getSize(void)
Definition: tvector.h:82
nr_double_t * x
Definition: spline.h:65
nr_double_t dn
Definition: spline.h:70