23 #if defined (HAVE_CONFIG_H) 76 double y =
x - root[j];
78 dif3[
i] =
y * dif3[
i] + 3.0 * dif2[
i];
79 dif2[
i] =
y * dif2[
i] + 2.0 * dif1[
i];
80 dif1[
i] =
y * dif1[
i];
148 double alpha,
double beta,
double *dif1,
double *dif2,
149 double *dif3,
double *root)
151 assert (n0 == 0 || n0 == 1);
152 assert (n1 == 0 || n1 == 1);
161 double ab = alpha + beta;
162 double ad = beta - alpha;
163 double ap = beta * alpha;
165 dif1[0] = (ad / (ab + 2.0) + 1.0) / 2.0;
173 double z = ab + 2 * z1;
175 dif1[
i] = (ab * ad / z / (z + 2.0) + 1.0) / 2.0;
178 dif2[
i] = (ab + ap + z1) / z / z / (z + 1.0);
182 double y = z1 * (ab + z1);
184 dif2[
i] =
y / z / (z - 1.0);
209 double xp = (dif1[j] -
x) * xn - dif2[j] * xd;
210 double xp1 = (dif1[j] -
x) * xn1 - dif2[j] * xd1 - xn;
224 zc -= z / (
x - root[j-1]);
236 if (
std::abs (z) <= 100 * std::numeric_limits<double>::epsilon ())
241 x += std::sqrt (std::numeric_limits<double>::epsilon ());
257 dif (nt, root, dif1, dif2, dif3);
311 double *dif2,
double *dif3,
double *root,
double *vect)
313 assert (n0 == 0 || n0 == 1);
314 assert (n1 == 0 || n1 == 1);
320 assert (
id == 1 ||
id == 2 ||
id == 3);
323 assert (
i >= 0 &&
i < nt);
335 vect[
i] = dif2[
i] / dif1[
i] / 2.0;
337 vect[
i] = dif3[
i] / dif1[
i] / 3.0;
341 double y = root[
i] - root[j];
343 vect[j] = dif1[
i] / dif1[j] /
y;
346 vect[j] = vect[j] * (dif2[
i] / dif1[
i] - 2.0 /
y);
358 double ax =
x * (1.0 -
x);
364 ax = ax / (1.0 -
x) / (1.0 -
x);
366 vect[j] = ax / (dif1[j] * dif1[j]);
372 vect[j] = vect[j] /
y;
381 (*current_liboctave_error_handler) (
"CollocWt: fatal error '%s'", msg);
388 error (
"CollocWt: left bound greater than right bound");
399 error (
"CollocWt: right bound less than left bound");
411 double wid =
rb -
lb;
414 error (
"CollocWt: width less than or equal to zero");
421 error (
"CollocWt: total number of collocation points less than zero");
447 error (
"jcobi: newton iteration failed");
485 if (
a.left_included ())
486 os <<
"left boundary is included\n";
488 os <<
"left boundary is not included\n";
490 if (
a.right_included ())
491 os <<
"right boundary is included\n";
493 os <<
"right boundary is not included\n";
497 os <<
a.Alpha <<
' ' <<
a.Beta <<
"\n\n"
CollocWt & set_right(double val)
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
CollocWt & set_left(double val)
identity matrix If supplied two scalar respectively For allows like xample val
const T * fortran_vec(void) const
std::ostream & operator<<(std::ostream &os, const CollocWt &a)
static void dif(octave_idx_type nt, double *root, double *dif1, double *dif2, double *dif3)
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
void error(const char *msg)
static void dfopr(octave_idx_type n, octave_idx_type n0, octave_idx_type n1, octave_idx_type i, octave_idx_type id, double *dif1, double *dif2, double *dif3, double *root, double *vect)
the element is set to zero In other the statement xample y
Vector representing the dimensions (size) of an Array.
octave_idx_type inc_right
F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE * x
static bool jcobi(octave_idx_type n, octave_idx_type n0, octave_idx_type n1, double alpha, double beta, double *dif1, double *dif2, double *dif3, double *root)
void resize(octave_idx_type n, const double &rfv=0)