While I was preparing this post, I found this C# function based on LU matrix decomposition to solve a system of linear equations. I converted it to VBA Function and decided to share it with you. I found this method simpler and quicker than the matrix inverse method.
This function accepts two ranges as parameters and returns the solutions on the worksheet.
Function in a standard module
Option Explicit
Option Base 1
Public Function SolveLinearEquation(LeftPart As Range, RightPart As Range) As Variant()
Dim L() As Variant, U() As Variant, sum As Double, n As Integer, i As Integer, j As Integer, k As Integer, Result() As Variant
Dim tempresult() As Variant, temp2() As Variant
n = LeftPart.Columns.Count
sum = 0
ReDim L(n, n)
ReDim U(n, n)
ReDim tempresult(n)
ReDim Result(n)
'' LU decomposition
For i = 1 To n
For j = i To n
sum = 0
For k = 1 To i
sum = sum + L(i, k) * U(k, j)
Next k
U(i, j) = LeftPart(i, j) - sum
Next j
For j = i To n
If j = i Then
L(i, i) = 1
Else
sum = 0
For k = 1 To i
sum = sum + L(j, k) * U(k, i)
Next k
L(j, i) = (1 / U(i, i)) * (LeftPart(j, i) - sum)
End If
Next j
Next i
''Solving LY =rightpart
For i = 1 To n
sum = 0
For k = 1 To i
sum = sum + L(i, k) * tempresult(k)
Next k
tempresult(i) = RightPart(i) - sum
Next i
''Back substitution UX=Y
For i = n To 1 Step -1
sum = 0
For k = 1 To n
sum = sum + U(i, k) * Result(k)
Next k
Result(i) = (1 / U(i, i)) * (tempresult(i) - sum)
Next i
'' Return the solutions
SolveLinearEquation = WorksheetFunction.Transpose(Result)
End Function
You can copy and paste the above code in a standard module in your workbook and just call the function in a cell to solve your system of linear equations. Don’t forget to use Ctrl+Shift+Enter to enter your function.
How to use this function
Look at this video.
I hope you will find this custom function useful and if you have any suggestions, please leave a comment or email me at info@civilconstructiontools.com.
Hi Rodrigue,
This matrix decomposition by triangulation comes handy for ranks over 3, although it degrades computation performance as the rank increases (it turns to be O(n³)), and some other solvers as Conjugate gradient are better suited when the rank reaches some hundreds. For information purposes only, for rank 3 (and for rank 2 also), there is a direct function to solve the inverse, as explained here https://ardoris.wordpress.com/2008/07/18/general-formula-for-the-inverse-of-a-3×3-matrix/
Kind regards mate
Here is some code Conjugate gradient from my VBA arsenal:
'CONJUGATE GRADIENT METHOD FOR SOLVING EQUATIONS
Public Sub sConjugateGradient()
Dim a() As Double
Dim B() As Double
Dim X() As Double
'-------------------------------
Dim arrA As Variant
Dim arrB As Variant
Dim row As Long
Dim col As Long
With Application
.Calculation = xlCalculationManual
.ScreenUpdating = False
End With
'arrA = Range(Cells(1, 1), Cells(xxx, xxx)).Value2
'arrB = Range(Cells(1, 1), Cells(xxx, xxx)).Value2
ReDim a(LBound(arrA, 1) To UBound(arrA, 1), LBound(arrA, 2) To UBound(arrA, 2))
ReDim B(LBound(arrB, 1) To UBound(arrB, 1))
For row = LBound(arrA, 1) To UBound(arrA, 1)
For col = LBound(arrA, 2) To UBound(arrA, 2)
a(row, col) = arrA(row, col)
Next col
B(row) = arrB(row, 1)
Next row
'-------------------------------
X() = fConjugateGradient_Solver(a(), B())
With Application
.Calculation = xlCalculationAutomatic
.ScreenUpdating = True
End With
Stop
End Sub
Public Function fConjugateGradient_Solver(ByRef a() As Double, _
ByRef B() As Double, _
Optional ByVal dbTolerance As Double = 0.000001, _
Optional ByVal iter_Max As Long = 10) As Double()
' Fletcher-Reeves algorithm (array should be symmetric)
Dim retValue As VbMsgBoxResult
Dim i As Long
Dim j As Long
Dim n As Long
Dim X() As Double
Dim g() As Double
Dim D() As Double
Dim ad() As Double
Dim iter As Long
Dim dad As Double
Dim c As Double
Dim aL As Double
Dim gg1 As Double
Dim gg2 As Double
Dim GG_Previous As Double
Dim bt As Double
'Dim dbTolerance As Double
dbTolerance = 0.000001
'Ensure matrix is square
If (UBound(a, 1) - LBound(a, 1)) _
(UBound(a, 2) - LBound(a, 2)) Then
retValue = fMsgBox("Array [A] is not square", vbCritical)
Exit Function
End If
' Solutions array
n = UBound(a, 1) - LBound(a, 1) + 1
ReDim X(LBound(a, 1) To UBound(a, 1))
' Create auxiliary arrays
ReDim g(LBound(a, 1) To UBound(a, 1))
ReDim D(LBound(a, 1) To UBound(a, 1))
ReDim ad(LBound(a, 1) To UBound(a, 1))
For i = LBound(a, 1) To UBound(a, 1)
'X(I) = 0 'Initialize array
g(i) = -B(i)
D(i) = B(i)
Next i
gg1 = 0
For i = LBound(a, 1) To UBound(a, 1)
gg1 = gg1 + (g(i) * g(i))
Next i
GG_Previous = gg1
Do While gg1 > dbTolerance
' Convergence analysis
If gg1 <= GG_Previous Then
GG_Previous = gg1
Else
' it's not converging
'GoTo ErrConvergence
End If
iter = iter + 1
dad = 0
For i = LBound(a, 1) To UBound(a, 1)
c = 0
For j = LBound(a, 1) To UBound(a, 1)
c = c + (a(i, j) * D(j))
Next j
ad(i) = c
dad = dad + (c * D(i))
Next i
'If VBA.Abs(DAD) iter_Max Then GoTo ErrSlowConvergence ' it's not converging
Loop
' Free memory
Erase g(), D(), ad() ', A(), B()
fConjugateGradient_Solver = X()
ExitProc:
Exit Function
ErrConvergence:
retValue = MsgBox("Conjugate gradient solver it's not converging, array seems to be bad conditioned" & vbNewLine & _
"[" & gg1 & ">" & GG_Previous & "] on iteration [" & iter & "]", _
vbCritical)
GoTo ExitProc
ErrSlowConvergence:
retValue = MsgBox("Conjugate gradient solver it's converging slowly, array seems to be bad conditioned" & vbNewLine & _
"Convergence = [" & gg1 & "] on iteration [" & iter & "]", _
vbCritical)
GoTo ExitProc
End Function