How to find a if there is a zero or root between two points quickly
pnachtwey
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.
I can't be the first to have this problem. There must be an easier
way.
What I have found is that the smaller the interval between x0 and x1
the smaller the chance of a zero appearing between x0 and x1. Duh.
What would be really nice is to have symbolic answer. That way the
maximum x1 can be calculated that will keep the root outside the
'forbidden zone'.
I have wxMaxima and Mathcad and I can find the symbolic roots but the
answers are very messy. I can also find the root by using Newtons or
the secant method. The numerical approaches are not desirable because
the time to do the calculation is non-deterministic, I really don't
need to know where the root is and they don't provide a clue as to
how small x1 must be to keep the roots out side the 'forbidden zone'.
This is for an embedded real time application so speed and
deterministic calculation times are desired.
Thanks
Peter Nachtwey
Rob Johnson
I think what you want is the Sturm Chain for the polynomial.
See <http://en.wikipedia.org/wiki/Sturm%27s_theorem>.
Rob Johnson <r...@trash.whim.org>
take out the trash before replying
to view any ASCII art, display article in a monospaced font
sttsc...@tesco.net
If you want to locate real roots in an interval,
have a look at Sturm's algorithm. It is a clever application of
Euclid's algorithm to a polynomial and its derivative.
Daniel Lichtblau
In your "normal case" you already know two roots of the quartic are x0
and x1. If you can (computationally) deflate by (x-x0)*(x-x1) then the
quadratic formula should be a fast way to decide whether either
remaining root is in the interval. And checking the discriminant first
will quickly rule out the possibility, when it is negative.
Daniel Lichtblau
Wolfram Research
pnachtwey
> In article <d47c3576-ea88-4d27-be75-18dd8ddf6...@w24g2000prd.googlegroups.com>,
easy once one knows the algorithm. In my application I know the y(x1)
and y'(x1) so there was just one more equation to compute in the Strum
chain. I have done this symbolically using wxMaxima so determining if
there is a zero between x0 and x1 is easy. It was so easy that I just
had to make it more difficult and raise the high bar. Now I want to
find the largest value of x1 that will not have a zero between x0 and
x1. In my application the coefficients for f(x) are calculated from
the y(x0) and y(x1) its derivatives so finding this value of x1 is
getting messy.
I am trying a less general method where set the derivatives of y at x1
to 0. This appears to work well with fewer computations but the
intervals between x0 and x1 seem to be a little bigger and it seems
there must be checks for special cases.
BTW, wxMaxima got a significant face lift and is much more usable that
it has been in the past.
Peter Nachtwey
pnachtwey
What I am really doing is trying to calculate a motion profile for
'undefined' conditions where the goal
is to get to a ending position, velocity and acceleration and nicely
as possible. To avoid overshooting the velocity must not have a zero
between the current time, t0, and the ending time t1. Everything is
fine until I get to the Sturm chain calculations. The symbolic answer
is right. The numeric answer is wrong. One can see that %o3 and %o4
are the negative of each other but when I calculate the Sturm chain
numerically this is not the case. If I plug the symbolic formulas
into Mathcad the correct numeric solutions can be found. Mathcad
doesn't do the polynomial division which is why I use wxmaxima for
these types of calculations.
I have included print out. The wxm file can be found
ftp://ftp.deltacompsys.com/public/Maxima/Calc5thOrder.wxm
/*
Calculate the coefficients for a fifth order polynomial that goes
between
point 0 and point 1. The position, velocity and acceleration ( PVA )
are
provided at point 0 and point 1 along with a time between the two
points.
Solve for c3,c4,c5 given in the intial and final PVA and the time,
t01,
between the two points.
(%i1) remvalue(all)$ /* Undefine the variables and equations so the
solutions will be symbolic */
eq0: x0+v0*t01+(a0/2)*t01^2+c3*t01^3+c4*t01^4+c5*t01^5=x1$ /*
Position */
eq1: v0+5*c5*t01^4+4*c4*t01^3+3*c3*t01^2+a0*t01=v1$ /*
Velocity */
eq2: 20*c5*t01^3+12*c4*t01^2+6*c3*t01+a0=a1$ /*
Acceleration */
eq3: solve([eq0,eq1,eq2],[c3,c4,c5]);
grind(eq3);
(%o5) [[c3=(20*x1-20*x0-8*t01*v1-12*t01*v0+(a1-3*a0)*t01^2)/
(2*t01^3),c4=-(30*x1-30*x0-14*t01*v1-16*t01*v0+(2*a1-3*a0)*t01^2)/
(2*t01^4),c5=(12*x1-12*x0-6*t01*v1-6*t01*v0+(a1-a0)*t01^2)/(2*t01^5)]]
[[c3 = (20*x1-20*x0-8*t01*v1-12*t01*v0+(a1-3*a0)*t01^2)/(2*t01^3),
c4 = -(30*x1-30*x0-14*t01*v1-16*t01*v0+(2*a1-3*a0)*t01^2)/(2*t01^4),
c5 = (12*x1-12*x0-6*t01*v1-6*t01*v0+(a1-a0)*t01^2)/(2*t01^5)]]$
(%o6) done
/*
Compute the polynomial coefficients c3,c4 and c5
(%i7) x0: 0$ /* Initial position */
v0: 0$ /* Initial velocity */
a0: 1000$ /* Initial acceleration */
x1: 1/2$ /* Final position */
v1: 0$ /* Final velocity */
a1: 0$ /* Final acceleration */
t01: 1/4$ /* Time between point 0 and point 1 */
/* Calculate the coefficients c3,c4, and c5 */
remvalue(c3,c4,c5)$
c3=(c3: (20*x1-20*x0-8*t01*v1-12*t01*v0+(a1-3*a0)*t01^2)/(2*t01^3));
c4=(c4: -(30*x1-30*x0-14*t01*v1-16*t01*v0+(2*a1-3*a0)*t01^2)/
(2*t01^4));
c5=(c5: (12*x1-12*x0-6*t01*v1-6*t01*v0+(a1-a0)*t01^2)/(2*t01^5));
(%o15) c3=-5680
(%o16) c4=22080
(%o17) c5=-28928
/*
Calculate the position, velocity and acceleration as a function of
time and plot them.
>> pva(t):=[x0+v0*t+(a0/2)*t^2+c3*t^3+c4*t^4+c5*t^5,
v0+5*c5*t^4+4*c4*t^3+3*c3*t^2+a0*t,
20*c5*t^3+12*c4*t^2+6*c3*t+a0]$
/* Plot the position, velocity and acceleration for all the motion
segments */
wxplot2d(pva(t)[1],[t,0,t01],[xlabel,"Time"],[ylabel,"Position"],
[legend,"Postion"]);
wxplot2d(pva(t)[2],[t,0,t01],[xlabel,"Time"],[ylabel,"Velocity"],
[legend,"Velocity"]);
wxplot2d(pva(t)[3],[t,0,t01],[xlabel,"Time"],[ylabel,"Acceleration"],
[legend,"Acceleration"]);
/*
If the acceleration at t0 is 0 there will be a pole at t0.
If the velocity at t0 is 0there will be a another pole at t0.
If the acceleration at t1 is 0 there will be a pole at t1.
If the velocity at t1 is 0 there will be another pole at t1.
Divide by (t-t0) to simplify the problem
(%i18) t1: t01$
eq10: divide(v0+5*c5*t^4+4*c4*t^3+3*c3*t^2+a0*t-v1,(t-t1));
eq11: divide(eq10[1],(t-t1));
eq12: solve(eq11[1],t);
(%o19) [-144640*t^3+52160*t^2-4000*t,0]
(%o20) [16000*t-144640*t^2,0]
(%o21) [t=25/226,t=0]
/*
Sturm chain Calculations
Clear memory and shift-enter to calculate the symbolic equations.
(%i1) eq20: x0+v0*t+(a0/2)*t^2+c3*t^3+c4*t^4+c5*t^5-x1$ /* Position,
we do care about this now */
eq21: diff(eq20,t,1)-v1; /* Velocity, find
the zeros between 0 and t1 */
eq22: diff(eq21,t,1)-a1; /* Acceleration */
eq23: -remainder(eq21,eq22);
eq24: -remainder(eq22,eq23);
(%o2) -v1+v0+5*c5*t^4+4*c4*t^3+3*c3*t^2+a0*t
(%o3) 20*c5*t^3+12*c4*t^2+6*c3*t-a1+a0
(%o4) -(c5*(20*v0-20*v1)+(30*c3*c5-12*c4^2)*t^2+((5*a1+15*a0)
*c5-6*c3*c4)*t+(a1-a0)*c4)/(20*c5)
(%o5) -20*c5*t^3-12*c4*t^2-6*c3*t+a1-a0
/*
Evauluate the Sturm chain at t0
(%i12) t0: 0$
eq21,t:t0;
eq22,t:t0;
eq23,t:t0;
eq24,t:t0;
(%o13) v0-v1
(%o14) a0-a1
(%o15) -(c5*(20*v0-20*v1)+(a1-a0)*c4)/(20*c5)
(%o16) a1-a0
/*
Evaluate the Sturm chain at t1 which is equal to t01
since t0=0
(%i33) t1: t01$
eq21,t:t1;
eq22,t:t1;
eq23,t:t1;
eq24,t:t1;
(%o34) 0
(%o35) 0
(%o36) 0
(%o37) 0
Now I show the numerical calculations. I would expect %o51 to be
-1000.
(%i42) eq20: x0+v0*t+(a0/2)*t^2+c3*t^3+c4*t^4+c5*t^5-x1$ /*
Position, we do care about this now */
eq21: diff(eq20,t,1)-v1; /* Velocity, find
the zeros between 0 and t1 */
eq22: diff(eq21,t,1)-a1; /* Acceleration */
eq23: -remainder(eq21,eq22);
eq24: -remainder(eq22,eq23);
(%o43) -144640*t^4+88320*t^3-17040*t^2+1000*t
(%o44) -578560*t^3+264960*t^2-34080*t+1000
(%o45) -(359760*t^2-124440*t+8625)/226
(%o46) -(4983300000*t-1245825000)/2247001
/*
Evauluate the Sturm chain at t0
(%i47) t0: 0$
eq21,t:t0;
eq22,t:t0;
eq23,t:t0;
eq24,t:t0;
(%o48) 0
(%o49) 1000
(%o50) -8625/226
(%o51) 1245825000/2247001
/*
Evaluate the Sturm chain at t1 which is equal to t01
since t0=0
(%i52) t1: t01$
eq21,t:t1;
eq22,t:t1;
eq23,t:t1;
eq24,t:t1;
(%o53) 0
(%o54) 0
(%o55) 0
(%o56) 0
(%i57) kill(all);
(%o0) done
Any help is appreciated. I am not really stuck I just want to know if
this is a bug or I am doing something wrong. Any other suggestions on
better ways to use wxMaxima appreciated too.
Danial Lichtblau, I am looking at your root idea because you are right
that most of the time the solution is trivial. I can tell how many
roots are at the end points by the initial and final conditions. It is
relatively easy to to make special case code for condition. If only
the initial velocity is non-zero then there is only one non-zero
root. This case happens a lot but usually a decel or accel parameter
determine the time between the points so the very simple cases are
handled in a different way. However, if the initial acceleration and
velocity are non zero and final velocity and acceleration are 0 then
your idea can save a lot of time. The solution still needs to be
symbolic not in terms of c3, c4 and c5 but in terms the initial and
final conditions.
The goal is to find the optimal time for t1 that isn't too long that
the motion profile overshoots and yet doesn't exceed the maximum
velocity and acceleration specified by the user. Believe it or not is
is the small moves where one is only going a new thousandths of an
inch and will never reach the maximum velocity and acceleration that
are a hard problems. Using lower order polynomials doesn't help
because two or there must be spliced together.
Thanks
Peter Nachtwey
Dana DeLouis
> answers are very messy.
Hi. Just throwing this out @ Mathematica.
This says there is one real root between 0 and 1.5 inclusive.
Don't know the underlying algorithm thou.
poly = Expand[(x - 5)*(x - 4)*(x - 3)*(x - 2)*(x - 1)]
x^5 - 15*x^4 + 85*x^3 - 225*x^2 + 274*x - 120
CountRoots[poly, {x, 0, 1.5}]
1
= = =
Dana DeLouis
Achava Nakhash, the Loving Snake
Perhaps I am being silly but if f(X_0) AND f(x_1) have opposite signs,
then there is a root between them. If they have the same sign, there
is no root between them. If f(x_0) = 0 and f(x_1) != 0, then f(x) has
a root between implies the derivative of f(x) has a root between
them. "This can rule out a root even though it can't show that there
is one.
Obviousloy this approach doesn't completely solve the problem which is
not very difficult as others have pointed out since you have a
quadratic. However it is computationally extremely rapid and can
quickly evaluate many of your cases.
Hth,\
Achava
Virgil
<0aaf3214-5dd3-41c2...@r37g2000prr.googlegroups.com>,
Your first conclusion is correct but you second is not.
For example, f(x) = x^2 -1 is positive at x = 2 and x = -2 but has two
roots between them.
pnachtwey
> > I need to determine if there is a zero or root between two points
> > answers are very messy.
>
> Hi. Just throwing this out @ Mathematica.
> This says there is one real root between 0 and 1.5 inclusive.
> Don't know the underlying algorithm thou.
>
> poly = Expand[(x - 5)*(x - 4)*(x - 3)*(x - 2)*(x - 1)]
>
> x^5 - 15*x^4 + 85*x^3 - 225*x^2 + 274*x - 120
>
> CountRoots[poly, {x, 0, 1.5}]
>
wondering why the results for the Sturm chain are different
symbolically and numerically.
Peter Nachtwey
pnachtwey
The velocity has a root between t0 and t1. If the difference between
t0 and t1 is too big there will be a zero. I am trying to find what
that difference can be given the initial and final position, velocity
and acceleration.
>
> Obviousloy this approach doesn't completely solve the problem which is
> not very difficult as others have pointed out since you have a
> quadratic. However it is computationally extremely rapid and can
> quickly evaluate many of your cases.
Two of the coefficients are known because we know the initial velocity
and acceleration. c3,c4, and c5 are known in terms of the initial and
final position, velocity and acceleration.
Peter Nachtwey
Achava Nakhash, the Loving Snake
> In article
> <0aaf3214-5dd3-41c2-b79c-4465b53af...@r37g2000prr.googlegroups.com>,
Yes, of course. I should avoid posting when I am up past my bedtime.
However, one could do a quick seive operation with a small partition
to rule out simple cases, and that might be a useful preprocessor.
Regards,
achava
Virgil
<2632d3c3-48f8-4ba8...@p2g2000prf.googlegroups.com>,
Also, for differentiable functions, if there are at least two roots in
an interval, the derivative must have a root in the interval between
them.
Peter Spellucci
In article <450b6ef8-eb19-46c3...@p36g2000prp.googlegroups.com>,
pnachtwey <pnac...@gmail.com> writes:
>On Feb 1, 9:41=A0pm, Dana DeLouis <delo...@bellsouth.net> wrote:
>> =A0> I need to determine if there is a zero or root between two points
>> =A0> I have wxMaxima and Mathcad and I can find the symbolic roots but th=
>e
>> =A0> answers are very messy.
>>
>> Hi. =A0Just throwing this out @ Mathematica.
>> This says there is one real root between 0 and 1.5 inclusive.
>> Don't know the underlying algorithm thou.
>>
>>
>> x^5 - 15*x^4 + 85*x^3 - 225*x^2 + 274*x - 120
>>
>> CountRoots[poly, {x, 0, 1.5}]
>>
>This is for a real time embedded application. Right now I am
>wondering why the results for the Sturm chain are different
>symbolically and numerically.
>
>Peter Nachtwey
because the Sturm chain algorithm (Euclid) is numerically unstable.
symbolically done, everything is o.k., but the roundoff effects
are very strong having the effect that, interpreted in backward analysis,
you deal with a considerably different polynomial.
hth
peter
Chip Eastham
Spellucci) wrote:
> In article <450b6ef8-eb19-46c3-87c8-12662f525...@p36g2000prp.googlegroups.com>,
However the example he presents involves integer
coefficients and integer roots, well separated,
without multiplicity > 1.
I've not plowed through his "numeric" code, but
because it is significantly more complicated
than the one "symbolic" call to CountRoots,
it has a much greater likelihood of a bug.
regards, chip
pnachtwey
Symbolically my Sturm chain is:
(%o14) v0-v1
(%o15) a0-a1
(%o16) -(c5*(20*v0-20*v1)+(a1-a0)*c4)/(20*c5)
(%o17) a1-a0
Numerically my Sturm chain is:
(%o50) 0
(%o51) 1000
(%o52) -8625/226
(%o53) 1245825000/2247001
for the numbers I provided above one can set that %o51 and %o53 should
have opposite signs but they don't numerically.
Peter Nachtwey
JCH
> > I have wxMaxima and Mathcad and I can find the symbolic roots but the
> > answers are very messy.
>
> Hi. Just throwing this out @ Mathematica.
> This says there is one real root between 0 and 1.5 inclusive.
> Don't know the underlying algorithm thou.
>
> poly = Expand[(x - 5)*(x - 4)*(x - 3)*(x - 2)*(x - 1)]
>
> x^5 - 15*x^4 + 85*x^3 - 225*x^2 + 274*x - 120
>
> CountRoots[poly, {x, 0, 1.5}]
>
> 1
>
--
Regards/Grüße Jan C. Hoffmann
Microsoft-kompatibel/optimiert für IE7+OE7
hagman
> "Dana DeLouis" <delo...@bellsouth.net> schrieb im Newsbeitragnews:1Gvhl.1029$vb....@bignews4.bellsouth.net...
>
>
>
> > > I need to determine if there is a zero or root between two points
> > > I have wxMaxima and Mathcad and I can find the symbolic roots but the
> > > answers are very messy.
>
> > Hi. Just throwing this out @ Mathematica.
> > This says there is one real root between 0 and 1.5 inclusive.
> > Don't know the underlying algorithm thou.
>
> > poly = Expand[(x - 5)*(x - 4)*(x - 3)*(x - 2)*(x - 1)]
>
> > x^5 - 15*x^4 + 85*x^3 - 225*x^2 + 274*x - 120
>
> > CountRoots[poly, {x, 0, 1.5}]
>
> > 1
>
> Are you sure?
>
> y = x^5 - 15*x^4 + 85*x^3 - 225*x^2 + 274*x - 120
>
> has 5 roots:
>
> x1=1
> x2=2
> x3=3
> x4=4
> x5=5
Erm, yes, wasn't this just his point?
He explicitly constructed poly to have these five roots
and indeed exactly one of these is between 0 and 1.5
Back to the topic I guess that Mathematica uses Sturm chains as well
for CountRoots - what could be better? It involves only derivatives
of polynomials and polynomial division
>
> * See also
> *http://home.arcor.de/janch/janch/_news/20090207-polynomial/
>
> --
> Regards/Grüße Jan C. Hoffmann
> Microsoft-kompatibel/optimiert für IE7+OE7
Before optimizing for a single user agent, I suggest getting rid of
the
more than 100 validation errors:
http://validator.w3.org/check?verbose=1&uri=http%3A%2F%2Fhome.arcor.de%2Fjanch%2Fjanch%2F_news%2F20090207-polynomial%2F
and
http://schneegans.de/sv/?url=http%3A%2F%2Fhome.arcor.de%2Fjanch%2Fjanch%2F_news%2F20090207-polynomial%2F
;)
Dana DeLouis
the link below.
I put something together with Mathematica.
Since I wasn't sure how long the "chain" was, I fixed-pointed it until
it ended with Indeterminate (0/0).
(Quiet was used so as to not display error messages)
Not tested very well, but I do see that it counts double root only once.
http://mathworld.wolfram.com/SturmFunction.html
SturmChain[equ_, var_] := Module[{},
v = {equ, D[equ, var]};
v = FixedPointList[FullSimplify[
{#[[2]],
#[[2]]* PolynomialQuotient[#[[1]], #[[2]], var] - #[[1]]}] &,
v] // Quiet;
v = DeleteCases[v, {_, Indeterminate}];
v[[All, 1]]
]
CountRoot[equ_, {x_, a_, b_}] := Module[{g},
g = SturmChain[equ, x];
g = {g /. x -> a, g /. x -> b};
g = (Count[Abs[Differences[Sign[#1]]], 2] &) /@ g;
First[Abs[Differences[g]]]
]
I counted a difference of -1 to +1 =2 as a change in sign.
There may be a better way.
poly1 = -1 - 3 x + x^5;
This is from the example on the above link:
SturmChain[poly1, x]
{
-1 - 3 x + x^5,
-3 + 5 x^4,
1 + (12 x)/5,
59083/20736
}
CountRoot[poly1, {x, -2, 0}]
2
CountRoot[poly1, {x, 0, 2}]
1
Here's an example of a double root.
For some reason, this algorithm counts 3 once.
The Mathematica built-in function counts 3 twice.
poly2 = Expand[(x - 5)*(x - 3)*(x - 3)*(x - 2)*(x - 1)]
-90 + 213 x - 184 x^2 + 74 x^3 - 14 x^4 + x^5
This is ok:
CountRoot[poly2, {x, 0, 2.5}]
2
This counts 3 once!
CountRoot[poly2, {x, 2.5, 3.5}]
1
The Mathematica built-in function counts 3 twice:
CountRoots[poly2, {x, 2.5, 3.5}]
2
I don't know how to improve on this.
Dana DeLouis
"CountRoots" from Mathematica is that maybe the op can find a similar
function in the math programs he is using:
> I have wxMaxima and Mathcad ...
Dana
rjf
in wxmaxima (a free open-source program that runs on windows,
linux, ...):
(%i3) ? nroots;
-- Function: nroots (<p>, <low>, <high>)
Returns the number of real roots of the real univariate
polynomial
<p> in the half-open interval `(<low>, <high>]'. The endpoints
of
the interval may be `minf' or `inf'. infinity and plus infinity.
`nroots' uses the method of Sturm sequences.
(%i1) p: x^10 - 2*x^4 + 1/2$
(%i2) nroots (p, -6, 9.1);
(%o2) 4
....
This program has been in Macsyma/Maxima probably for decades. If you
know in
advance the degree of the polynomial and most of the coefficients, and
some
other information and you are operating under some real-time
constraint, you
might not want to use this program. It presumably will tell you the
answer,
though the example uses a floating-point number which is a
questionable strategy.
9.1 is slightly less than 91/100. Closer to 9.09999847.
But the polynomial p given in the example has roots nowhere near 9.1
RJF
pnachtwey
written on a CAS before being implemented in C. The first thing I did
after being enlightened about the Sturm chain was to do a search an I
found the same example on the Wolfram site
(%i1) f0: x^5-3*x-1;
f1: diff(f0,x,1);
f2: -remainder(f0,f1);
f3: -remainder(f1,f2);
(%o1) x^5-3*x-1
(%o2) 5*x^4-3
(%o3) (12*x+5)/5
(%o4) 59083/20736
That was easy so I raised the high bar. Now the goal is to find the
maximum interval that has no zeros between the beginning and time.
wxMaxima appears to have a bug so that the numerical solution did not
work so I copied the equations to Mathcad and tried to solve for the
interval. The solution was very messy so I doubt this method can be
used in an embedded application. I have another solution that is easy
to implement and executes quickly but isn't optimal.
Peter Nachtwey
Now determine if there is a zero between t0 and t1 in the velocity.
This would indicate that the motion profile will overshoot the set
point.
First start with the position polynomial.
The position polynomial is only used to generate the velocity and
acceleration polynomial.
Only two divides are required to get a constant quotient.
(%i15) eq0: x0+v0*t+(a0/2)*t^2+c3*t^3+c4*t^4+c5*t^5-x1$ /* Position
*/
eq1: diff(eq0,t,1)-v1; /* Velocity, find
the zeros between 0 and t1 */
eq2: diff(eq1,t,1)-a1; /* Acceleration */
eq3: -remainder(eq1,eq2);
eq4: -remainder(eq2,eq3);
(%o16) -v1+v0+5*c5*t^4+4*c4*t^3+3*c3*t^2+a0*t
(%o17) 20*c5*t^3+12*c4*t^2+6*c3*t-a1+a0
(%o18) -(c5*(20*v0-20*v1)+(30*c3*c5-12*c4^2)*t^2+((5*a1+15*a0)
*c5-6*c3*c4)*t+(a1-a0)*c4)/(20*c5)
(%o19) -20*c5*t^3-12*c4*t^2-6*c3*t+a1-a0
hagman
> On Feb 7, 6:51 am, Dana DeLouis <delo...@bellsouth.net> wrote:> I forgot to mention that the only reason I mention this function
> > "CountRoots" from Mathematica is that maybe the op can find a similar
> > function in the math programs he is using:
>
> > > I have wxMaxima and Mathcad ...
>
> > Dana
>
> This is for an embedded application. Much of the code is first
> written on a CAS before being implemented in C. The first thing I did
> after being enlightened about the Sturm chain was to do a search an I
> found the same example on the Wolfram site
>
> (%i1) f0: x^5-3*x-1;
> f1: diff(f0,x,1);
> f2: -remainder(f0,f1);
> f3: -remainder(f1,f2);
> (%o1) x^5-3*x-1
> (%o2) 5*x^4-3
> (%o3) (12*x+5)/5
> (%o4) 59083/20736
>
> That was easy so I raised the high bar. Now the goal is to find the
> maximum interval that has no zeros between the beginning and time.
Erm, doesn't this simply mean you need to find the first root?
Phil Carmody
>
> "Dana DeLouis" <del...@bellsouth.net> schrieb im Newsbeitrag
> news:1Gvhl.1029$vb....@bignews4.bellsouth.net...
>> > I need to determine if there is a zero or root between two points
>> > I have wxMaxima and Mathcad and I can find the symbolic roots but the
>> > answers are very messy.
>>
>> Hi. Just throwing this out @ Mathematica.
>> This says there is one real root between 0 and 1.5 inclusive.
>> Don't know the underlying algorithm thou.
>>
>> poly = Expand[(x - 5)*(x - 4)*(x - 3)*(x - 2)*(x - 1)]
>>
>> x^5 - 15*x^4 + 85*x^3 - 225*x^2 + 274*x - 120
>>
>> CountRoots[poly, {x, 0, 1.5}]
>>
>> 1
>
>
> Are you sure?
>
> y = x^5 - 15*x^4 + 85*x^3 - 225*x^2 + 274*x - 120
>
> has 5 roots:
>
> x1=1
> x2=2
> x3=3
> x4=4
> x5=5
And how many of those are between 0 and 1.5?
Phil
--
I tried the Vista speech recognition by running the tutorial. I was
amazed, it was awesome, recognised every word I said. Then I said the
wrong word ... and it typed the right one. It was actually just
detecting a sound and printing the expected word! -- pbhj on /.
JCH
>>
>> "Dana DeLouis" <del...@bellsouth.net> schrieb im Newsbeitrag
>> news:1Gvhl.1029$vb....@bignews4.bellsouth.net...
>>> > I need to determine if there is a zero or root between two points
>>> > I have wxMaxima and Mathcad and I can find the symbolic roots but the
>>> > answers are very messy.
>>>
>>> Hi. Just throwing this out @ Mathematica.
>>> This says there is one real root between 0 and 1.5 inclusive.
>>> Don't know the underlying algorithm thou.
>>>
>>> poly = Expand[(x - 5)*(x - 4)*(x - 3)*(x - 2)*(x - 1)]
>>>
>>> x^5 - 15*x^4 + 85*x^3 - 225*x^2 + 274*x - 120
>>>
>>> CountRoots[poly, {x, 0, 1.5}]
>>>
>>> 1
>>
>>
>> Are you sure?
>>
>> y = x^5 - 15*x^4 + 85*x^3 - 225*x^2 + 274*x - 120
>>
>> has 5 roots:
>>
>> x1=1
>> x2=2
>> x3=3
>> x4=4
>> x5=5
>
> And how many of those are between 0 and 1.5?
>
c_1 = 274
c_2 = -225
c_3 = 85
c_4 = -15
c_5 = 1
c_6 = 0
c_7 = 0
c_8 = 0
c_9 = 0
Gordon Sande
The book "The Numerical Treatment of a Single Nonlinear Equation"
by Alston Householder in 1970 is really about solving polynomial roots.
It gives the usual Sturm's theorem as well as other facts. Like alternate
constructions of a Sturm's sequence as the GCD of the polynomial and its
derivative is the common but not the only possibility. It also gives the
Budan and Fourier theorem which is easier to compute but less precise
than the Sturm's theorem. Lots of stuff that never makes it into a book
with a single chapter on polynomial roots. Google throws up quite a bit
of interesting stuff for "Budan Fourier" that does not seem to show for
"Sturm theorem" or "Strum sequence".
Phil Carmody
Are you sure?
JCH
Dana DeLouis
Oh! I see. That's the long way of saying the Polynomial Remainder. :>(
ref:
http://mathworld.wolfram.com/SturmFunction.html
SturmChain[poly_, var_] :=
Module[{v, x},
x = var;
v = {poly, D[poly, x]};
While[Exponent[Last[v], x] > 0,
AppendTo[v, Simplify[-PolynomialRemainder[v[[-2]], v[[-1]], x]]]];
v]
equ = -1 - 3 x + x^5;
ans1 = SturmChain[equ, x]
{-1 - 3 x + x^5, -3 + 5 x^4, 1 + (12 x)/5, 59083/20736}
Symbolically...
equ = a x^5 + b x + c;
SturmChain[equ, x]
{c + b x + a x^5,
b + 5 a x^4,
-c - (4 b x)/5,
-b - (3125 a c^4)/(256 b^4)}
If I plug in the same values...
ans2 = % /. {a -> 1, b -> -3, c -> -1}
{-1 - 3 x + x^5, -3 + 5 x^4, 1 + (12 x)/5, 59083/20736}
I get the same answer.
ans1 == ans2
True
A full 5th degree symbolic polynomial was a very long mess!!
Dana
<snip>
Chip Eastham
> If solving as non-linear equation and using start
> value = 0 then you will find the 1st solution for
> sure for this example.
Perhaps what you have in mind is that Newton's
method, due to the 2nd derivative being negative
on [0,1], will converge to the "1st solution" if
the initial guess is x = 0. In general it's
possible for Newton's method to "overshoot" the
nearest root...
regards, chip
Dana DeLouis
I used Mathematica.
My function...
SturmChain[poly_, var_] := Module[{v, x},
x = var;
v = {poly, D[poly, x]};
While[Exponent[Last[v], x] > 0,
AppendTo[v,
Simplify[-PolynomialRemainder @@ Join[Take[v, -2], {x}]]
]];
v]
When the chain does not end in 0, there are no double roots:
sc = SturmChain[-1 - 3 x + x^5, x]
{-1 - 3 x + x^5, -3 + 5 x^4, 1 + (12 x)/5, 59083/20736}
I use expand so as to see the obvious solutions...
equ = (x - 1) (x - 2) (x - 3) (x - 3) (x - 4) // Expand
-72 + 174 x - 155 x^2 + 65 x^3 - 13 x^4 + x^5
sc = SturmChain[equ, x]
{-72 + 174 x - 155 x^2 + 65 x^3 - 13 x^4 + x^5,
174 - 310 x + 195 x^2 - 52 x^3 + 5 x^4,
2/25 (-231 + 275 x - 105 x^2 + 13 x^3),
25/169 (219 - 166 x + 31 x^2), (12168 (-3 + x))/24025, 0}
The chain ends in zero, so there are multiple solutions.
A plot shows all the equations go thru zero at x=3.
Plot[sc, {x, 1, 4}];
They are all zero!
sc /. x -> 3
{0, 0, 0, 0, 0, 0}
We use the second to last entry to determine the double roots.
sc[[-2]]
(12168 (-3 + x))/24025
From inspection, x=3 sets this to 0...we add 1 to the count, and
conclude there are two 3's.
equ = (x - 1) (x - 2) (x - 3) (x - 3) (x - 3) // Expand
-54 + 135 x - 126 x^2 + 56 x^3 - 12 x^4 + x^5
sc = SturmChain[equ, x]
{-54 + 135 x - 126 x^2 + 56 x^3 - 12 x^4 + x^5,
135 - 252 x + 168 x^2 - 48 x^3 + 5 x^4, 2/25 (-3 + x)^2 (-15 + 8 x),
75/64 (-3 + x)^2, 0}
Ends in zero, so multiple roots...
sc[[-2]]
75/64 (-3 + x)^2
From inspection, there are 2+1=3 3's.
Here's a more complicated one...
equ = (x - 2) (x - 2) (x - 4) (x - 4) (x - 4) // Expand
-256 + 448 x - 304 x^2 + 100 x^3 - 16 x^4 + x^5
sc = SturmChain[equ, x]
{-256 + 448 x - 304 x^2 + 100 x^3 - 16 x^4 + x^5,
448 - 608 x + 300 x^2 - 64 x^3 + 5 x^4, 24/25 (-4 + x)^2 (-2 + x), 0}
sc[[-2]]
24/25 (-4 + x)^2 (-2 + x)
From inspection, there are 2+1=3 4's, and 1+1=2 2's
Basically, we are solving the equation, and adding 1 to the count of
each solution.
Solve[% == 0] // Tally
{{{x -> 2}, 1}, {{x -> 4}, 2}}
Interesting subject... :>)
Dana DeLouis
pnachtwey
See an engineers approach to the problem:
(%i6) eq0: x0+v0*t+(a0/2)*t^2+c3*t^3+c4*t^4+c5*t^5-x1$ /* Position,
this is just a starting point */
eq1: diff(eq0,t,1); /* Velocity, find
the zeros between 0 and t1 */
eq3: -remainder(eq1,eq2);
eq4: -remainder(eq2,eq3);
(%o7) v0+5*c5*t^4+4*c4*t^3+3*c3*t^2+a0*t
(%o8) 20*c5*t^3+12*c4*t^2+6*c3*t+a0
(%o9) -(20*c5*v0+(30*c3*c5-12*c4^2)*t^2+(15*a0*c5-6*c3*c4)*t-a0*c4)/
(20*c5)
(%o10) -20*c5*t^3-12*c4*t^2-6*c3*t-a0
I don't care about the position. I included it just to diff() to
generate the formulas for the velocity and acceleration so the Sturm
chain starts with the velocity which is 4th order. Also the real
problem has conditions or constraints to my solution doesn't have to
work for all cases, just the cases I am interested in.
Now look at (%o7), this will evaluate to v0 when t= and to v1 when
t=1.
(%o8) Also evaluates to a0 and a1 and (%o10) evaluates to -a0 and -a1.
I v0,v1,a0 and a1 are known parameters that don't need to be
calculated. This leaves the (%o9) term. Evaluating (%o9) at t=0 and
t=t1 is not a problem. The first goal of the using the Sturm chain to
find if the is a zero between t=0 and t=t1 is easy.
I have found that if I change my formula just a little bit that if I
subtract v1 from the velocity equation and subtract a1 from the
acceleration term then the Sturm chain at time t1 evaluates to
[0;0;0;0] when v1=0 and a1=0 which is the normal case. Now I only
need to valuate the Sturm chain at t0 and look for a sign change
within the four values at t0.
%i14) eq0: x0+v0*t+(a0/2)*t^2+c3*t^3+c4*t^4+c5*t^5-x1$ /* Position,
this is just a starting point */
eq1: diff(eq0,t,1)-v1; /* Velocity, find
the zeros between 0 and t1 */
eq2: diff(eq1,t,1)-a1; /* Acceleration */
eq3: -remainder(eq1,eq2);
eq4: -remainder(eq2,eq3);
(%o15) -v1+v0+5*c5*t^4+4*c4*t^3+3*c3*t^2+a0*t
(%o16) 20*c5*t^3+12*c4*t^2+6*c3*t-a1+a0
(%o17) -(c5*(20*v0-20*v1)+(30*c3*c5-12*c4^2)*t^2+((5*a1+15*a0)
*c5-6*c3*c4)*t+(a1-a0)*c4)/(20*c5)
(%o18) -20*c5*t^3-12*c4*t^2-6*c3*t+a1-a0
That simplifies to
(%i19) eq1,t:0;
eq2,t:0;
eq3,t:0;
eq4,t:0;
(%o19) v0-v1
(%o20) a0-a1
(%o21) -(c5*(20*v0-20*v1)+(a1-a0)*c4)/(20*c5)
(%o22) a1-a0
No square root or divide required!!! I don't care about the real
value of (%o21), just the sign of the numerator and the denominator.
Remember this is for an embedded system where divisions and square
roots are time sinks. Evaluating this Sturm chain symbolically ahead
time is trivial and evaluating at run time is faster than factoring
out two roots and solving the quadratic equation which requires a
square root function.
Now for the hard part. As t1 gets smaller any root between t0 and t1
will move toward t1 or t1 will move to the root so that the root is no
longer between t0 and t1. I want to calculate t1 where root has just
moved to t1 and not between t0 and t1. Now things get messy because
c3,c4 and c5 are defined in terms of interval between t0 and t1 so
substitutions must be made for c3,c4,and c5 and then solve for t1.
This is the new problem, I raised the high bar.
Peter Nachtwey
JCH
Peter Spellucci
In article <b84bffca-423a-413a...@w39g2000prb.googlegroups.com>,
pnachtwey <pnac...@gmail.com> writes:
>On Feb 8, 8:30=A0am, Dana DeLouis <delo...@bellsouth.net> wrote:
>> =A0> =A0 #[[2]]* PolynomialQuotient[#[[1]], #[[2]], var] - #[[1]]}] &,
>>
>> Oh! =A0I see. =A0That's the long way of saying the Polynomial Remainder. =
>=A0:>(
>>
>> ref:http://mathworld.wolfram.com/SturmFunction.html
>>
>> SturmChain[poly_, var_] :=3D
>> =A0 Module[{v, x},
>> x =3D var;
>> v =3D {poly, D[poly, x]};
>> =A0 While[Exponent[Last[v], x] > 0,
>> =A0 =A0 AppendTo[v, Simplify[-PolynomialRemainder[v[[-2]], v[[-1]], x]]]]=
>;
>> v]
>>
>> equ =3D -1 - 3 x + x^5;
>>
>> ans1 =3D SturmChain[equ, x]
>>
>> {-1 - 3 x + x^5, -3 + 5 x^4, 1 + (12 x)/5, 59083/20736}
>>
>> Symbolically...
>>
>>
>> SturmChain[equ, x]
>>
>> {c + b x + a x^5,
>> =A0 -c - (4 b x)/5,
>> =A0 -b - (3125 a c^4)/(256 b^4)}
>>
>> If I plug in the same values...
>>
>>
>> {-1 - 3 x + x^5, -3 + 5 x^4, 1 + (12 x)/5, 59083/20736}
>>
>> I get the same answer.
>>
>>
>> True
>>
>> A full 5th degree symbolic polynomial was a very long mess!!
>>
>> Dana
>>
>> <snip>
>The Sturm chain for my problem is not a mess. That is the simple part.
>See an engineers approach to the problem:
>(%i6) eq0: x0+v0*t+(a0/2)*t^2+c3*t^3+c4*t^4+c5*t^5-x1$ /* Position,
>this is just a starting point */
>eq1: diff(eq0,t,1); /* Velocity, find
>the zeros between 0 and t1 */
>eq2: diff(eq1,t,1); /* Acceleration */
>eq3: -remainder(eq1,eq2);
>eq4: -remainder(eq2,eq3);
>
>(%o7) v0+5*c5*t^4+4*c4*t^3+3*c3*t^2+a0*t
>(%o8) 20*c5*t^3+12*c4*t^2+6*c3*t+a0
>(%o9) -(20*c5*v0+(30*c3*c5-12*c4^2)*t^2+(15*a0*c5-6*c3*c4)*t-a0*c4)/
>(20*c5)
>(%o10) -20*c5*t^3-12*c4*t^2-6*c3*t-a0
>
>I don't care about the position. I included it just to diff() to
>generate the formulas for the velocity and acceleration so the Sturm
>chain starts with the velocity which is 4th order. Also the real
>problem has conditions or constraints to my solution doesn't have to
>work for all cases, just the cases I am interested in.
>
>t=3D1.
>(%o8) Also evaluates to a0 and a1 and (%o10) evaluates to -a0 and -a1.
>I v0,v1,a0 and a1 are known parameters that don't need to be
>t=3Dt1 is not a problem. The first goal of the using the Sturm chain to
>find if the is a zero between t=3D0 and t=3Dt1 is easy.
>
>I have found that if I change my formula just a little bit that if I
>subtract v1 from the velocity equation and subtract a1 from the
>acceleration term then the Sturm chain at time t1 evaluates to
doesn't that mean that you want the smallest strictly positive root
of a (real) polynomial?
then look here:
http://numawww.mathematik.tu-darmstadt.de:80
in the section
"nichtlineare gls"
the algorithm
"reelle nullstellen reeller polynome.
(uses the complete hornerscheme for taylor-development at a point, then
maintains the quadratic but computes a bound for the second derivative on a
prospective interval such that its zeros underestimate the zeros of
the original polynomial, then computes the zeros of the quadratic and proceeds.)
hth
peter
Chip Eastham
[in part]
> Now things get messy because c3,c4 and c5 are
> defined in terms of interval between t0 and t1
> so substitutions must be made for c3,c4,and c5
> and then solve for t1.
> This is the new problem, I raised the high bar.
Hi, Peter:
You seem to be saying that the coefficients, and
thus the locations of roots, depend on interval
[t0,t1]. I suspect this is why you've not
responded to the point raised by hagman and
Peter Spelluci, that for a given polynomial,
the longest interval starting at the origin
which does not contain a root would end at
the location of the smallest positive real
root (the finding of which is probably
easier than the Sturm chain analysis).
So if you are dealing with a succession of
(quartic) polynomials, I'll think you'll
have to explain how the polynomials depend
on [t0,t1] if you hope for a more targeted
suggestion about computing.
regards, chip
pnachtwey
should be.
(%i5) eq0: x0+v0*t+(a0/2)*t^2+c3*t^3+c4*t^4+c5*t^5-x1$ /* Position,
this is just a starting point */
eq1: diff(eq0,t,1); /* Velocity, find
the zeros between 0 and t1 */
eq2: diff(eq1,t,1); /* Acceleration */
eq4: -remainder(eq2,eq3,t);
(%o6) v0+5*c5*t^4+4*c4*t^3+3*c3*t^2+a0*t
(%o7) 20*c5*t^3+12*c4*t^2+6*c3*t+a0
(%o8) -(20*c5*v0+(30*c3*c5-12*c4^2)*t^2+(15*a0*c5-6*c3*c4)*t-a0*c4)/
(20*c5)
(%o9) (t*((1000*c3*c5^3-400*c4^2*c5^2)
*v0-375*a0^2*c5^3+(700*a0*c3*c4-450*c3^3)
*c5^2+(120*c3^2*c4^2-160*a0*c4^3)*c5)+
(-500*a0*c5^3+800*c3*c4*c5^2-240*c4^3*c5)*v0+(25*a0^2*c4-75*a0*c3^2)
*c5^2+20*a0*c3*c4^2*c5)/(75*c3^2*c5^2-60*c3*c4^2*c5+12*c4^4)
My error was not specifying the variable being divided and wxMaxima is
confused. This also explains why I was getting different symbolic and
numeric results. Notice that I specified a third parameter in my
remainder function. Thanks for being the doubter.
The Sturm chain at t1 is still [0;0;0;0;0]
However, the Sturm chain at t0 is now more complicated because the 4th
term needs to be calculated and is messy. I am going to investigate
other suggestions. Its time to have another German lesson.
Peter Nachtwey
JCH
----- Original Message -----
From: "pnachtwey" <pnac...@gmail.com>
Newsgroups: sci.math,sci.math.num-analysis,sci.math.symbolic
Sent: Tuesday, February 10, 2009 7:20 AM
Subject: Re: How to find a if there is a zero or root between two points quickly
[...]
However, the Sturm chain at t0 is now more complicated because the 4th
term needs to be calculated and is messy. I am going to investigate
other suggestions.
*** JCH
-3 + 5 * x_1 ^ 4 = 0 y' (x_1) = 0
20 * x_2 ^ 3 - 5 = 0 y''(x_2) = 5
x_1 = -0,880111737
x_2 = 0,629960525
c_1 = -3
c_2 = 0
c_3 = 0
c_4 = 0
c_6 = 0
c_7 = 0
c_8 = 0
c_9 = 0
b = 0
b = 0
b = 0
b = 1,32835511
b = -1,32835511
Bart Goddard
be75-18d...@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.
x-x0 and x-x1, leaving a quadratic, which is easily solved
for the other two roots.
Peter Spellucci
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 (schw...@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
* <math.h>, 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 <math.h>
#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