Eigen value for Characteristic equation

[Rskm]

Hi, I wish to know if I can solve characteristic equation (A-lambda. I)=0 using autoit. The size of the matrix is  72x72. thanks

[Werty]

Maybe this is what you are looking for...

 

[AspirinJunkie]

With >>this UDF<<, this could be implemented as follows, for example:

expandcollapse popup

#include "LinearAlgebra.au3"

; Generate test data: 72x72 Matrix as a normal 2D AutoIt-Array
Global Const $N = 72
Global $aRawData[$N][$N]

; create a symmetric matrix (A = A^T) to ensure real eigenvalues.
Global $fVal
For $i = 0 To $N - 1
    For $j = $i To $N - 1
        $fVal = Random(-10.0, 10.0)
        If $i = $j Then $fVal += $N ; Increase diagonal elements (to ensure numerical stability)
        
        $aRawData[$i][$j] = $fVal
        $aRawData[$j][$i] = $fVal ; Symmetry
    Next
Next
; convert AutoIt array to UDF Matrix type
Global $mA = _la_fromArray($aRawData)

; display the input matrix
_la_display($mA, "input matrix A")

; Calculate the eigenvalues for A (corresponding to the solution of the characteristic equation det(A-λI) = 0).

; _la_eigen parameters:
; Param 1: The Matrix
; Param 2: True = Treat as symmetric (uses optimized algorithm) - _la_eigen() also outputs imaginary eigenvalue parts for non-symmetric matrices if required.
; Param 3: True = Compute Eigenvectors as well (needed for verification)
Global $mResult = _la_eigen($mA, True, True)
If @error Then Exit ConsoleWrite("Error during calculation! @error: " & @error & @CRLF)

; 4. Output Results
; For symmetric matrices, eigenvalues are stored in $mResult.R (Real part).
; $mResult.I (Imaginary part) is not used/empty for symmetric matrices.

; Display the first 10 Eigenvalues
Global $aLambda = _la_toArray($mResult.R)
ConsoleWrite(@CRLF & "--- First 10 Eigenvalues (Lambda) ---" & @CRLF)
For $k = 0 To 9
    ConsoleWrite(StringFormat("Lambda[%d] = %.5g", $k, $aLambda[$k]) & @CRLF)
Next


; Verification
; Definition: A * v = lambda * v  <=>  (A - lambda * I) * v = 0
; We check this for the first eigenvalue (Lambda[0]) and its eigenvector (v[0]).

ConsoleWrite(@CRLF & "--- Verification for Lambda[0] ---" & @CRLF)

Global $fLambda0 = $aLambda[0]    ; First eigenvalue
Global $v0       = $mResult.V[0]  ; Corresponding eigenvector

; Calculate: Residual = (A * v) - (lambda * v)
Global $mRes   = _la_sub(_la_mul($mA, $v0), _la_scale($v0, $fLambda0))

; Calculate the Euclidean norm of the residual vector.
; Ideally, this should be very close to 0 (e.g., 1e-14).
Global $fError = _la_norm($mRes)

ConsoleWrite(StringFormat("Check ||A*v - lambda*v|| = %.5e", $fError) & @CRLF)

If $fError < 1e-9 Then
    ConsoleWrite("Result: SUCCESS (Residual is close to zero)" & @CRLF)
Else
    ConsoleWrite("Result: WARNING (Residual is high)" & @CRLF)
EndIf

 

Edited Tuesday at 08:17 AM by AspirinJunkie

[Rskm]

Hi,  Sorry, m equation is NOT of the form  |A.I - lambda.I|  , instead it is |A - lambda.B|= 0 where  A and B are two independent matrices (72x72).    Is it possible to spit the eigen vectors too ?  thanks 

[AspirinJunkie]

2 hours ago, Rskm said:

m equation is NOT of the form  |A.I - lambda.I|  , instead it is |A - lambda.B|= 0 where  A and B are two independent matrices (72x72).    Is it possible to spit the eigen vectors too ? 

That is actually a rather different task now.
If I am not overlooking anything, it should suffice to solve the equation system B*C=A for C and then determine the eigenvalues/vectors for C. It would then look like this, for example:

expandcollapse popup

#include "LinearAlgebra.au3"

; generate test data (two 72x72 Matrices)
Global Const $N = 72
Global $aRawA[$N][$N], $aRawB[$N][$N]

