Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.forth > #12717

Re: evplace (eigenvalue placement)

Show all headers | View raw

On Mon, 04 Jun 2012 20:47:47 +0200, Marcel Hendrix wrote:

> The subproblem that caused me some unexpected headache was calculation
> of the polynomial coefficients of p(z) when given the matrix of
> eigenvalues P. (Either real or complex conjugated pairs.)
> 
> The definition of p(z) = (z-p1)*(z-p2) ... (z-pn), with pi the
> eigenvalues stored in the vector P.
> E.g. when P = [0.7784+0i, 0.8055+0.1543i, 0.8055-0.1543i]
> p(z) = z^3 - 2.3894 z^2 + 1.926641 z - 0.523582.

Here's a solution using Higher Order Functions with Lazy Evaluation. 
The trick is to correctly represent the polynomial as an lazy infinite 
list. 

\ Multiply a polynomial by ( Z - c )
: *(z-c) { c } ( poly c -- poly' )
   dup [: c -1 * * ;] map 0 cons  \ multiply by -c and demote the poly. 
   ['] +  zipwith ;                

\ Calculates the characteristic polynomial with N factored terms
: char-poly ( c1 c2..cN N -- poly )
  0 0 ... 1 cons          \ this is just 1,0,0...                  
  swap dup 0 do
        >R swap *(z-c) R>
  loop
  1+ take ;

0.7784
0.8055::0.1543    \ ie, 0.8055 + 0.1543i in my Forth. 
0.8055::-0.1543   
3
 
char-poly .list

1 
{-2.38939999999999996838084825867554172873497009277343750 , 0E-53i} 
{1.926641139999999947150968182540964665451900138091569717566046569985
6271903911419940413907170295715332031250 , 0E-106i} 
{-0.52358199521599997754971805363766210317773486756087243684688815687
00625016519242871603082500336281690305497623505329510273931248320877784863114356994628906250 , 
0E-157i} ok

Cheers,
Arnold

Back to comp.lang.forth | Previous | Next — Previous in thread | Find similar | Unroll thread

Thread

evplace (eigenvalue placement) mhx@iae.nl (Marcel Hendrix) - 2012-06-04 20:47 +0200
  Re: evplace (eigenvalue placement) Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-06-04 21:08 -0700
    Re: evplace (eigenvalue placement) mhx@iae.nl (Marcel Hendrix) - 2012-06-05 20:39 +0200
  Re: evplace (eigenvalue placement) Arnold Doray <invalid@invalid.com> - 2012-06-05 16:15 +0000

csiph-web