Mercurial > octave
annotate liboctave/numeric/CollocWt.cc @ 25438:cb1606f78f6b
prefer <istream>, <ostream>, or <iosfwd> to <iostream> where possible
Using <iostream> brings with it a static initializer for the std::cin,
std::cout, and std::cerr streams. In most cases they are not needed
and should be avoided if possible.
Files affected:
build-aux/mk-opts.pl
libgui/qterminal/libqterminal/win32/QWinTerminalImpl.cpp
libinterp/corefcn/__dsearchn__.cc
libinterp/corefcn/c-file-ptr-stream.cc
libinterp/corefcn/c-file-ptr-stream.h
libinterp/corefcn/daspk.cc
libinterp/corefcn/dasrt.cc
libinterp/corefcn/dassl.cc
libinterp/corefcn/defaults.cc
libinterp/corefcn/defun.cc
libinterp/corefcn/file-io.cc
libinterp/corefcn/ft-text-renderer.cc
libinterp/corefcn/gl-render.cc
libinterp/corefcn/help.cc
libinterp/corefcn/ls-ascii-helper.cc
libinterp/corefcn/ls-hdf5.cc
libinterp/corefcn/ls-hdf5.h
libinterp/corefcn/ls-mat-ascii.cc
libinterp/corefcn/ls-mat4.cc
libinterp/corefcn/ls-mat5.cc
libinterp/corefcn/ls-oct-binary.cc
libinterp/corefcn/ls-oct-text.cc
libinterp/corefcn/lsode.cc
libinterp/corefcn/oct-iostrm.cc
libinterp/corefcn/oct-procbuf.cc
libinterp/corefcn/oct-stdstrm.h
libinterp/corefcn/procstream.cc
libinterp/corefcn/procstream.h
libinterp/corefcn/quad.cc
libinterp/corefcn/symscope.h
libinterp/corefcn/symtab.h
libinterp/corefcn/toplev.cc
libinterp/corefcn/urlwrite.cc
libinterp/corefcn/utils.cc
libinterp/corefcn/zfstream.cc
libinterp/dldfcn/__ode15__.cc
libinterp/dldfcn/convhulln.cc
libinterp/octave-value/ov-base-diag.cc
libinterp/octave-value/ov-base-int.cc
libinterp/octave-value/ov-base-mat.cc
libinterp/octave-value/ov-base-scalar.cc
libinterp/octave-value/ov-base-sparse.cc
libinterp/octave-value/ov-base.cc
libinterp/octave-value/ov-bool-mat.cc
libinterp/octave-value/ov-bool-sparse.cc
libinterp/octave-value/ov-bool.cc
libinterp/octave-value/ov-cell.cc
libinterp/octave-value/ov-ch-mat.cc
libinterp/octave-value/ov-class.cc
libinterp/octave-value/ov-colon.cc
libinterp/octave-value/ov-complex.cc
libinterp/octave-value/ov-cs-list.cc
libinterp/octave-value/ov-cx-mat.cc
libinterp/octave-value/ov-cx-sparse.cc
libinterp/octave-value/ov-fcn-handle.cc
libinterp/octave-value/ov-fcn-inline.cc
libinterp/octave-value/ov-float.cc
libinterp/octave-value/ov-flt-complex.cc
libinterp/octave-value/ov-flt-cx-mat.cc
libinterp/octave-value/ov-flt-re-mat.cc
libinterp/octave-value/ov-int16.cc
libinterp/octave-value/ov-int32.cc
libinterp/octave-value/ov-int64.cc
libinterp/octave-value/ov-int8.cc
libinterp/octave-value/ov-java.cc
libinterp/octave-value/ov-range.cc
libinterp/octave-value/ov-re-mat.cc
libinterp/octave-value/ov-re-sparse.cc
libinterp/octave-value/ov-scalar.cc
libinterp/octave-value/ov-str-mat.cc
libinterp/octave-value/ov-struct.cc
libinterp/octave-value/ov-typeinfo.cc
libinterp/octave-value/ov-uint16.cc
libinterp/octave-value/ov-uint32.cc
libinterp/octave-value/ov-uint64.cc
libinterp/octave-value/ov-uint8.cc
libinterp/octave.cc
libinterp/parse-tree/bp-table.cc
libinterp/parse-tree/lex.h
libinterp/parse-tree/profiler.cc
libinterp/parse-tree/pt-arg-list.cc
libinterp/parse-tree/pt-array-list.cc
libinterp/parse-tree/pt-assign.cc
libinterp/parse-tree/pt-cell.cc
libinterp/parse-tree/pt-const.cc
libinterp/parse-tree/pt-eval.cc
libinterp/parse-tree/pt-exp.cc
libinterp/parse-tree/pt-fcn-handle.cc
libinterp/parse-tree/pt-jit.cc
libinterp/parse-tree/pt-pr-code.cc
libinterp/parse-tree/pt-tm-const.cc
libinterp/parse-tree/pt.cc
liboctave/array/Array.cc
liboctave/array/CColVector.cc
liboctave/array/CDiagMatrix.cc
liboctave/array/CMatrix.cc
liboctave/array/CNDArray.cc
liboctave/array/CRowVector.cc
liboctave/array/CSparse.cc
liboctave/array/DiagArray2.cc
liboctave/array/MArray.cc
liboctave/array/Range.cc
liboctave/array/Sparse.cc
liboctave/array/boolMatrix.cc
liboctave/array/boolSparse.cc
liboctave/array/chMatrix.cc
liboctave/array/dColVector.cc
liboctave/array/dDiagMatrix.cc
liboctave/array/dMatrix.cc
liboctave/array/dNDArray.cc
liboctave/array/dRowVector.cc
liboctave/array/dSparse.cc
liboctave/array/fCColVector.cc
liboctave/array/fCDiagMatrix.cc
liboctave/array/fCMatrix.cc
liboctave/array/fCNDArray.cc
liboctave/array/fCRowVector.cc
liboctave/array/fColVector.cc
liboctave/array/fDiagMatrix.cc
liboctave/array/fMatrix.cc
liboctave/array/fNDArray.cc
liboctave/array/fRowVector.cc
liboctave/array/idx-vector.cc
liboctave/numeric/CollocWt.cc
liboctave/numeric/eigs-base.cc
liboctave/system/file-ops.cc
liboctave/system/oct-time.cc
liboctave/util/cmd-hist.cc
liboctave/util/data-conv.cc
liboctave/util/data-conv.h
liboctave/util/file-info.cc
liboctave/util/lo-utils.cc
liboctave/util/lo-utils.h
liboctave/util/quit.cc
liboctave/util/str-vec.cc
liboctave/util/url-transfer.cc
liboctave/util/url-transfer.h
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Thu, 07 Jun 2018 10:11:54 -0400 |
parents | db31e068f4db |
children | 00f796120a6d |
rev | line source |
---|---|
3 | 1 /* |
2 | |
25054
6652d3823428
maint: Update copyright dates in all source files.
John W. Eaton <jwe@octave.org>
parents:
24534
diff
changeset
|
3 Copyright (C) 1993-2018 John W. Eaton |
3 | 4 |
5 This file is part of Octave. | |
6 | |
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23807
diff
changeset
|
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 |
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23807
diff
changeset
|
9 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
|
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 |
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23807
diff
changeset
|
19 <https://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 | |
23475
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23449
diff
changeset
|
27 #include <cassert> |
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23449
diff
changeset
|
28 #include <cmath> |
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23449
diff
changeset
|
29 |
23443
3f1bf237908b
maint: Eliminate <cfloat.h> header from liboctave files.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
30 #include <limits> |
25438
cb1606f78f6b
prefer <istream>, <ostream>, or <iosfwd> to <iostream> where possible
John W. Eaton <jwe@octave.org>
parents:
25247
diff
changeset
|
31 #include <ostream> |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
32 |
23475
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23449
diff
changeset
|
33 #include "Array.h" |
3 | 34 #include "CollocWt.h" |
227 | 35 #include "lo-error.h" |
23475
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23449
diff
changeset
|
36 #include "lo-mappers.h" |
3 | 37 |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
38 // 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
|
39 // 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
|
40 // Equation Models by Polynomial Approximation, Prentice-Hall (1978) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
41 // pages 418-420. |
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 // Translated to C++ by jwe. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
44 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
45 // 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
|
46 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
47 // n0 (alpha,beta) n1 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
48 // p (x) = (x) * p (x) * (1 - x) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
49 // nt n |
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 // 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
|
52 // 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
|
53 // nt = n + n0 + n1 |
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 // 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
|
56 // They may be computed using jcobi. |
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 static void |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
59 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
|
60 double *dif3) |
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 // Evaluate derivatives of node polynomial using recursion formulas. |
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 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
|
65 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
66 double x = root[i]; |
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 dif1[i] = 1.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
69 dif2[i] = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
70 dif3[i] = 0.0; |
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 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
|
73 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
74 if (j != i) |
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 double y = x - root[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 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
|
79 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
|
80 dif1[i] = y * dif1[i]; |
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 } |
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 } |
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 // Compute the zeros of the Jacobi polynomial. |
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 // (alpha,beta) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
89 // p (x) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
90 // n |
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 // 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
|
93 // polynomial |
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 // n0 (alpha,beta) n1 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
96 // p (x) = (x) * p (x) * (1 - x) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
97 // nt n |
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 // at the interpolation points. |
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 // 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
|
102 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
103 // Input parameters: |
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 // 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
|
106 // |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
107 // 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
|
108 // of interior interpolation points) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
109 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
110 // 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
|
111 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
112 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
113 // n0 = 0 ==> x = 0 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
114 // n0 = 1 ==> x = 0 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
115 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
116 // 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
|
117 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
118 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
119 // n1 = 0 ==> x = 1 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
120 // n1 = 1 ==> x = 1 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
121 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
122 // 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
|
123 // polynomial |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
124 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
125 // 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
|
126 // polynomial |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
127 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
128 // 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
|
129 // and Michelsen, pages 57 to 59. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
130 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
131 // Output parameters: |
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 // root : one dimensional vector containing on exit the |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
134 // 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
|
135 // interpolation routine |
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 // dif1 : one dimensional vector containing the first derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
138 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
139 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
140 // dif2 : one dimensional vector containing the second derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
141 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
142 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
143 // dif3 : one dimensional vector containing the third derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
144 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
145 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
146 static bool |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
147 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
|
148 double alpha, double beta, double *dif1, double *dif2, |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
149 double *dif3, double *root) |
3 | 150 { |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
151 assert (n0 == 0 || n0 == 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
152 assert (n1 == 0 || n1 == 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 octave_idx_type nt = n + n0 + n1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
155 |
25247
db31e068f4db
Rewrite incorrect assert statement in colloc calculation (bug #53653)
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
156 assert (nt >= 1); |
10603
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 // -- first evaluation of coefficients in recursion formulas. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
159 // -- recursion coefficients are stored in dif1 and dif2. |
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 double ab = alpha + beta; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
162 double ad = beta - alpha; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
163 double ap = beta * alpha; |
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 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
|
166 dif2[0] = 0.0; |
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 if (n >= 2) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
169 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
170 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
|
171 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
172 double z1 = i; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
173 double z = ab + 2 * z1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
174 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
175 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
|
176 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
177 if (i == 1) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
178 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
|
179 else |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
180 { |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
181 z *= z; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
182 double y = z1 * (ab + z1); |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
183 y *= (ap + y); |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
184 dif2[i] = y / z / (z - 1.0); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
185 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
186 } |
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 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
189 // 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
|
190 // determined roots. |
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 double x = 0.0; |
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 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
|
195 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
196 bool done = false; |
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 int k = 0; |
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 while (! done) |
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 double xd = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
203 double xn = 1.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
204 double xd1 = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
205 double xn1 = 0.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
206 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
207 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
|
208 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
209 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
|
210 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
|
211 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
212 xd = xn; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
213 xd1 = xn1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
214 xn = xp; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
215 xn1 = xp1; |
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 |
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
218 double zc = 1.0; |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
219 double z = xn / xn1; |
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 if (i != 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
222 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
223 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
|
224 zc -= z / (x - root[j-1]); |
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 |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
227 z /= zc; |
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
228 x -= z; |
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 // 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
|
231 // enough in all cases. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
232 |
21782
2aef506f3fec
use namespace for lo-mappers.h functions
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
233 if (++k > 100 || octave::math::isnan (z)) |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
234 return false; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
235 |
15220
61822c866ba1
use std::numeric_limits<T>::epsilon in C++ code
John W. Eaton <jwe@octave.org>
parents:
15018
diff
changeset
|
236 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
|
237 done = true; |
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 root[i] = x; |
23649
aabf6cd222ac
Use sqrt from C++ std library in liboctave.
Rik <rik@octave.org>
parents:
23475
diff
changeset
|
241 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
|
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 // 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
|
245 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
246 if (n0 != 0) |
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 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
|
249 root[i] = root[i-1]; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
250 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
251 root[0] = 0.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 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
254 if (n1 != 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
255 root[nt-1] = 1.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 dif (nt, root, dif1, dif2, dif3); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
258 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
259 return true; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
260 } |
3 | 261 |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
262 // Compute derivative weights for orthogonal collocation. |
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 // See Villadsen and Michelsen, pages 133-134, 419. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
265 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
266 // Input parameters: |
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 // 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
|
269 // |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
270 // 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
|
271 // of interior interpolation points) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
272 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
273 // 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
|
274 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
275 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
276 // n0 = 0 ==> x = 0 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
277 // n0 = 1 ==> x = 0 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
278 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
279 // 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
|
280 // interpolation point |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
281 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
282 // n1 = 0 ==> x = 1 is not included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
283 // n1 = 1 ==> x = 1 is included |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
284 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
285 // 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
|
286 // calculated |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
287 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
288 // id : indicator |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
289 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
290 // id = 1 ==> first derivative weights are computed |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
291 // id = 2 ==> second derivative weights are computed |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
292 // id = 3 ==> gaussian weights are computed (in this |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
293 // case, the value of i is irrelevant) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
294 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
295 // Output parameters: |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
296 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
297 // dif1 : one dimensional vector containing the first derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
298 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
299 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
300 // dif2 : one dimensional vector containing the second derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
301 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
302 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
303 // dif3 : one dimensional vector containing the third derivative |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
304 // of the node polynomial at the zeros |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
305 // |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
306 // vect : one dimensional vector of computed weights |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
307 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
308 static void |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
309 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
|
310 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
|
311 double *dif2, double *dif3, double *root, double *vect) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
312 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
313 assert (n0 == 0 || n0 == 1); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
314 assert (n1 == 0 || n1 == 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 octave_idx_type nt = n + n0 + n1; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
317 |
25247
db31e068f4db
Rewrite incorrect assert statement in colloc calculation (bug #53653)
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
318 assert (nt >= 1); |
10603
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 assert (id == 1 || id == 2 || id == 3); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
321 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
322 if (id != 3) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
323 assert (i >= 0 && i < nt); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
324 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
325 // Evaluate discretization matrices and Gaussian quadrature weights. |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
326 // Quadrature weights are normalized to sum to one. |
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 (id != 3) |
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 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
|
331 { |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
332 if (j == i) |
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 if (id == 1) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
335 vect[i] = dif2[i] / dif1[i] / 2.0; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
336 else |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
337 vect[i] = dif3[i] / dif1[i] / 3.0; |
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 else |
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 double y = root[i] - root[j]; |
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 vect[j] = dif1[i] / dif1[j] / y; |
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 if (id == 2) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
346 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
|
347 } |
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 } |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
350 else |
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 double y = 0.0; |
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 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
|
355 { |
21568
3d60ed163b70
maint: Eliminate bad spacing around '='.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
356 double x = root[j]; |
10603
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 double ax = x * (1.0 - x); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
359 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
360 if (n0 == 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
361 ax = ax / x / x; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
362 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
363 if (n1 == 0) |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
364 ax = ax / (1.0 - x) / (1.0 - x); |
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 vect[j] = ax / (dif1[j] * dif1[j]); |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
367 |
20230
e914b5399c67
Use in-place operators in C++ code where possible.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
368 y += vect[j]; |
10603
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 |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
371 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
|
372 vect[j] = vect[j] / y; |
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
373 } |
3 | 374 } |
375 | |
376 // Error handling. | |
377 | |
378 void | |
23449
c763214a8260
maint: Use convention 'int *x' for naming pointers.
Rik <rik@octave.org>
parents:
23443
diff
changeset
|
379 CollocWt::error (const char *msg) |
3 | 380 { |
20726
25d676f9619c
Preface error() messages with name of function when possible.
Rik <rik@octave.org>
parents:
20230
diff
changeset
|
381 (*current_liboctave_error_handler) ("CollocWt: fatal error '%s'", msg); |
3 | 382 } |
383 | |
384 CollocWt& | |
385 CollocWt::set_left (double val) | |
386 { | |
387 if (val >= rb) | |
21055
5e00ed38a58b
maint: Replace if/error/else paradigm with just if/error.
Rik <rik@octave.org>
parents:
20726
diff
changeset
|
388 error ("CollocWt: left bound greater than right bound"); |
3 | 389 |
390 lb = val; | |
391 initialized = 0; | |
392 return *this; | |
393 } | |
394 | |
395 CollocWt& | |
396 CollocWt::set_right (double val) | |
397 { | |
398 if (val <= lb) | |
21055
5e00ed38a58b
maint: Replace if/error/else paradigm with just if/error.
Rik <rik@octave.org>
parents:
20726
diff
changeset
|
399 error ("CollocWt: right bound less than left bound"); |
3 | 400 |
401 rb = val; | |
402 initialized = 0; | |
403 return *this; | |
404 } | |
405 | |
406 void | |
407 CollocWt::init (void) | |
408 { | |
1360 | 409 // Check for possible errors. |
3 | 410 |
411 double wid = rb - lb; | |
412 if (wid <= 0.0) | |
227 | 413 { |
20726
25d676f9619c
Preface error() messages with name of function when possible.
Rik <rik@octave.org>
parents:
20230
diff
changeset
|
414 error ("CollocWt: width less than or equal to zero"); |
227 | 415 return; |
416 } | |
3 | 417 |
5275 | 418 octave_idx_type nt = n + inc_left + inc_right; |
1870 | 419 |
3 | 420 if (nt < 0) |
21055
5e00ed38a58b
maint: Replace if/error/else paradigm with just if/error.
Rik <rik@octave.org>
parents:
20726
diff
changeset
|
421 error ("CollocWt: total number of collocation points less than zero"); |
3 | 422 else if (nt == 0) |
423 return; | |
424 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
425 Array<double> dif1 (dim_vector (nt, 1)); |
1938 | 426 double *pdif1 = dif1.fortran_vec (); |
427 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
428 Array<double> dif2 (dim_vector (nt, 1)); |
1938 | 429 double *pdif2 = dif2.fortran_vec (); |
430 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
431 Array<double> dif3 (dim_vector (nt, 1)); |
1938 | 432 double *pdif3 = dif3.fortran_vec (); |
433 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
434 Array<double> vect (dim_vector (nt, 1)); |
1938 | 435 double *pvect = vect.fortran_vec (); |
3 | 436 |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
437 r.resize (nt, 1); |
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
438 q.resize (nt, 1); |
3 | 439 A.resize (nt, nt); |
440 B.resize (nt, nt); | |
441 | |
442 double *pr = r.fortran_vec (); | |
443 | |
1360 | 444 // Compute roots. |
3 | 445 |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
446 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
|
447 error ("jcobi: newton iteration failed"); |
3 | 448 |
5275 | 449 octave_idx_type id; |
3 | 450 |
1360 | 451 // First derivative weights. |
3 | 452 |
453 id = 1; | |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
454 for (octave_idx_type i = 0; i < nt; i++) |
3 | 455 { |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
456 dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect); |
3 | 457 |
5275 | 458 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
|
459 A(i,j) = vect(j); |
3 | 460 } |
461 | |
1360 | 462 // Second derivative weights. |
3 | 463 |
464 id = 2; | |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
465 for (octave_idx_type i = 0; i < nt; i++) |
3 | 466 { |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
467 dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect); |
3 | 468 |
5275 | 469 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
|
470 B(i,j) = vect(j); |
3 | 471 } |
472 | |
1360 | 473 // Gaussian quadrature weights. |
3 | 474 |
475 id = 3; | |
476 double *pq = q.fortran_vec (); | |
10603
d909c4c14b63
convert villad functions to C++
John W. Eaton <jwe@octave.org>
parents:
10350
diff
changeset
|
477 dfopr (n, inc_left, inc_right, id, id, pdif1, pdif2, pdif3, pr, pq); |
3 | 478 |
479 initialized = 1; | |
480 } | |
481 | |
3504 | 482 std::ostream& |
483 operator << (std::ostream& os, const CollocWt& a) | |
3 | 484 { |
485 if (a.left_included ()) | |
486 os << "left boundary is included\n"; | |
487 else | |
488 os << "left boundary is not included\n"; | |
489 | |
490 if (a.right_included ()) | |
491 os << "right boundary is included\n"; | |
492 else | |
493 os << "right boundary is not included\n"; | |
494 | |
495 os << "\n"; | |
496 | |
23807
336f89b6208b
Use character literals 'c' rather than string literals "c" when possible.
Rik <rik@octave.org>
parents:
23649
diff
changeset
|
497 os << a.Alpha << ' ' << a.Beta << "\n\n" |
3 | 498 << a.r << "\n\n" |
499 << a.q << "\n\n" | |
500 << a.A << "\n" | |
501 << a.B << "\n"; | |
502 | |
503 return os; | |
504 } |