; create a symmetric matrices A and B to ensure real eigenvalues.
; to make B invertible we set much higher values on it`s diagonal
For $i = 0 To $N - 1
    For $j = $i To $N - 1
        ; Fill A
        Local $fValA = Random(-10.0, 10.0)
        $aRawA[$i][$j] = $fValA
        $aRawA[$j][$i] = $fValA
        
        ; Fill B
        Local $fValB = Random(-5.0, 5.0)
        If $i = $j Then $fValB += ($N * 2)
        $aRawB[$i][$j] = $fValB
        $aRawB[$j][$i] = $fValB
    Next
Next

; convert AutoIt arrays to UDF Matrix type
Global $mA = _la_fromArray($aRawA)
Global $mB = _la_fromArray($aRawB)

; solve B * C = A   to C
Global $mC = _la_solve($mB, $mA)
If @error Then Exit ConsoleWrite("Error: Matrix B is likely singular (not invertible)." & @CRLF)

; _la_eigen parameters:
; Param 1: The Matrix
; Param 2 (False): Treat as non-symmetric (C is not symmetric)
; Param 3 (True):  Compute Eigenvectors
Global $mResult = _la_eigen($mC, False, True)
If @error Then Exit ConsoleWrite("Error during eigen calculation! @error: " & @error & @CRLF)
    
; Check if we have complex eigenvalues (Imaginary part != 0)
Local $fMaxImag = _la_amax($mResult.I)
If Abs($fMaxImag) > 1e-9 Then ConsoleWrite("Warning: Complex eigenvalues detected." & @CRLF)

; print out the first 10 eigenvalues
Global $aLambda = _la_toArray($mResult.R)
ConsoleWrite(@CRLF & "--- First 10 Generalized Eigenvalues (Lambda) ---" & @CRLF)
For $k = 0 To 9
    ConsoleWrite(StringFormat("Lambda[%d] = %.5g", $k, $aLambda[$k]) & @CRLF)
Next

; display the eigenvectors (here only the first vector)
Global $v0 = $mResult.VR[0]
_la_display($v0, "first Eigenvector")

When maximum performance and numerical stability are essential, there is even a low-level LAPACK function for this case: DGGEV (solves the generalized eigenvalue problem).
I had not yet implemented this in LAPACK.au3. I have now done so and will add it to the UDF when the opportunity arises.
This would allow the problem to be solved directly (with higher performance and stability).:

expandcollapse popup

#include "LinearAlgebra.au3"

; generate test data (two 72x72 Matrices)
Global Const $N = 72
Global $aRawA[$N][$N], $aRawB[$N][$N]

; create a symmetric matrices A and B to ensure real eigenvalues.
; to make B invertible we set much higher values on it`s diagonal
For $i = 0 To $N - 1
    For $j = $i To $N - 1
        ; Fill A
        Local $fValA = Random(-10.0, 10.0)
        $aRawA[$i][$j] = $fValA
        $aRawA[$j][$i] = $fValA
        
        ; Fill B
        Local $fValB = Random(-5.0, 5.0)
        If $i = $j Then $fValB += ($N * 2)
        $aRawB[$i][$j] = $fValB
        $aRawB[$j][$i] = $fValB
    Next
Next

; convert AutoIt arrays to UDF Matrix type
Global $mA = _la_fromArray($aRawA), _
       $mB = _la_fromArray($aRawB)

; call the low-level LAPACK function DGGEV for solving the generalized eigenvalues problem
; "N" = No Left Vectors, "V" = Compute Right Vectors
Global $mResult = _lp_ggev($mA, $mB, "N", "V")

If @error Then Exit ConsoleWrite("Error in _lp_ggev: " & @error & @CRLF)

; Process and Output Eigenvalues
; Generalized eigenvalues are calculated as: lambda = (AlphaR + i * AlphaI) / Beta

Global $aAlphaR = _la_toArray($mResult.ALPHAR), _
       $aAlphaI = _la_toArray($mResult.ALPHAI), _
       $aBeta   = _la_toArray($mResult.BETA), _
       $fReal, $fImag, $fBeta, $fLamR, $fLamI, $k

; print out the first 10 eigenvalues
ConsoleWrite(@CRLF & "--- First 10 Generalized Eigenvalues (Lambda) ---" & @CRLF)
For $k = 0 To 9
    $fReal = $aAlphaR[$k]
    $fImag = $aAlphaI[$k]
    $fBeta = $aBeta[$k]
    
    ; Avoid division by zero (infinite eigenvalue)
    If Abs($fBeta) < 1e-15 Then
        ConsoleWrite(StringFormat("Lambda[%d] = Infinite", $k) & @CRLF)
    Else
        $fLamR = $fReal / $fBeta
        $fLamI = $fImag / $fBeta
        
        If Abs($fLamI) < 1e-9 Then ; non imaginary
            ConsoleWrite(StringFormat("Lambda[%d] = %.4g", $k, $fLamR) & @CRLF)
        Else ; imaginary
            ConsoleWrite(StringFormat("Lambda[%d] = %.4g + %.4gi", $k, $fLamR, $fLamI) & @CRLF)
        EndIf
    EndIf
Next

; CAUTION with DGGEV and Complex Eigenvalues:
; If the eigenvalue is real, the eigenvector is in the corresponding column.
; If eigenvalues k and k+1 are a complex conjugate pair, the eigenvectors are stored
; as v(k) = VR(:,k) + i*VR(:,k+1) and v(k+1) = VR(:,k) - i*VR(:,k+1).


; display the eigenvectors (here only the first vector)
Local $v0 = _la_getColumn($mResult.VR, 0)
_la_display($v0, "first Eigenvector")

 

Edited Tuesday at 12:39 PM by AspirinJunkie
I have now also pushed the _lp_ggev() function to GitHub. I have therefore removed it from the example here.