Mercurial > octave
annotate liboctave/numeric/CollocWt.cc @ 27919:1891570abac8
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2020.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Mon, 06 Jan 2020 22:29:51 -0500 |
parents | b442ec6dda5c |
children | bd51beb6205e |
rev | line source |
---|---|
3 | 1 /* |
2 | |
27919
1891570abac8
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
27918
diff
changeset
|
3 Copyright (C) 1993-2020 The Octave Project Developers |
27918
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
4 |
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
5 See the file COPYRIGHT.md in the top-level directory of this distribution |
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
6 or <https://octave.org/COPYRIGHT.html/>. |
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
7 |
3 | 8 |
9 This file is part of Octave. | |
10 | |
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23807
diff
changeset
|
11 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
|
12 under the terms of the GNU General Public License as published by |
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23807
diff
changeset
|
13 the Free Software Foundation, either version 3 of the License, or |
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
14 (at your option) any later version. |
3 | 15 |
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
16 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
|
17 WITHOUT ANY WARRANTY; without even the implied warranty of |
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22402
diff
changeset
|
19 GNU General Public License for more details. |
3 | 20 |
21 You should have received a copy of the GNU General Public License | |
7016 | 22 along with Octave; see the file COPYING. If not, see |
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23807
diff
changeset
|
23 <https://www.gnu.org/licenses/>. |
3 | 24 |
25 */ | |
26 | |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21568
diff
changeset
|
27 #if defined (HAVE_CONFIG_H) |
21301
40de9f8f23a6
Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents:
21202
diff
changeset
|
28 # include "config.h" |
3 | 29 #endif |
30 | |
23475
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23449
diff
changeset
|
31 #include <cassert> |
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23449
diff
changeset
|
32 #include <cmath> |
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23449
diff
changeset
|
33 |
23443
3f1bf237908b
maint: Eliminate <cfloat.h> header from liboctave files.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
34 #include <limits> |
25438
cb1606f78f6b
prefer <istream>, <ostream>, or <iosfwd> to <iostream> where possible
John W. Eaton <jwe@octave.org>
parents:
25247
diff
changeset
|
35 #include <ostream> |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
36 |
23475
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23449
diff
changeset
|
37 #include "Array.h" |
3 | 38 #include "CollocWt.h" |
227 | 39 #include "lo-error.h" |
23475
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23449
diff
changeset
|
40 #include "lo-mappers.h" |
3 | 41 |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
42 // 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
|
43 // 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
|
44 // Equation Models by Polynomial Approximation, Prentice-Hall (1978) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
45 // pages 418-420. |
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 // Translated to C++ by jwe. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
48 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
49 // 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
|
50 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
51 // n0 (alpha,beta) n1 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
52 // p (x) = (x) * p (x) * (1 - x) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
53 // nt n |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
54 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
55 // 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
|
56 // 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
|
57 // nt = n + n0 + n1 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
58 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
59 // 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
|
60 // They may be computed using jcobi. |
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 static void |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
63 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
|
64 double *dif3) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
65 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
66 // Evaluate derivatives of node polynomial using recursion formulas. |
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 i = 0; i < nt; i++) |
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 double x = root[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 dif1[i] = 1.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
73 dif2[i] = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
74 dif3[i] = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
75 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
76 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
|
77 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
78 if (j != i) |
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 double y = x - root[j]; |
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 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
|
83 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
|
84 dif1[i] = y * dif1[i]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
85 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
86 } |
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 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
89 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
90 // Compute the zeros of the Jacobi polynomial. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
91 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
92 // (alpha,beta) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
93 // p (x) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
94 // n |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
95 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
96 // 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
|
97 // polynomial |
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 // n0 (alpha,beta) n1 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
100 // p (x) = (x) * p (x) * (1 - x) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
101 // nt n |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
102 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
103 // at the interpolation points. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
104 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
105 // 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
|
106 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
107 // Input parameters: |
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 // 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
|
110 // |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
111 // 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
|
112 // of interior interpolation points) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
113 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
114 // 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
|
115 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
116 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
117 // n0 = 0 ==> x = 0 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
118 // n0 = 1 ==> x = 0 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
119 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
120 // 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
|
121 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
122 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
123 // n1 = 0 ==> x = 1 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
124 // n1 = 1 ==> x = 1 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
125 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
126 // 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
|
127 // polynomial |
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 // 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
|
130 // polynomial |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
131 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
132 // 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
|
133 // and Michelsen, pages 57 to 59. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
134 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
135 // Output parameters: |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
136 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
137 // root : one dimensional vector containing on exit the |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
138 // 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
|
139 // interpolation routine |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
140 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
141 // dif1 : one dimensional vector containing the first derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
142 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
143 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
144 // dif2 : one dimensional vector containing the second derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
145 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
146 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
147 // dif3 : one dimensional vector containing the third derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
148 // of the node polynomial at the zeros |
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 static bool |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
151 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
|
152 double alpha, double beta, double *dif1, double *dif2, |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
153 double *dif3, double *root) |
3 | 154 { |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
155 assert (n0 == 0 || n0 == 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
156 assert (n1 == 0 || n1 == 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
157 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
158 octave_idx_type nt = n + n0 + n1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
159 |
25247
db31e068f4db
Rewrite incorrect assert statement in colloc calculation (bug #53653)
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
160 assert (nt >= 1); |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
161 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
162 // -- first evaluation of coefficients in recursion formulas. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
163 // -- recursion coefficients are stored in dif1 and dif2. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
164 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
165 double ab = alpha + beta; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
166 double ad = beta - alpha; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
167 double ap = beta * alpha; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
168 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
169 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
|
170 dif2[0] = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
171 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
172 if (n >= 2) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
173 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
174 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
|
175 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
176 double z1 = i; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
177 double z = ab + 2 * z1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
178 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
179 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
|
180 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
181 if (i == 1) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
182 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
|
183 else |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
184 { |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
185 z *= z; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
186 double y = z1 * (ab + z1); |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
187 y *= (ap + y); |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
188 dif2[i] = y / z / (z - 1.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 } |
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 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
193 // 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
|
194 // determined roots. |
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 double x = 0.0; |
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 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
|
199 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
200 bool done = false; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
201 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
202 int k = 0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
203 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
204 while (! done) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
205 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
206 double xd = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
207 double xn = 1.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
208 double xd1 = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
209 double xn1 = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
210 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
211 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
|
212 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
213 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
|
214 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
|
215 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
216 xd = xn; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
217 xd1 = xn1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
218 xn = xp; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
219 xn1 = xp1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
220 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
221 |
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
222 double zc = 1.0; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
223 double z = xn / xn1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
224 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
225 if (i != 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
226 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
227 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
|
228 zc -= z / (x - root[j-1]); |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
229 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
230 |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
231 z /= zc; |
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
232 x -= z; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
233 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
234 // 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
|
235 // enough in all cases. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
236 |
21782
2aef506f3fec
use namespace for lo-mappers.h functions
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
237 if (++k > 100 || octave::math::isnan (z)) |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
238 return false; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
239 |
15220
61822c866ba1
use std::numeric_limits<T>::epsilon in C++ code
John W. Eaton <jwe@octave.org>
parents:
15018
diff
changeset
|
240 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
|
241 done = true; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
242 } |
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 root[i] = x; |
23649
aabf6cd222ac
Use sqrt from C++ std library in liboctave.
Rik <rik@octave.org>
parents:
23475
diff
changeset
|
245 x += std::sqrt (std::numeric_limits<double>::epsilon ()); |
10603
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 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
248 // 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
|
249 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
250 if (n0 != 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
251 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
252 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
|
253 root[i] = root[i-1]; |
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 root[0] = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
256 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
257 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
258 if (n1 != 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
259 root[nt-1] = 1.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
260 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
261 dif (nt, root, dif1, dif2, dif3); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
262 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
263 return true; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
264 } |
3 | 265 |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
266 // Compute derivative weights for orthogonal collocation. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
267 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
268 // See Villadsen and Michelsen, pages 133-134, 419. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
269 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
270 // Input parameters: |
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 // 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
|
273 // |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
274 // 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
|
275 // of interior interpolation points) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
276 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
277 // 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
|
278 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
279 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
280 // n0 = 0 ==> x = 0 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
281 // n0 = 1 ==> x = 0 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
282 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
283 // 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
|
284 // interpolation point |
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 // n1 = 0 ==> x = 1 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
287 // n1 = 1 ==> x = 1 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
288 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
289 // 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
|
290 // calculated |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
291 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
292 // id : indicator |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
293 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
294 // id = 1 ==> first derivative weights are computed |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
295 // id = 2 ==> second derivative weights are computed |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
296 // id = 3 ==> gaussian weights are computed (in this |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
297 // case, the value of i is irrelevant) |
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 // Output parameters: |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
300 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
301 // dif1 : one dimensional vector containing the first derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
302 // of the node polynomial at the zeros |
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 // dif2 : one dimensional vector containing the second derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
305 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
306 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
307 // dif3 : one dimensional vector containing the third derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
308 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
309 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
310 // vect : one dimensional vector of computed weights |
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 static void |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
313 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
|
314 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
|
315 double *dif2, double *dif3, double *root, double *vect) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
316 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
317 assert (n0 == 0 || n0 == 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
318 assert (n1 == 0 || n1 == 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
319 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
320 octave_idx_type nt = n + n0 + n1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
321 |
25247
db31e068f4db
Rewrite incorrect assert statement in colloc calculation (bug #53653)
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
322 assert (nt >= 1); |
10603
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 assert (id == 1 || id == 2 || 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 if (id != 3) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
327 assert (i >= 0 && i < nt); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
328 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
329 // Evaluate discretization matrices and Gaussian quadrature weights. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
330 // Quadrature weights are normalized to sum to one. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
331 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
332 if (id != 3) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
333 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
334 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
|
335 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
336 if (j == i) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
337 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
338 if (id == 1) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
339 vect[i] = dif2[i] / dif1[i] / 2.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
340 else |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
341 vect[i] = dif3[i] / dif1[i] / 3.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
342 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
343 else |
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 double y = root[i] - root[j]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
346 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
347 vect[j] = dif1[i] / dif1[j] / y; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
348 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
349 if (id == 2) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
350 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
|
351 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
352 } |
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 else |
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 double y = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
357 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
358 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
|
359 { |
21568
3d60ed163b70
maint: Eliminate bad spacing around '='.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
360 double x = root[j]; |
10603
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 double ax = x * (1.0 - x); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
363 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
364 if (n0 == 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
365 ax = ax / x / x; |
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 if (n1 == 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
368 ax = ax / (1.0 - x) / (1.0 - x); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
369 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
370 vect[j] = ax / (dif1[j] * dif1[j]); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
371 |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
372 y += vect[j]; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
373 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
374 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
375 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
|
376 vect[j] = vect[j] / y; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
377 } |
3 | 378 } |
379 | |
380 // Error handling. | |
381 | |
382 void | |
23449
c763214a8260
maint: Use convention 'int *x' for naming pointers.
Rik <rik@octave.org>
parents:
23443
diff
changeset
|
383 CollocWt::error (const char *msg) |
3 | 384 { |
20726
25d676f9619c
Preface error() messages with name of function when possible.
Rik <rik@octave.org>
parents:
20230
diff
changeset
|
385 (*current_liboctave_error_handler) ("CollocWt: fatal error '%s'", msg); |
3 | 386 } |
387 | |
388 CollocWt& | |
389 CollocWt::set_left (double val) | |
390 { | |
391 if (val >= rb) | |
21055
5e00ed38a58b
maint: Replace if/error/else paradigm with just if/error.
Rik <rik@octave.org>
parents:
20726
diff
changeset
|
392 error ("CollocWt: left bound greater than right bound"); |
3 | 393 |
394 lb = val; | |
395 initialized = 0; | |
396 return *this; | |
397 } | |
398 | |
399 CollocWt& | |
400 CollocWt::set_right (double val) | |
401 { | |
402 if (val <= lb) | |
21055
5e00ed38a58b
maint: Replace if/error/else paradigm with just if/error.
Rik <rik@octave.org>
parents:
20726
diff
changeset
|
403 error ("CollocWt: right bound less than left bound"); |
3 | 404 |
405 rb = val; | |
406 initialized = 0; | |
407 return *this; | |
408 } | |
409 | |
410 void | |
411 CollocWt::init (void) | |
412 { | |
1360 | 413 // Check for possible errors. |
3 | 414 |
415 double wid = rb - lb; | |
416 if (wid <= 0.0) | |
227 | 417 { |
20726
25d676f9619c
Preface error() messages with name of function when possible.
Rik <rik@octave.org>
parents:
20230
diff
changeset
|
418 error ("CollocWt: width less than or equal to zero"); |
227 | 419 return; |
420 } | |
3 | 421 |
5275 | 422 octave_idx_type nt = n + inc_left + inc_right; |
1870 | 423 |
3 | 424 if (nt < 0) |
21055
5e00ed38a58b
maint: Replace if/error/else paradigm with just if/error.
Rik <rik@octave.org>
parents:
20726
diff
changeset
|
425 error ("CollocWt: total number of collocation points less than zero"); |
3 | 426 else if (nt == 0) |
427 return; | |
428 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
429 Array<double> dif1 (dim_vector (nt, 1)); |
1938 | 430 double *pdif1 = dif1.fortran_vec (); |
431 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
432 Array<double> dif2 (dim_vector (nt, 1)); |
1938 | 433 double *pdif2 = dif2.fortran_vec (); |
434 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
435 Array<double> dif3 (dim_vector (nt, 1)); |
1938 | 436 double *pdif3 = dif3.fortran_vec (); |
437 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
438 Array<double> vect (dim_vector (nt, 1)); |
1938 | 439 double *pvect = vect.fortran_vec (); |
3 | 440 |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
441 r.resize (nt, 1); |
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
442 q.resize (nt, 1); |
3 | 443 A.resize (nt, nt); |
444 B.resize (nt, nt); | |
445 | |
446 double *pr = r.fortran_vec (); | |
447 | |
1360 | 448 // Compute roots. |
3 | 449 |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
450 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
|
451 error ("jcobi: newton iteration failed"); |
3 | 452 |
5275 | 453 octave_idx_type id; |
3 | 454 |
1360 | 455 // First derivative weights. |
3 | 456 |
457 id = 1; | |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
458 for (octave_idx_type i = 0; i < nt; i++) |
3 | 459 { |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
460 dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect); |
3 | 461 |
5275 | 462 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
|
463 A(i,j) = vect(j); |
3 | 464 } |
465 | |
1360 | 466 // Second derivative weights. |
3 | 467 |
468 id = 2; | |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
469 for (octave_idx_type i = 0; i < nt; i++) |
3 | 470 { |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
471 dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect); |
3 | 472 |
5275 | 473 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
|
474 B(i,j) = vect(j); |
3 | 475 } |
476 | |
1360 | 477 // Gaussian quadrature weights. |
3 | 478 |
479 id = 3; | |
480 double *pq = q.fortran_vec (); | |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
481 dfopr (n, inc_left, inc_right, id, id, pdif1, pdif2, pdif3, pr, pq); |
3 | 482 |
483 initialized = 1; | |
484 } | |
485 | |
3504 | 486 std::ostream& |
487 operator << (std::ostream& os, const CollocWt& a) | |
3 | 488 { |
489 if (a.left_included ()) | |
490 os << "left boundary is included\n"; | |
491 else | |
492 os << "left boundary is not included\n"; | |
493 | |
494 if (a.right_included ()) | |
495 os << "right boundary is included\n"; | |
496 else | |
497 os << "right boundary is not included\n"; | |
498 | |
499 os << "\n"; | |
500 | |
23807
336f89b6208b
Use character literals 'c' rather than string literals "c" when possible.
Rik <rik@octave.org>
parents:
23649
diff
changeset
|
501 os << a.Alpha << ' ' << a.Beta << "\n\n" |
3 | 502 << a.r << "\n\n" |
503 << a.q << "\n\n" | |
504 << a.A << "\n" | |
505 << a.B << "\n"; | |
506 | |
507 return os; | |
508 } |