NeillSmith

187 Reputation

4 Badges

19 years, 216 days

MaplePrimes Activity


These are answers submitted by NeillSmith

Thanks for the upper bound on the real positive root, and thanks for the tip on how to preserve formatting in my posts.

 I have been experimenting with numerical solutions to x^4 - a*x -1 = 0 in Matlab. The Matlab ROOTS command returns all four roots in around 0.01 seconds over a very wide range of positive values for a. The Matlab command FZERO using x = 1 + 0.5 * a^(1/3) as its intial guess takes about 0.06 seconds to converge to the one real positive root over a wide range of values for a.

The following iteration

x[i+1] = (a*x[i] + 1)^(1/4)

starting with x[1] = 1 + 0.5 * a^(1/3) converges super fast with 10 iterations taking about 0.0003 seconds. The gain per iteration is not quite one significant digit until machine precision is reached. A semilogy plot of x[i+1]-x[i] vs i is essentially a straight line.

Neill S.

Besides working with the general analytic quartic solution, I also started down the follwing path.

restart;
> assume(a > 0, r1, 'real', r2, 'real', r3, 'real', i3, 'real');
> interface(showassumed = 0);
                                      0
Consider the quartic polynomial
> P1 := x^4-a*x-1 = 0;
                               4             
                              x  - a x - 1 = 0
If a = 0, it becomes
> P2 := subs(a = 0, P1);
                                  4       
                                 x  - 1 = 0
with roots
> solve(P2, x);
                                1, -1, I, -I
                               
For a > 0 it has a positive real root r1 > +1, a negative real root r2 in the range -1/a  <  r2  <  0 and a complex conjugate root pair r3 +/- I*i3.
Building the polynomial from these roots we have

> P3 := (x-r1)*(x-r2)*(x-r3-I*i3)*(x-r3+I*i3) = 0;

            (x - r1) (x - r2) (x - r3 - I i3) (x - r3 + I i3) = 0
           
> P3 := collect(simplify(expand(P3)), [x, r1, r2, r3, i3]);

           4                               3   /                              2      2                 \  2
P3 := x  + (-2 r3 - r2 - r1) x  + \(2 r3 + r2) r1 + r3  + i3  + 2 r2 r3/ x

      //    2                    2\         /    2     2\     \       /  2       2\         
   + \\-r3  - 2 r2 r3 - i3 / r1 + \-r3  - i3 / r2/ x + \r3  + i3 / r2 r1 = 0
  
Equating coefficients between P1 and P3 gives

> coeff(lhs(P1), x, 3) = coeff(lhs(P3), x, 3);

                             0 = -2 r3 - r2 - r1
                            
> coeff(lhs(P1), x, 2) = coeff(lhs(P3), x, 2);

                                                            2     2         
                  0 = (2 r3 + r2)  r1     + r3  + i3  + 2 r2 r3
                 
> coeff(lhs(P1), x, 1) = coeff(lhs(P3), x, 1);

                      /    2                     2\         /    2     2 \  
              -a = \-r3  - 2 r2 r3 - i3  / r1 + \-r3  - i3  / r2
             
> coeff(lhs(P1), x, 0) = coeff(lhs(P3), x, 0);

                                   /   2      2\     
                           -1 = \r3  + i3  / r2 r1
>

I believe my statements above about the four roots of the polynomial x^4 - a * x -1 = 0 if a > 0 is correct, but I can't prove them. Assuming I am correct about the properties of these four roots, I have then created four equations in the four real unknowns r1, r2, r3 and i3. I've been playing with those four equations above but have gotten nowhere.

By the way, how does one paste a Maple worksheet into this Web comment block and preserve the spacing of the output lines? I always have to edit my posts after I paste them in to get brackets and exponents to line back up in their proper positions.

Neill S.

Thank you very much.

 

Neill Smith

Page 1 of 1