Group Abstract Group Abstract

Message Boards Message Boards

2
|
2.1K Views
|
1 Reply
|
3 Total Likes
View groups...
Share
Share this post:

Suggestion: function for Whittaker's root series formula

Posted 1 year ago
POSTED BY: Raul Prisacariu

Interesting method I guess. It is different from the one I'll show, but at least the similarities of being iterative and geometrically convergent.

The (reasonably) well-known "power method" can be used to find the largest root. It amounts to finding the largest eigenvalue of the characteristic matrix (constructed to have the original polynomial, or its negative, as its characteristic polynomial). One starts with a random vector and iterates newvec = mat.oldvec; newvec=newvec/Norm[newvec] and then takes the ratio of the last iterate vector to the matrix times that vector. Finding the smallest root amounts to using "inverse" iterations, that is, multiplying by the matrix inverse instead of the matrix. (In any practical setting, one would use a solver rather than explicit inversion though.)

Here is an example from the paper.

poly = x^3 - 4*x^2 - 321*x + 20;

We construct the matrix and verify that the characteristic polynomial is as claimed.

mat = {{0, 0, -20}, {1, 0, 321}, {0, 1, 4}};
cp = CharacteristicPolynomial[mat, x]

Out[292]= -20 + 321 x + 4 x^2 - x^3mat = {{0, 0, -20}, {1, 0, 321}, {0, 1, 4}};

For illustrative purposes I construct the inverse directly, then start with a random vector, and do 10 iterations.

invmat = Inverse[mat];
SeedRandom[1234];
x0 = RandomReal[1, 3];
Do[nx = invmat . x0;
  x0 = nx/Norm[nx], {10}];

Dividing as described above shows that we have an eigenvalue as the common ratio 9so it converged nicely).

mat . x0/x0

Out[297]= {0.0622577, 0.0622577, 0.0622577}

Now verify it is the smallest root.

In[300]:= NSolve[poly == 0, x]

Out[300]= {{x -> -16.0623}, {x -> 0.0622577}, {x -> 20.}}
POSTED BY: Daniel Lichtblau
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard