Rskm Posted February 9 Posted February 9 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 Posted February 10 Posted February 10 Maybe this is what you are looking for... Some guy's script + some other guy's script = my script!
AspirinJunkie Posted February 10 Posted February 10 (edited) 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 February 10 by AspirinJunkie
Rskm Posted February 10 Author Posted February 10 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 Posted February 10 Posted February 10 (edited) 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 February 10 by AspirinJunkie I have now also pushed the _lp_ggev() function to GitHub. I have therefore removed it from the example here.
Rskm Posted February 16 Author Posted February 16 Hi, I tried the above two methods and I get error that the Matrix B is not invertible. I am deriving matrices 'Kgredm' and 'Mgredm' - both 60x60 and symmetric and invertible, based on some calculations and assigning them to $mA and $mB respectively. These matrices are invertible when I try them in Mathcad and the eigen values are what I was expecting. However, the autoit script is not proceeding beyond line 456. Kindly give me your opinion.
Rskm Posted February 16 Author Posted February 16 (edited) 2 hours ago, Rskm said: Hi, I tried the above two methods and I get error that the Matrix B is not invertible. I am deriving matrices 'Kgredm' and 'Mgredm' - both 60x60 and symmetric and invertible, based on some calculations and assigning them to $mA and $mB respectively. These matrices are invertible when I try them in Mathcad and the eigen values are what I was expecting. However, the autoit script is not proceeding beyond line 456. Kindly give me your opinion. I am attaching the Kregdm and Mgredm matrices in excel format. I actually tried with simple 3 x 3 matrices and I get the same error for that too. K-lambdaM=0.xlsx Edited February 16 by Rskm
AspirinJunkie Posted February 17 Posted February 17 Either your problem lies in reading the matrices. Because it works wonderfully for me with your matrices: expandcollapse popup#include "LinearAlgebra.au3"; convert AutoIt arrays to UDF Matrix type Global $mK = _la_fromArray($aRawK), _ $mM = _la_fromArray($aRawM) ; display the matrices _la_display($mK, "K") _la_display($mM, "M") ; 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($mK, $mM, "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 ; display the eigenvectors (here only the first vector) Local $v0 = _la_getColumn($mResult.VR, 0) _la_display($v0, "first Eigenvector") Or - as I suspect - you don't have the corresponding OpenBLAS DLL file stored in the path (see the UDF description here in the forum or on GitHub).
Recommended Posts
Create an account or sign in to comment
You need to be a member in order to leave a comment
Create an account
Sign up for a new account in our community. It's easy!
Register a new accountSign in
Already have an account? Sign in here.
Sign In Now