23 #if defined (HAVE_CONFIG_H)
73 double y = x - root[j];
75 dif3[
i] = y * dif3[
i] + 3.0 * dif2[
i];
76 dif2[
i] = y * dif2[
i] + 2.0 * dif1[
i];
77 dif1[
i] = y * dif1[
i];
145 double alpha,
double beta,
double *dif1,
double *dif2,
146 double *dif3,
double *root)
148 assert (n0 == 0 || n0 == 1);
149 assert (n1 == 0 || n1 == 1);
158 double ab = alpha + beta;
159 double ad = beta - alpha;
160 double ap = beta * alpha;
162 dif1[0] = (ad / (ab + 2.0) + 1.0) / 2.0;
170 double z = ab + 2 * z1;
172 dif1[
i] = (ab * ad / z / (z + 2.0) + 1.0) / 2.0;
175 dif2[
i] = (ab + ap + z1) / z / z / (z + 1.0);
179 double y = z1 * (ab + z1);
181 dif2[
i] = y / z / (z - 1.0);
206 double xp = (dif1[j] -
x) * xn - dif2[j] * xd;
207 double xp1 = (dif1[j] -
x) * xn1 - dif2[j] * xd1 - xn;
221 zc -= z / (x - root[j-1]);
233 if (
std::abs (z) <= 100 * std::numeric_limits<double>::epsilon ())
238 x += sqrt (std::numeric_limits<double>::epsilon ());
254 dif (nt, root, dif1, dif2, dif3);
308 double *dif2,
double *dif3,
double *root,
double *vect)
310 assert (n0 == 0 || n0 == 1);
311 assert (n1 == 0 || n1 == 1);
317 assert (
id == 1 ||
id == 2 ||
id == 3);
320 assert (i >= 0 && i < nt);
332 vect[
i] = dif2[
i] / dif1[
i] / 2.0;
334 vect[
i] = dif3[
i] / dif1[
i] / 3.0;
338 double y = root[
i] - root[j];
340 vect[j] = dif1[
i] / dif1[j] /
y;
343 vect[j] = vect[j] * (dif2[
i] / dif1[
i] - 2.0 /
y);
355 double ax = x * (1.0 -
x);
361 ax = ax / (1.0 -
x) / (1.0 - x);
363 vect[j] = ax / (dif1[j] * dif1[j]);
369 vect[j] = vect[j] / y;
378 (*current_liboctave_error_handler) (
"CollocWt: fatal error '%s'", msg);
385 error (
"CollocWt: left bound greater than right bound");
396 error (
"CollocWt: right bound less than left bound");
408 double wid =
rb -
lb;
411 error (
"CollocWt: width less than or equal to zero");
418 error (
"CollocWt: total number of collocation points less than zero");
444 error (
"jcobi: newton iteration failed");
474 dfopr (
n,
inc_left, inc_right,
id,
id, pdif1, pdif2, pdif3, pr, pq);
483 os <<
"left boundary is included\n";
485 os <<
"left boundary is not included\n";
488 os <<
"right boundary is included\n";
490 os <<
"right boundary is not included\n";
CollocWt & set_right(double val)
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
CollocWt & set_left(double val)
octave_idx_type right_included(void) const
identity matrix If supplied two scalar respectively For allows like xample val
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)
=val(i)}if ode{val(i)}occurs in table i
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
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
octave_idx_type left_included(void) const
const T * fortran_vec(void) const
Vector representing the dimensions (size) of an Array.
octave_idx_type inc_right
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x
static bool jcobi(octave_idx_type n, octave_idx_type n0, octave_idx_type n1, double alpha, double beta, double *dif1, double *dif2, double *dif3, double *root)
void resize(octave_idx_type n, const double &rfv=0)