GNU Octave  4.2.1 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-2017 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 #if defined (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;
179  double y = z1 * (ab + z1);
180  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 -= z / (x - root[j-1]);
222  }
223
224  z /= zc;
225  x -= z;
226
227  // Famous last words: 100 iterations should be more than
228  // enough in all cases.
229
230  if (++k > 100 || octave::math::isnan (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 += 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 += 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) ("CollocWt: fatal error '%s'", msg);
379 }
380
381 CollocWt&
383 {
384  if (val >= rb)
385  error ("CollocWt: left bound greater than right bound");
386
387  lb = val;
388  initialized = 0;
389  return *this;
390 }
391
392 CollocWt&
394 {
395  if (val <= lb)
396  error ("CollocWt: right bound less than left bound");
397
398  rb = val;
399  initialized = 0;
400  return *this;
401 }
402
403 void
405 {
406  // Check for possible errors.
407
408  double wid = rb - lb;
409  if (wid <= 0.0)
410  {
411  error ("CollocWt: width less than or equal to zero");
412  return;
413  }
414
416
417  if (nt < 0)
418  error ("CollocWt: total number of collocation points less than zero");
419  else if (nt == 0)
420  return;
421
422  Array<double> dif1 (dim_vector (nt, 1));
423  double *pdif1 = dif1.fortran_vec ();
424
425  Array<double> dif2 (dim_vector (nt, 1));
426  double *pdif2 = dif2.fortran_vec ();
427
428  Array<double> dif3 (dim_vector (nt, 1));
429  double *pdif3 = dif3.fortran_vec ();
430
431  Array<double> vect (dim_vector (nt, 1));
432  double *pvect = vect.fortran_vec ();
433
434  r.resize (nt, 1);
435  q.resize (nt, 1);
436  A.resize (nt, nt);
437  B.resize (nt, nt);
438
439  double *pr = r.fortran_vec ();
440
441  // Compute roots.
442
443  if (! jcobi (n, inc_left, inc_right, Alpha, Beta, pdif1, pdif2, pdif3, pr))
444  error ("jcobi: newton iteration failed");
445
447
448  // First derivative weights.
449
450  id = 1;
451  for (octave_idx_type i = 0; i < nt; i++)
452  {
453  dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect);
454
455  for (octave_idx_type j = 0; j < nt; j++)
456  A(i,j) = vect(j);
457  }
458
459  // Second derivative weights.
460
461  id = 2;
462  for (octave_idx_type i = 0; i < nt; i++)
463  {
464  dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect);
465
466  for (octave_idx_type j = 0; j < nt; j++)
467  B(i,j) = vect(j);
468  }
469
471
472  id = 3;
473  double *pq = q.fortran_vec ();
474  dfopr (n, inc_left, inc_right, id, id, pdif1, pdif2, pdif3, pr, pq);
475
476  initialized = 1;
477 }
478
479 std::ostream&
480 operator << (std::ostream& os, const CollocWt& a)
481 {
482  if (a.left_included ())
483  os << "left boundary is included\n";
484  else
485  os << "left boundary is not included\n";
486
487  if (a.right_included ())
488  os << "right boundary is included\n";
489  else
490  os << "right boundary is not included\n";
491
492  os << "\n";
493
494  os << a.Alpha << " " << a.Beta << "\n\n"
495  << a.r << "\n\n"
496  << a.q << "\n\n"
497  << a.A << "\n"
498  << a.B << "\n";
499
500  return os;
501 }
uint32_t id
Definition: graphics.cc:11587
double Alpha
Definition: CollocWt.h:176
CollocWt & set_right(double val)
Definition: CollocWt.cc:393
double Beta
Definition: CollocWt.h:177
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:145
CollocWt & set_left(double val)
Definition: CollocWt.cc:382
octave_idx_type right_included(void) const
Definition: CollocWt.h:145
identity matrix If supplied two scalar respectively For allows like xample val
Definition: data.cc:5068
bool isnan(double x)
Definition: lo-mappers.cc:347
for large enough k
Definition: lu.cc:606
std::ostream & operator<<(std::ostream &os, const CollocWt &a)
Definition: CollocWt.cc:480
void init(void)
Definition: CollocWt.cc:404
bool initialized
Definition: CollocWt.h:185
static void dif(octave_idx_type nt, double *root, double *dif1, double *dif2, double *dif3)
Definition: CollocWt.cc:56
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
done
Definition: syscalls.cc:248
double rb
Definition: CollocWt.h:174
ColumnVector r
Definition: CollocWt.h:179
Matrix B
Definition: CollocWt.h:183
void error(const char *msg)
Definition: CollocWt.cc:376
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
ColumnVector q
Definition: CollocWt.h:180
static void dfopr(octave_idx_type n, octave_idx_type n0, octave_idx_type n1, octave_idx_type i, octave_idx_type id, double *dif1, double *dif2, double *dif3, double *root, double *vect)
Definition: CollocWt.cc:306
double lb
Definition: CollocWt.h:173
the element is set to zero In other the statement xample y
Definition: data.cc:5342
octave_idx_type n
Definition: CollocWt.h:168
octave_idx_type left_included(void) const
Definition: CollocWt.h:144
const T * fortran_vec(void) const
Definition: Array.h:584
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
Matrix A
Definition: CollocWt.h:182
octave_idx_type inc_right
Definition: CollocWt.h:171
octave_idx_type inc_left
Definition: CollocWt.h:170
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x
static bool jcobi(octave_idx_type n, octave_idx_type n0, octave_idx_type n1, double alpha, double beta, double *dif1, double *dif2, double *dif3, double *root)
Definition: CollocWt.cc:144
void resize(octave_idx_type n, const double &rfv=0)
Definition: dColVector.h:108