GNU Octave  3.8.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
CollocWt.cc
Go to the documentation of this file.
1 /*
2
3 Copyright (C) 1993-2013 John W. Eaton
4
5 This file is part of Octave.
6
7 Octave is free software; you can redistribute it and/or modify it
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
20
21 */
22
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26
27 #include <iostream>
28
29 #include <cfloat>
30
31 #include "CollocWt.h"
32 #include "f77-fcn.h"
33 #include "lo-error.h"
34
35 // The following routines jcobi, dif, and dfopr are based on the code
36 // found in Villadsen, J. and M. L. Michelsen, Solution of Differential
37 // Equation Models by Polynomial Approximation, Prentice-Hall (1978)
38 // pages 418-420.
39 //
40 // Translated to C++ by jwe.
41
42 // Compute the first three derivatives of the node polynomial.
43 //
44 // n0 (alpha,beta) n1
45 // p (x) = (x) * p (x) * (1 - x)
46 // nt n
47 //
48 // at the interpolation points. Each of the parameters n0 and n1
49 // may be given the value 0 or 1. The total number of points
50 // nt = n + n0 + n1
51 //
52 // The values of root must be known before a call to dif is possible.
53 // They may be computed using jcobi.
54
55 static void
56 dif (octave_idx_type nt, double *root, double *dif1, double *dif2,
57  double *dif3)
58 {
59  // Evaluate derivatives of node polynomial using recursion formulas.
60
61  for (octave_idx_type i = 0; i < nt; i++)
62  {
63  double x = root[i];
64
65  dif1[i] = 1.0;
66  dif2[i] = 0.0;
67  dif3[i] = 0.0;
68
69  for (octave_idx_type j = 0; j < nt; j++)
70  {
71  if (j != i)
72  {
73  double y = x - root[j];
74
75  dif3[i] = y * dif3[i] + 3.0 * dif2[i];
76  dif2[i] = y * dif2[i] + 2.0 * dif1[i];
77  dif1[i] = y * dif1[i];
78  }
79  }
80  }
81 }
82
83 // Compute the zeros of the Jacobi polynomial.
84 //
85 // (alpha,beta)
86 // p (x)
87 // n
88 //
89 // Use dif to compute the derivatives of the node
90 // polynomial
91 //
92 // n0 (alpha,beta) n1
93 // p (x) = (x) * p (x) * (1 - x)
94 // nt n
95 //
96 // at the interpolation points.
97 //
98 // See Villadsen and Michelsen, pages 131-132 and 418.
99 //
100 // Input parameters:
101 //
102 // nd : the dimension of the vectors dif1, dif2, dif3, and root
103 //
104 // n : the degree of the jacobi polynomial, (i.e. the number
105 // of interior interpolation points)
106 //
107 // n0 : determines whether x = 0 is included as an
108 // interpolation point
109 //
110 // n0 = 0 ==> x = 0 is not included
111 // n0 = 1 ==> x = 0 is included
112 //
113 // n1 : determines whether x = 1 is included as an
114 // interpolation point
115 //
116 // n1 = 0 ==> x = 1 is not included
117 // n1 = 1 ==> x = 1 is included
118 //
119 // alpha : the value of alpha in the description of the jacobi
120 // polynomial
121 //
122 // beta : the value of beta in the description of the jacobi
123 // polynomial
124 //
125 // For a more complete explanation of alpha an beta, see Villadsen
126 // and Michelsen, pages 57 to 59.
127 //
128 // Output parameters:
129 //
130 // root : one dimensional vector containing on exit the
131 // n + n0 + n1 zeros of the node polynomial used in the
132 // interpolation routine
133 //
134 // dif1 : one dimensional vector containing the first derivative
135 // of the node polynomial at the zeros
136 //
137 // dif2 : one dimensional vector containing the second derivative
138 // of the node polynomial at the zeros
139 //
140 // dif3 : one dimensional vector containing the third derivative
141 // of the node polynomial at the zeros
142
143 static bool
145  double alpha, double beta, double *dif1, double *dif2,
146  double *dif3, double *root)
147 {
148  assert (n0 == 0 || n0 == 1);
149  assert (n1 == 0 || n1 == 1);
150
151  octave_idx_type nt = n + n0 + n1;
152
153  assert (nt > 1);
154
155 // -- first evaluation of coefficients in recursion formulas.
156 // -- recursion coefficients are stored in dif1 and dif2.
157
158  double ab = alpha + beta;
159  double ad = beta - alpha;
160  double ap = beta * alpha;
161
162  dif1[0] = (ad / (ab + 2.0) + 1.0) / 2.0;
163  dif2[0] = 0.0;
164
165  if (n >= 2)
166  {
167  for (octave_idx_type i = 1; i < n; i++)
168  {
169  double z1 = i;
170  double z = ab + 2 * z1;
171
172  dif1[i] = (ab * ad / z / (z + 2.0) + 1.0) / 2.0;
173
174  if (i == 1)
175  dif2[i] = (ab + ap + z1) / z / z / (z + 1.0);
176  else
177  {
178  z = z * z;
179  double y = z1 * (ab + z1);
180  y = y * (ap + y);
181  dif2[i] = y / z / (z - 1.0);
182  }
183  }
184  }
185
186  // Root determination by Newton method with suppression of previously
187  // determined roots.
188
189  double x = 0.0;
190
191  for (octave_idx_type i = 0; i < n; i++)
192  {
193  bool done = false;
194
195  int k = 0;
196
197  while (! done)
198  {
199  double xd = 0.0;
200  double xn = 1.0;
201  double xd1 = 0.0;
202  double xn1 = 0.0;
203
204  for (octave_idx_type j = 0; j < n; j++)
205  {
206  double xp = (dif1[j] - x) * xn - dif2[j] * xd;
207  double xp1 = (dif1[j] - x) * xn1 - dif2[j] * xd1 - xn;
208
209  xd = xn;
210  xd1 = xn1;
211  xn = xp;
212  xn1 = xp1;
213  }
214
215  double zc = 1.0;
216  double z = xn / xn1;
217
218  if (i != 0)
219  {
220  for (octave_idx_type j = 1; j <= i; j++)
221  zc = zc - z / (x - root[j-1]);
222  }
223
224  z = z / zc;
225  x = x - z;
226
227  // Famous last words: 100 iterations should be more than
228  // enough in all cases.
229
230  if (++k > 100 || xisnan (z))
231  return false;
232
233  if (std::abs (z) <= 100 * std::numeric_limits<double>::epsilon ())
234  done = true;
235  }
236
237  root[i] = x;
238  x = x + sqrt (std::numeric_limits<double>::epsilon ());
239  }
240
241  // Add interpolation points at x = 0 and/or x = 1.
242
243  if (n0 != 0)
244  {
245  for (octave_idx_type i = n; i > 0; i--)
246  root[i] = root[i-1];
247
248  root[0] = 0.0;
249  }
250
251  if (n1 != 0)
252  root[nt-1] = 1.0;
253
254  dif (nt, root, dif1, dif2, dif3);
255
256  return true;
257 }
258
259 // Compute derivative weights for orthogonal collocation.
260 //
261 // See Villadsen and Michelsen, pages 133-134, 419.
262 //
263 // Input parameters:
264 //
265 // nd : the dimension of the vectors dif1, dif2, dif3, and root
266 //
267 // n : the degree of the jacobi polynomial, (i.e. the number
268 // of interior interpolation points)
269 //
270 // n0 : determines whether x = 0 is included as an
271 // interpolation point
272 //
273 // n0 = 0 ==> x = 0 is not included
274 // n0 = 1 ==> x = 0 is included
275 //
276 // n1 : determines whether x = 1 is included as an
277 // interpolation point
278 //
279 // n1 = 0 ==> x = 1 is not included
280 // n1 = 1 ==> x = 1 is included
281 //
282 // i : the index of the node for which the weights are to be
283 // calculated
284 //
285 // id : indicator
286 //
287 // id = 1 ==> first derivative weights are computed
288 // id = 2 ==> second derivative weights are computed
289 // id = 3 ==> gaussian weights are computed (in this
290 // case, the value of i is irrelevant)
291 //
292 // Output parameters:
293 //
294 // dif1 : one dimensional vector containing the first derivative
295 // of the node polynomial at the zeros
296 //
297 // dif2 : one dimensional vector containing the second derivative
298 // of the node polynomial at the zeros
299 //
300 // dif3 : one dimensional vector containing the third derivative
301 // of the node polynomial at the zeros
302 //
303 // vect : one dimensional vector of computed weights
304
305 static void
307  octave_idx_type i, octave_idx_type id, double *dif1,
308  double *dif2, double *dif3, double *root, double *vect)
309 {
310  assert (n0 == 0 || n0 == 1);
311  assert (n1 == 0 || n1 == 1);
312
313  octave_idx_type nt = n + n0 + n1;
314
315  assert (nt > 1);
316
317  assert (id == 1 || id == 2 || id == 3);
318
319  if (id != 3)
320  assert (i >= 0 && i < nt);
321
322  // Evaluate discretization matrices and Gaussian quadrature weights.
323  // Quadrature weights are normalized to sum to one.
324
325  if (id != 3)
326  {
327  for (octave_idx_type j = 0; j < nt; j++)
328  {
329  if (j == i)
330  {
331  if (id == 1)
332  vect[i] = dif2[i] / dif1[i] / 2.0;
333  else
334  vect[i] = dif3[i] / dif1[i] / 3.0;
335  }
336  else
337  {
338  double y = root[i] - root[j];
339
340  vect[j] = dif1[i] / dif1[j] / y;
341
342  if (id == 2)
343  vect[j] = vect[j] * (dif2[i] / dif1[i] - 2.0 / y);
344  }
345  }
346  }
347  else
348  {
349  double y = 0.0;
350
351  for (octave_idx_type j = 0; j < nt; j++)
352  {
353  double x = root[j];
354
355  double ax = x * (1.0 - x);
356
357  if (n0 == 0)
358  ax = ax / x / x;
359
360  if (n1 == 0)
361  ax = ax / (1.0 - x) / (1.0 - x);
362
363  vect[j] = ax / (dif1[j] * dif1[j]);
364
365  y = y + vect[j];
366  }
367
368  for (octave_idx_type j = 0; j < nt; j++)
369  vect[j] = vect[j] / y;
370  }
371 }
372
373 // Error handling.
374
375 void
376 CollocWt::error (const char* msg)
377 {
378  (*current_liboctave_error_handler) ("fatal CollocWt error: %s", msg);
379 }
380
381 CollocWt&
382 CollocWt::set_left (double val)
383 {
384  if (val >= rb)
385  {
386  error ("left bound greater than right bound");
387  return *this;
388  }
389
390  lb = val;
391  initialized = 0;
392  return *this;
393 }
394
395 CollocWt&
397 {
398  if (val <= lb)
399  {
400  error ("right bound less than left bound");
401  return *this;
402  }
403
404  rb = val;
405  initialized = 0;
406  return *this;
407 }
408
409 void
411 {
412  // Check for possible errors.
413
414  double wid = rb - lb;
415  if (wid <= 0.0)
416  {
417  error ("width less than or equal to zero");
418  return;
419  }
420
422
423  if (nt < 0)
424  {
425  error ("total number of collocation points less than zero");
426  return;
427  }
428  else if (nt == 0)
429  return;
430
431  Array<double> dif1 (dim_vector (nt, 1));
432  double *pdif1 = dif1.fortran_vec ();
433
434  Array<double> dif2 (dim_vector (nt, 1));
435  double *pdif2 = dif2.fortran_vec ();
436
437  Array<double> dif3 (dim_vector (nt, 1));
438  double *pdif3 = dif3.fortran_vec ();
439
440  Array<double> vect (dim_vector (nt, 1));
441  double *pvect = vect.fortran_vec ();
442
443  r.resize (nt, 1);
444  q.resize (nt, 1);
445  A.resize (nt, nt);
446  B.resize (nt, nt);
447
448  double *pr = r.fortran_vec ();
449
450  // Compute roots.
451
452  if (! jcobi (n, inc_left, inc_right, Alpha, Beta, pdif1, pdif2, pdif3, pr))
453  {
454  error ("jcobi: newton iteration failed");
455  return;
456  }
457
458  octave_idx_type id;
459
460  // First derivative weights.
461
462  id = 1;
463  for (octave_idx_type i = 0; i < nt; i++)
464  {
465  dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect);
466
467  for (octave_idx_type j = 0; j < nt; j++)
468  A(i,j) = vect(j);
469  }
470
471  // Second derivative weights.
472
473  id = 2;
474  for (octave_idx_type i = 0; i < nt; i++)
475  {
476  dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect);
477
478  for (octave_idx_type j = 0; j < nt; j++)
479  B(i,j) = vect(j);
480  }
481
483
484  id = 3;
485  double *pq = q.fortran_vec ();
486  dfopr (n, inc_left, inc_right, id, id, pdif1, pdif2, pdif3, pr, pq);
487
488  initialized = 1;
489 }
490
491 std::ostream&
492 operator << (std::ostream& os, const CollocWt& a)
493 {
494  if (a.left_included ())
495  os << "left boundary is included\n";
496  else
497  os << "left boundary is not included\n";
498
499  if (a.right_included ())
500  os << "right boundary is included\n";
501  else
502  os << "right boundary is not included\n";
503
504  os << "\n";
505
506  os << a.Alpha << " " << a.Beta << "\n\n"
507  << a.r << "\n\n"
508  << a.q << "\n\n"
509  << a.A << "\n"
510  << a.B << "\n";
511
512  return os;
513 }