Mercurial > octave
annotate liboctave/numeric/CollocWt.cc @ 23443:3f1bf237908b
maint: Eliminate <cfloat.h> header from liboctave files.
* build-aux/mk-opts.pl: Remove "#include <cfloat>" from script which creates *-opts.h
* Array-util.cc, Array-util.h, Array.cc, Array.h, CColVector.h,
CDiagMatrix.h, CMatrix.cc, CMatrix.h, CNDArray.cc, CNDArray.h, CRowVector.h,
CSparse.cc, CSparse.h, MSparse.h, Matrix.h, Range.cc, Sparse.h, boolMatrix.h,
boolNDArray.h, boolSparse.h, chMatrix.h, chNDArray.h, dColVector.h,
dDiagMatrix.h, dMatrix.cc, dMatrix.h, dNDArray.cc, dNDArray.h, dRowVector.h,
dSparse.cc, dSparse.h, dim-vector.cc, dim-vector.h, fCColVector.h,
fCDiagMatrix.h, fCMatrix.cc, fCMatrix.h, fCNDArray.cc, fCNDArray.h,
fCRowVector.h, fColVector.h, fDiagMatrix.h, fMatrix.cc, fMatrix.h, fNDArray.cc,
fNDArray.h, fRowVector.h, int16NDArray.h, int32NDArray.h, int64NDArray.h,
int8NDArray.h, intNDArray.h, uint16NDArray.h, uint32NDArray.h, uint64NDArray.h,
uint8NDArray.h, CollocWt.cc, DASPK.cc, DASPK.h, DASRT.cc, DASRT.h, DASSL.cc,
DASSL.h, LSODE.cc, LSODE.h, Quad.h, eigs-base.cc, eigs-base.h, lo-mappers.cc,
lo-mappers.h, oct-norm.cc, lo-ieee.cc, lo-utils.cc, lo-utils.h, oct-refcount.h:
Remove "#include <cfloat>". Add "#include <limits>" if necessary. Remove
redundant #includes. Alphabetize #includes.
author | Rik <rik@octave.org> |
---|---|
date | Wed, 26 Apr 2017 09:57:28 -0700 |
parents | 092078913d54 |
children | c763214a8260 |
rev | line source |
---|---|
3 | 1 /* |
2 | |
23219
3ac9f9ecfae5
maint: Update copyright dates.
John W. Eaton <jwe@octave.org>
parents:
23083
diff
changeset
|
3 Copyright (C) 1993-2017 John W. Eaton |
3 | 4 |
5 This file is part of Octave. | |
6 | |
7 Octave is free software; you can redistribute it and/or modify it | |
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
8 under the terms of the GNU General Public License as published by |
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
9 the Free Software Foundation; either version 3 of the License, or |
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
10 (at your option) any later version. |
3 | 11 |
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
12 Octave is distributed in the hope that it will be useful, but |
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
13 WITHOUT ANY WARRANTY; without even the implied warranty of |
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
15 GNU General Public License for more details. |
3 | 16 |
17 You should have received a copy of the GNU General Public License | |
7016 | 18 along with Octave; see the file COPYING. If not, see |
19 <http://www.gnu.org/licenses/>. | |
3 | 20 |
21 */ | |
22 | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21568
diff
changeset
|
23 #if defined (HAVE_CONFIG_H) |
21301
40de9f8f23a6
Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents:
21202
diff
changeset
|
24 # include "config.h" |
3 | 25 #endif |
26 | |
3503 | 27 #include <iostream> |
23443
3f1bf237908b
maint: Eliminate <cfloat.h> header from liboctave files.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
28 #include <limits> |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
29 |
3 | 30 #include "CollocWt.h" |
1847 | 31 #include "f77-fcn.h" |
227 | 32 #include "lo-error.h" |
3 | 33 |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
34 // The following routines jcobi, dif, and dfopr are based on the code |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
35 // found in Villadsen, J. and M. L. Michelsen, Solution of Differential |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
36 // Equation Models by Polynomial Approximation, Prentice-Hall (1978) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
37 // pages 418-420. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
38 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
39 // Translated to C++ by jwe. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
40 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
41 // Compute the first three derivatives of the node polynomial. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
42 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
43 // n0 (alpha,beta) n1 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
44 // p (x) = (x) * p (x) * (1 - x) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
45 // nt n |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
46 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
47 // at the interpolation points. Each of the parameters n0 and n1 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
48 // may be given the value 0 or 1. The total number of points |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
49 // nt = n + n0 + n1 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
50 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
51 // The values of root must be known before a call to dif is possible. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
52 // They may be computed using jcobi. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
53 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
54 static void |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
55 dif (octave_idx_type nt, double *root, double *dif1, double *dif2, |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
56 double *dif3) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
57 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
58 // Evaluate derivatives of node polynomial using recursion formulas. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
59 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
60 for (octave_idx_type i = 0; i < nt; i++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
61 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
62 double x = root[i]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
63 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
64 dif1[i] = 1.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
65 dif2[i] = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
66 dif3[i] = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
67 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
68 for (octave_idx_type j = 0; j < nt; j++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
69 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
70 if (j != i) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
71 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
72 double y = x - root[j]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
73 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
74 dif3[i] = y * dif3[i] + 3.0 * dif2[i]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
75 dif2[i] = y * dif2[i] + 2.0 * dif1[i]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
76 dif1[i] = y * dif1[i]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
77 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
78 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
79 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
80 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
81 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
82 // Compute the zeros of the Jacobi polynomial. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
83 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
84 // (alpha,beta) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
85 // p (x) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
86 // n |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
87 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
88 // Use dif to compute the derivatives of the node |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
89 // polynomial |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
90 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
91 // n0 (alpha,beta) n1 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
92 // p (x) = (x) * p (x) * (1 - x) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
93 // nt n |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
94 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
95 // at the interpolation points. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
96 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
97 // See Villadsen and Michelsen, pages 131-132 and 418. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
98 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
99 // Input parameters: |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
100 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
101 // nd : the dimension of the vectors dif1, dif2, dif3, and root |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
102 // |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
103 // n : the degree of the jacobi polynomial, (i.e., the number |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
104 // of interior interpolation points) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
105 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
106 // n0 : determines whether x = 0 is included as an |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
107 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
108 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
109 // n0 = 0 ==> x = 0 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
110 // n0 = 1 ==> x = 0 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
111 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
112 // n1 : determines whether x = 1 is included as an |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
113 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
114 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
115 // n1 = 0 ==> x = 1 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
116 // n1 = 1 ==> x = 1 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
117 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
118 // alpha : the value of alpha in the description of the jacobi |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
119 // polynomial |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
120 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
121 // beta : the value of beta in the description of the jacobi |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
122 // polynomial |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
123 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
124 // For a more complete explanation of alpha an beta, see Villadsen |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
125 // and Michelsen, pages 57 to 59. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
126 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
127 // Output parameters: |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
128 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
129 // root : one dimensional vector containing on exit the |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
130 // n + n0 + n1 zeros of the node polynomial used in the |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
131 // interpolation routine |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
132 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
133 // dif1 : one dimensional vector containing the first derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
134 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
135 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
136 // dif2 : one dimensional vector containing the second derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
137 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
138 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
139 // dif3 : one dimensional vector containing the third derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
140 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
141 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
142 static bool |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
143 jcobi (octave_idx_type n, octave_idx_type n0, octave_idx_type n1, |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
144 double alpha, double beta, double *dif1, double *dif2, |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
145 double *dif3, double *root) |
3 | 146 { |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
147 assert (n0 == 0 || n0 == 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
148 assert (n1 == 0 || n1 == 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
149 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
150 octave_idx_type nt = n + n0 + n1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
151 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
152 assert (nt > 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
153 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
154 // -- first evaluation of coefficients in recursion formulas. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
155 // -- recursion coefficients are stored in dif1 and dif2. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
156 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
157 double ab = alpha + beta; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
158 double ad = beta - alpha; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
159 double ap = beta * alpha; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
160 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
161 dif1[0] = (ad / (ab + 2.0) + 1.0) / 2.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
162 dif2[0] = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
163 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
164 if (n >= 2) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
165 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
166 for (octave_idx_type i = 1; i < n; i++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
167 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
168 double z1 = i; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
169 double z = ab + 2 * z1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
170 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
171 dif1[i] = (ab * ad / z / (z + 2.0) + 1.0) / 2.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
172 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
173 if (i == 1) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
174 dif2[i] = (ab + ap + z1) / z / z / (z + 1.0); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
175 else |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
176 { |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
177 z *= z; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
178 double y = z1 * (ab + z1); |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
179 y *= (ap + y); |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
180 dif2[i] = y / z / (z - 1.0); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
181 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
182 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
183 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
184 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
185 // Root determination by Newton method with suppression of previously |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
186 // determined roots. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
187 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
188 double x = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
189 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
190 for (octave_idx_type i = 0; i < n; i++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
191 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
192 bool done = false; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
193 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
194 int k = 0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
195 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
196 while (! done) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
197 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
198 double xd = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
199 double xn = 1.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
200 double xd1 = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
201 double xn1 = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
202 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
203 for (octave_idx_type j = 0; j < n; j++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
204 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
205 double xp = (dif1[j] - x) * xn - dif2[j] * xd; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
206 double xp1 = (dif1[j] - x) * xn1 - dif2[j] * xd1 - xn; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
207 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
208 xd = xn; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
209 xd1 = xn1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
210 xn = xp; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
211 xn1 = xp1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
212 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
213 |
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
214 double zc = 1.0; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
215 double z = xn / xn1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
216 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
217 if (i != 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
218 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
219 for (octave_idx_type j = 1; j <= i; j++) |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
220 zc -= z / (x - root[j-1]); |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
221 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
222 |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
223 z /= zc; |
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
224 x -= z; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
225 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
226 // Famous last words: 100 iterations should be more than |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
227 // enough in all cases. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
228 |
21782
2aef506f3fec
use namespace for lo-mappers.h functions
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
229 if (++k > 100 || octave::math::isnan (z)) |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
230 return false; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
231 |
15220
61822c866ba1
use std::numeric_limits<T>::epsilon in C++ code
John W. Eaton <jwe@octave.org>
parents:
15018
diff
changeset
|
232 if (std::abs (z) <= 100 * std::numeric_limits<double>::epsilon ()) |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
233 done = true; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
234 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
235 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
236 root[i] = x; |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
237 x += sqrt (std::numeric_limits<double>::epsilon ()); |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
238 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
239 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
240 // Add interpolation points at x = 0 and/or x = 1. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
241 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
242 if (n0 != 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
243 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
244 for (octave_idx_type i = n; i > 0; i--) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
245 root[i] = root[i-1]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
246 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
247 root[0] = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
248 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
249 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
250 if (n1 != 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
251 root[nt-1] = 1.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
252 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
253 dif (nt, root, dif1, dif2, dif3); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
254 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
255 return true; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
256 } |
3 | 257 |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
258 // Compute derivative weights for orthogonal collocation. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
259 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
260 // See Villadsen and Michelsen, pages 133-134, 419. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
261 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
262 // Input parameters: |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
263 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
264 // nd : the dimension of the vectors dif1, dif2, dif3, and root |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
265 // |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
266 // n : the degree of the jacobi polynomial, (i.e., the number |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
267 // of interior interpolation points) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
268 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
269 // n0 : determines whether x = 0 is included as an |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
270 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
271 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
272 // n0 = 0 ==> x = 0 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
273 // n0 = 1 ==> x = 0 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
274 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
275 // n1 : determines whether x = 1 is included as an |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
276 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
277 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
278 // n1 = 0 ==> x = 1 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
279 // n1 = 1 ==> x = 1 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
280 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
281 // i : the index of the node for which the weights are to be |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
282 // calculated |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
283 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
284 // id : indicator |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
285 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
286 // id = 1 ==> first derivative weights are computed |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
287 // id = 2 ==> second derivative weights are computed |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
288 // id = 3 ==> gaussian weights are computed (in this |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
289 // case, the value of i is irrelevant) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
290 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
291 // Output parameters: |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
292 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
293 // dif1 : one dimensional vector containing the first derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
294 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
295 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
296 // dif2 : one dimensional vector containing the second derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
297 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
298 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
299 // dif3 : one dimensional vector containing the third derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
300 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
301 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
302 // vect : one dimensional vector of computed weights |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
303 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
304 static void |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
305 dfopr (octave_idx_type n, octave_idx_type n0, octave_idx_type n1, |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
306 octave_idx_type i, octave_idx_type id, double *dif1, |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
307 double *dif2, double *dif3, double *root, double *vect) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
308 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
309 assert (n0 == 0 || n0 == 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
310 assert (n1 == 0 || n1 == 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
311 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
312 octave_idx_type nt = n + n0 + n1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
313 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
314 assert (nt > 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
315 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
316 assert (id == 1 || id == 2 || id == 3); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
317 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
318 if (id != 3) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
319 assert (i >= 0 && i < nt); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
320 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
321 // Evaluate discretization matrices and Gaussian quadrature weights. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
322 // Quadrature weights are normalized to sum to one. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
323 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
324 if (id != 3) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
325 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
326 for (octave_idx_type j = 0; j < nt; j++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
327 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
328 if (j == i) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
329 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
330 if (id == 1) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
331 vect[i] = dif2[i] / dif1[i] / 2.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
332 else |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
333 vect[i] = dif3[i] / dif1[i] / 3.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
334 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
335 else |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
336 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
337 double y = root[i] - root[j]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
338 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
339 vect[j] = dif1[i] / dif1[j] / y; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
340 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
341 if (id == 2) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
342 vect[j] = vect[j] * (dif2[i] / dif1[i] - 2.0 / y); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
343 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
344 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
345 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
346 else |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
347 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
348 double y = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
349 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
350 for (octave_idx_type j = 0; j < nt; j++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
351 { |
21568
3d60ed163b70
maint: Eliminate bad spacing around '='.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
352 double x = root[j]; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
353 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
354 double ax = x * (1.0 - x); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
355 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
356 if (n0 == 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
357 ax = ax / x / x; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
358 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
359 if (n1 == 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
360 ax = ax / (1.0 - x) / (1.0 - x); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
361 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
362 vect[j] = ax / (dif1[j] * dif1[j]); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
363 |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
364 y += vect[j]; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
365 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
366 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
367 for (octave_idx_type j = 0; j < nt; j++) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
368 vect[j] = vect[j] / y; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
369 } |
3 | 370 } |
371 | |
372 // Error handling. | |
373 | |
374 void | |
375 CollocWt::error (const char* msg) | |
376 { | |
20726
25d676f9619c
Preface error() messages with name of function when possible.
Rik <rik@octave.org>
parents:
20230
diff
changeset
|
377 (*current_liboctave_error_handler) ("CollocWt: fatal error '%s'", msg); |
3 | 378 } |
379 | |
380 CollocWt& | |
381 CollocWt::set_left (double val) | |
382 { | |
383 if (val >= rb) | |
21055
5e00ed38a58b
maint: Replace if/error/else paradigm with just if/error.
Rik <rik@octave.org>
parents:
20726
diff
changeset
|
384 error ("CollocWt: left bound greater than right bound"); |
3 | 385 |
386 lb = val; | |
387 initialized = 0; | |
388 return *this; | |
389 } | |
390 | |
391 CollocWt& | |
392 CollocWt::set_right (double val) | |
393 { | |
394 if (val <= lb) | |
21055
5e00ed38a58b
maint: Replace if/error/else paradigm with just if/error.
Rik <rik@octave.org>
parents:
20726
diff
changeset
|
395 error ("CollocWt: right bound less than left bound"); |
3 | 396 |
397 rb = val; | |
398 initialized = 0; | |
399 return *this; | |
400 } | |
401 | |
402 void | |
403 CollocWt::init (void) | |
404 { | |
1360 | 405 // Check for possible errors. |
3 | 406 |
407 double wid = rb - lb; | |
408 if (wid <= 0.0) | |
227 | 409 { |
20726
25d676f9619c
Preface error() messages with name of function when possible.
Rik <rik@octave.org>
parents:
20230
diff
changeset
|
410 error ("CollocWt: width less than or equal to zero"); |
227 | 411 return; |
412 } | |
3 | 413 |
5275 | 414 octave_idx_type nt = n + inc_left + inc_right; |
1870 | 415 |
3 | 416 if (nt < 0) |
21055
5e00ed38a58b
maint: Replace if/error/else paradigm with just if/error.
Rik <rik@octave.org>
parents:
20726
diff
changeset
|
417 error ("CollocWt: total number of collocation points less than zero"); |
3 | 418 else if (nt == 0) |
419 return; | |
420 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
421 Array<double> dif1 (dim_vector (nt, 1)); |
1938 | 422 double *pdif1 = dif1.fortran_vec (); |
423 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
424 Array<double> dif2 (dim_vector (nt, 1)); |
1938 | 425 double *pdif2 = dif2.fortran_vec (); |
426 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
427 Array<double> dif3 (dim_vector (nt, 1)); |
1938 | 428 double *pdif3 = dif3.fortran_vec (); |
429 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
430 Array<double> vect (dim_vector (nt, 1)); |
1938 | 431 double *pvect = vect.fortran_vec (); |
3 | 432 |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
433 r.resize (nt, 1); |
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
434 q.resize (nt, 1); |
3 | 435 A.resize (nt, nt); |
436 B.resize (nt, nt); | |
437 | |
438 double *pr = r.fortran_vec (); | |
439 | |
1360 | 440 // Compute roots. |
3 | 441 |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
442 if (! jcobi (n, inc_left, inc_right, Alpha, Beta, pdif1, pdif2, pdif3, pr)) |
21056
d48fdf3a8c0c
maint: one more addition to cset 5e00ed38a58b.
Rik <rik@octave.org>
parents:
21055
diff
changeset
|
443 error ("jcobi: newton iteration failed"); |
3 | 444 |
5275 | 445 octave_idx_type id; |
3 | 446 |
1360 | 447 // First derivative weights. |
3 | 448 |
449 id = 1; | |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
450 for (octave_idx_type i = 0; i < nt; i++) |
3 | 451 { |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
452 dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect); |
3 | 453 |
5275 | 454 for (octave_idx_type j = 0; j < nt; j++) |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
455 A(i,j) = vect(j); |
3 | 456 } |
457 | |
1360 | 458 // Second derivative weights. |
3 | 459 |
460 id = 2; | |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
461 for (octave_idx_type i = 0; i < nt; i++) |
3 | 462 { |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
463 dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect); |
3 | 464 |
5275 | 465 for (octave_idx_type j = 0; j < nt; j++) |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
466 B(i,j) = vect(j); |
3 | 467 } |
468 | |
1360 | 469 // Gaussian quadrature weights. |
3 | 470 |
471 id = 3; | |
472 double *pq = q.fortran_vec (); | |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
473 dfopr (n, inc_left, inc_right, id, id, pdif1, pdif2, pdif3, pr, pq); |
3 | 474 |
475 initialized = 1; | |
476 } | |
477 | |
3504 | 478 std::ostream& |
479 operator << (std::ostream& os, const CollocWt& a) | |
3 | 480 { |
481 if (a.left_included ()) | |
482 os << "left boundary is included\n"; | |
483 else | |
484 os << "left boundary is not included\n"; | |
485 | |
486 if (a.right_included ()) | |
487 os << "right boundary is included\n"; | |
488 else | |
489 os << "right boundary is not included\n"; | |
490 | |
491 os << "\n"; | |
492 | |
493 os << a.Alpha << " " << a.Beta << "\n\n" | |
494 << a.r << "\n\n" | |
495 << a.q << "\n\n" | |
496 << a.A << "\n" | |
497 << a.B << "\n"; | |
498 | |
499 return os; | |
500 } |