XPost: sci.math, sci.math.num-analysis
From: spellucci@fb04633.mathematik.tu-darmstadt.de
Bart Goddard writes:
>pnachtwey wrote in news:d47c3576-ea88-4d27-
>be75-18dd8ddf6033@w24g2000prd.googlegroups.com:
>
>> I need to determine if there is a zero or root between two points
>for
>> a polynomial y=f(x). f(x) is a fourth order polynomial and the
>> 'forbidden zone' is between x0 to x1 where x0 is always at 0. I
>> don't need to know where the root is. I just need to quickly find
>if
>> there is a root in the 'forbidden zone'. What makes this a little
>> harder is that there will normally be roots at x0 and x1. I can
>find
>> the number of positive roots by counting the sign changes of the
>> coefficients of the polynomial. I can offset the time by x1-eps
>and
>> find the number of roots again. If the count is different then I
>know
>> a root exists between 0 and x1.
>
>If x0 and x1 are, in fact, roots, you can divide g(x) by
>x-x0 and x-x1, leaving a quadratic, which is easily solved
>for the other two roots.
for a fourth order polynomial there is an ''direct''
solution formula for all roots. well, this is not quite
super w.r.t. to ist roundoff behaviour, but fir this purpose it
should work. In netlib/toms there is a code (pdf file as copy of
a CACM publication) for this.
and here a -code (frm my annotations)
/*
* Roots3And4.c
*
* Utility functions to find cubic and quartic roots,
* coefficients are passed like this:
*
* c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0
*
* The functions return the number of non-complex roots and
* put the values into the s array.
*
* Author: Jochen Schwarze (schwarze@isa.de)
*
* Jan 26, 1990 Version for Graphics Gems
* Oct 11, 1990 Fixed sign problem for negative q's in SolveQuartic
* (reported by Mark Podlipec),
* Old-style function definitions,
* IsZero() as a macro
* Nov 23, 1990 Some systems do not declare acos() and cbrt() in
* , though the functions exist in the library.
* If large coefficients are used, EQN_EPS should be
* reduced considerably (e.g. to 1E-30), results will be
* correct but multiple roots might be reported more
* than once.
*/
#include
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
extern double sqrt(), cbrt(), cos(), acos();
/* epsilon surrounding for near zero values */
#define EQN_EPS 1e-9
#define IsZero(x) ((x) > -EQN_EPS && (x) < EQN_EPS)
#ifdef NOCBRT
#define cbrt(x) ((x) > 0.0 ? pow((double)(x), 1.0/3.0) : \
((x) < 0.0 ? -pow((double)-(x), 1.0/3.0) : 0.0))
#endif
int SolveQuadric(c, s)
double c[ 3 ];
double s[ 2 ];
{
double p, q, D;
/* normal form: x^2 + px + q = 0 */
p = c[ 1 ] / (2 * c[ 2 ]);
q = c[ 0 ] / c[ 2 ];
D = p * p - q;
if (IsZero(D))
{
s[ 0 ] = - p;
return 1;
}
else if (D < 0)
{
return 0;
}
else if (D > 0)
{
double sqrt_D = sqrt(D);
s[ 0 ] = sqrt_D - p;
s[ 1 ] = - sqrt_D - p;
return 2;
}
}
int SolveCubic(c, s)
double c[ 4 ];
double s[ 3 ];
{
int i, num;
double sub;
double A, B, C;
double sq_A, p, q;
double cb_p, D;
/* normal form: x^3 + Ax^2 + Bx + C = 0 */
A = c[ 2 ] / c[ 3 ];
B = c[ 1 ] / c[ 3 ];
C = c[ 0 ] / c[ 3 ];
/* substitute x = y - A/3 to eliminate quadric term:
x^3 +px + q = 0 */
sq_A = A * A;
p = 1.0/3 * (- 1.0/3 * sq_A + B);
q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);
/* use Cardano's formula */
cb_p = p * p * p;
D = q * q + cb_p;
if (IsZero(D))
{
if (IsZero(q)) /* one triple solution */
{
s[ 0 ] = 0;
num = 1;
}
else /* one single and one double solution */
{
double u = cbrt(-q);
s[ 0 ] = 2 * u;
s[ 1 ] = - u;
num = 2;
}
}
else if (D < 0) /* Casus irreducibilis: three real solutions */
{
double phi = 1.0/3 * acos(-q / sqrt(-cb_p));
double t = 2 * sqrt(-p);
s[ 0 ] = t * cos(phi);
s[ 1 ] = - t * cos(phi + M_PI / 3);
s[ 2 ] = - t * cos(phi - M_PI / 3);
num = 3;
}
else /* one real solution */
{
double sqrt_D = sqrt(D);
double u = cbrt(sqrt_D - q);
double v = - cbrt(sqrt_D + q);
s[ 0 ] = u + v;
num = 1;
}
/* resubstitute */
sub = 1.0/3 * A;
for (i = 0; i < num; ++i)
s[ i ] -= sub;
return num;
}
int SolveQuartic(c, s)
double c[ 5 ];
double s[ 4 ];
{
double coeffs[ 4 ];
double z, u, v, sub;
double A, B, C, D;
double sq_A, p, q, r;
int i, num;
/* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */
A = c[ 3 ] / c[ 4 ];
B = c[ 2 ] / c[ 4 ];
C = c[ 1 ] / c[ 4 ];
D = c[ 0 ] / c[ 4 ];
/* substitute x = y - A/4 to eliminate cubic term:
x^4 + px^2 + qx + r = 0 */
sq_A = A * A;
p = - 3.0/8 * sq_A + B;
q = 1.0/8 * sq_A * A - 1.0/2 * A * B + C;
r = - 3.0/256*sq_A*sq_A + 1.0/16*sq_A*B - 1.0/4*A*C + D;
if (IsZero(r))
{
/* no absolute term: y(y^3 + py + q) = 0 */
coeffs[ 0 ] = q;
coeffs[ 1 ] = p;
coeffs[ 2 ] = 0;
coeffs[ 3 ] = 1;
num = SolveCubic(coeffs, s);
s[ num++ ] = 0;
}
else
{
/* solve the resolvent cubic ... */
coeffs[ 0 ] = 1.0/2 * r * p - 1.0/8 * q * q;
coeffs[ 1 ] = - r;
coeffs[ 2 ] = - 1.0/2 * p;
coeffs[ 3 ] = 1;
(void) SolveCubic(coeffs, s);
/* ... and take the one real solution ... */
z = s[ 0 ];
/* ... to build two quadric equations */
u = z * z - r;
v = 2 * z - p;
if (IsZero(u))
u = 0;
else if (u > 0)
u = sqrt(u);
else
return 0;
if (IsZero(v))
v = 0;
else if (v > 0)
v = sqrt(v);
else
return 0;
coeffs[ 0 ] = z - u;
coeffs[ 1 ] = q < 0 ? -v : v;
coeffs[ 2 ] = 1;
num = SolveQuadric(coeffs, s);
coeffs[ 0 ]= z + u;
coeffs[ 1 ] = q < 0 ? v : -v;
coeffs[ 2 ] = 1;
num += SolveQuadric(coeffs, s + num);
}
/* resubstitute */
sub = 1.0/4 * A;
for (i = 0; i < num; ++i)
s[ i ] -= sub;
return num;
}
hope it works
Peter
--- SoupGate-Win32 v1.05
* Origin: you cannot sedate... all the things you hate (1:229/2)
|
