VBA Matrix Operations Help Needed.

**dennygibson** · 02-13-2019, 02:29 PM

Hi … here's a snippet from a macro I'm working on … I can't get the last statement (multiply a matrix by the inverse of another matrix) Everything else works but the last statement:

Option Explicit
Option Base 1
Sub Coefficients()
Dim C() As Variant
Dim D() As Variant
Dim F() As Variant 
ReDim C(1 To 2, 1 To 2)
ReDim D(1 To 2, 1 To 1)
ReDim F(1 To 2, 1 To 1)
ReDim Cinv(1 To 2, 1 To 2)
'Set initial conditions
        W1 = 0
        Del1 = 0.05
'    Swirl Chamber Coefficient Assembly
        C(1, 1) = Del1 * W1
        C(1, 2) = Del1 ^ 2
        C(2, 1) = ((A2 - B2) * Del1 * W1 ^ 2)
        C(2, 2) = ((A2 - 2 * B2 + 1) * Del1 ^ 2 * W1)
        D(1, 1) = B1 / A1
        D(2, 1) = C2 * W1
        ActiveSheet.Range("T5") = C(1, 1)
        ActiveSheet.Range("U5") = C(1, 2)
        ActiveSheet.Range("T6") = C(2, 1)
        ActiveSheet.Range("U6") = C(2, 2)
        ActiveSheet.Range("X5") = D(1, 1)
        ActiveSheet.Range("X6") = D(2, 1)
'    Get F
'    Note: we have {C}dx/dt={D} ... need: dx/dt={F}=[C_inverse]{D}
        F = WorksheetFunction.MMult(WorksheetFunction.MInverse(C), D) --- this is the statement that bombs!!!
End Sub

**Paul_Hossler** · 02-13-2019, 02:54 PM

Welcome to the forum

Take a minute to read the links and the FAQ in my sigature, especiallt the part about using code tags and including a small sample workbook with the macros and some data that shows the issue or desired result

As a GUESS …

Option Explicit
Option Base 1
Sub Coefficients()
Dim C() As Variant
Dim D() As Variant
Dim F As Variant '--------------------------------------------- make F a Variant, not an array of Variants
ReDim C(1 To 2, 1 To 2)
ReDim D(1 To 2, 1 To 1)
' ReDim F(1 To 2, 1 To 1) ------------------------------------------------------------------- no need to ReDim since above change
ReDim Cinv(1 To 2, 1 To 2)
'Set initial conditions
        W1 = 0
        Del1 = 0.05
'    Swirl Chamber Coefficient Assembly
        C(1, 1) = Del1 * W1
        C(1, 2) = Del1 ^ 2
        C(2, 1) = ((A2 - B2) * Del1 * W1 ^ 2)
        C(2, 2) = ((A2 - 2 * B2 + 1) * Del1 ^ 2 * W1)
        D(1, 1) = B1 / A1
        D(2, 1) = C2 * W1
        ActiveSheet.Range("T5") = C(1, 1)
        ActiveSheet.Range("U5") = C(1, 2)
        ActiveSheet.Range("T6") = C(2, 1)
        ActiveSheet.Range("U6") = C(2, 2)
        ActiveSheet.Range("X5") = D(1, 1)
        ActiveSheet.Range("X6") = D(2, 1)
'    Get F
'    Note: we have {C}dx/dt={D} ... need: dx/dt={F}=[C_inverse]{D}
        F = WorksheetFunction.MMult(WorksheetFunction.MInverse(C), D) --- this is the statement that bombs!!!
End Sub

**MagPower** · 02-13-2019, 03:19 PM

How are A2, B2, C2, defined, and what are their values?

**dennygibson** · 02-13-2019, 04:46 PM

Originally Posted by MagPower

How are A2, B2, C2, defined, and what are their values?

I knew someone would ask … I didn't include them because the whole code is fairly long (and getting longer). Here's the hole thing so far:

Option Explicit
Option Base 1
Sub Coefficients()
Dim ShMain As Worksheet
Dim Rswrl As Variant
Dim Lswrl As Variant
Dim Ror As Variant
Dim Alpha As Variant
Dim rho_liq As Variant
Dim visc As Variant
Dim rho_air As Variant
Dim surfT_liquid As Variant
Dim Omega As Variant
Dim Re As Variant
Dim Ain As Variant
Dim V1 As Variant
Dim Ri As Variant
Dim R As Variant
Dim psi As Variant
Dim psipsi As Variant
Dim a As Variant
Dim b As Variant
Dim h As Variant
Dim sum As Variant
Dim First As Variant
Dim Last As Variant
Dim A1 As Variant
Dim A2 As Variant
Dim B1 As Variant
Dim B2 As Variant
Dim C2 As Variant
Dim A3 As Variant
Dim B3 As Variant
Dim B33 As Variant
Dim C3 As Variant
Dim B4 As Variant
Dim A4 As Variant
Dim W1 As Variant
Dim Del1 As Variant
Dim C As Variant
Dim D As Variant
Dim F As Variant
Dim Pi As Variant
ReDim C(1 To 2, 1 To 2)
ReDim D(1 To 2, 1 To 1)
ReDim F(1 To 2, 1 To 1)
ReDim Cinv(1 To 2, 1 To 2)
Dim npts As Long
Dim k As Long
' Input Values:
' Geometry:
    Rswrl = ActiveSheet.Range("H18").Value / 2
    Lswrl = ActiveSheet.Range("H14").Value
    Ror = ActiveSheet.Range("H20").Value / 2
    Alpha = ActiveSheet.Range("H21").Value
    Ri = ActiveSheet.Range("H17").Value
' Properties:
    rho_liq = ActiveSheet.Range("H8").Value
    visc = ActiveSheet.Range("H9").Value
    rho_air = ActiveSheet.Range("H10").Value
    surfT_liquid = ActiveSheet.Range("H11").Value
' Preliminaries
    Pi = 3.14159265
    Ain = ActiveSheet.Range("H23").Value
    V1 = ActiveSheet.Range("H24").Value
' Swirl Chamber
    Omega = V1 * (Rswrl - Ri)
    Re = Omega / visc
    R = Rswrl
'   Coefficients (constants - variable values must be added.
'    Coefficient A1
        npts = 100
        b = 1#
        a = 0#
        h = (b - a) / npts
    
'    First and Last Integration Points
        First = (2 * a - a ^ 2) - (2 * a - a ^ 2) ^ 2
        Last = (2 * b - b ^ 2) - (2 * b - b ^ 2) ^ 2
        sum = First
'    Odd terms
        For k = 1 To npts - 1 Step 2
            psi = 2 * k * h - (k * h) * (k * h)
            sum = sum + 4 * (psi - psi * psi)
        Next k
'    Even terms
        For k = 2 To npts - 1 Step 2
            psi = 2 * k * h - (k * h) * (k * h)
            sum = sum + 2 * (psi - psi * psi)
        Next k
 '   Add final value
        sum = sum + Last
        A1 = h / 3 * sum
        ActiveSheet.Range("R5") = A1
    
'    Coefficient A2
        npts = 50
        b = 1#
        a = 0#
        h = (b - a) / npts
'    First and Last Integration Points
        First = (2 * a - a ^ 2)
        Last = (2 * b - b ^ 2)
        sum = First
'    Odd terms
        For k = 1 To npts - 1 Step 2
            psi = 2 * k * h - (k * h) * (k * h)
            sum = sum + 4 * psi
        Next k
'    Even terms
        For k = 2 To npts - 1 Step 2
            psi = 2 * k * h - (k * h) * (k * h)
            sum = sum + 2 * psi
        Next k
 '   Add final value
        sum = sum + Last
        A2 = h / 3 * sum
        ActiveSheet.Range("R6") = A2
 
 '   Coefficient B1 and C2
        b = 1#
        a = 0#
        B1 = (2 - 2 * b) - (2 - 2 * a)
        ActiveSheet.Range("R7") = B1
        C2 = B1
        ActiveSheet.Range("R9") = C2
 
 '   Coefficient B2
        npts = 100
        b = 1#
        a = 0#
        h = (b - a) / npts
'    First and Last Integration Points
        First = (2 * a - a ^ 2) ^ 2
        Last = (2 * b - b ^ 2) ^ 2
        sum = First
'    Odd terms
        For k = 1 To npts - 1 Step 2
            psipsi = (2 * k * h - (k * h) * (k * h)) * (2 * k * h - (k * h) * (k * h))
            sum = sum + 4 * psipsi
        Next k
'    Even terms
        For k = 2 To npts - 1 Step 2
            psipsi = (2 * k * h - (k * h) * (k * h)) * (2 * k * h - (k * h) * (k * h))
            sum = sum + 2 * psipsi
        Next k
 '   Add final value
        sum = sum + Last
        B2 = h / 3 * sum
        ActiveSheet.Range("R8") = B2
'Set initial conditions
        W1 = 0
        Del1 = 0.05
'    Swirl Chamber Coefficient Assembly
        C(1, 1) = Del1 * W1
        C(1, 2) = Del1 ^ 2
        C(2, 1) = ((A2 - B2) * Del1 * W1 ^ 2)
        C(2, 2) = ((A2 - 2 * B2 + 1) * Del1 ^ 2 * W1)
        D(1, 1) = B1 / A1
        D(2, 1) = C2 * W1
        ActiveSheet.Range("T5") = C(1, 1)
        ActiveSheet.Range("U5") = C(1, 2)
        ActiveSheet.Range("T6") = C(2, 1)
        ActiveSheet.Range("U6") = C(2, 2)
        ActiveSheet.Range("X5") = D(1, 1)
        ActiveSheet.Range("X6") = D(2, 1)
'    Get F
        F = WorksheetFunction.MMult(WorksheetFunction.MInverse(C), D)
        
        
' Spin Chamber
'    Coefficients
        A3 = A1
        B3 = A2
        B33 = B2
        C3 = B1
        B4 = C2
        A4 = A1
        ActiveSheet.Range("R10") = A3
        ActiveSheet.Range("R11") = B3
        ActiveSheet.Range("R12") = B33
        ActiveSheet.Range("R13") = C3
        ActiveSheet.Range("R14") = B4
        ActiveSheet.Range("R15") = A4

' Spray Orifice
        Omega = ActiveSheet.Range("H25")
        Re = ActiveSheet.Range("H26")
        R = Ror
    
' Solve Momentum Equations

        
' Zone 1
        
End Sub

**dennygibson** · 02-13-2019, 04:59 PM

Paul … thanks! I read it. Ok … I'm not familiar with [/CODE] tags rather than ask you what they do I'll go add to my code and see for myself. Also … I'm a bit confused … F is a 2x1 matrix.

**Paul_Hossler** · 02-13-2019, 05:27 PM

1. You can use the [#] icon to insert CODE and /CODE tags -- paste the macro between then

2. The reason I suggested attaching a workbook is you're getting values from the worksheet

    Rswrl = ActiveSheet.Range("H18").Value / 2
    Lswrl = ActiveSheet.Range("H14").Value
    Ror = ActiveSheet.Range("H20").Value / 2
    Alpha = ActiveSheet.Range("H21").Value
    Ri = ActiveSheet.Range("H17").Value
    
    ' Properties:
    rho_liq = ActiveSheet.Range("H8").Value
    visc = ActiveSheet.Range("H9").Value
    rho_air = ActiveSheet.Range("H10").Value
    surfT_liquid = ActiveSheet.Range("H11").Value

and I don't have any idea what a reasonable guess-timent is. Bottom right, [Go Advanced] and use the paperclip icon

3. When I ran the macro with my guesses, I go the error on MInverse, and a Watch window shows.

Capture.JPG

I'm fuzzy, but I don't think C is invertable (key word is 'fuzzy')

**dennygibson** · 02-13-2019, 05:32 PM

Paul … when I pasted my code, highlighted it, and pressed # … the code disappeared and left only # in its place. Also … your suggestion didn't fix it. Here's my code again with better formatting:

Dim visc As Variant
Dim rho_air As Variant
Dim surfT_liquid As Variant
Dim Omega As Variant
Dim Re As Variant
Dim Ain As Variant
Dim V1 As Variant
Dim Ri As Variant
Dim R As Variant
Dim psi As Variant
Dim psipsi As Variant
Dim a As Variant
Dim b As Variant
Dim h As Variant
Dim sum As Variant
Dim First As Variant
Dim Last As Variant
Dim A1 As Variant
Dim A2 As Variant
Dim B1 As Variant
Dim B2 As Variant
Dim C2 As Variant
Dim A3 As Variant
Dim B3 As Variant
Dim B33 As Variant
Dim C3 As Variant
Dim B4 As Variant
Dim A4 As Variant
Dim W1 As Variant
Dim Del1 As Variant
Dim C()As Variant
Dim D() As Variant
Dim F As Variant
Dim Pi As Variant
ReDim C(1 To 2, 1 To 2)
ReDim D(1 To 2, 1 To 1)
'ReDim F(1 To 2, 1 To 1)
Dim npts As Long
Dim k As Long
' Input Values:
' Geometry:
Rswrl = ActiveSheet.Range("H18").Value / 2
Lswrl = ActiveSheet.Range("H14").Value
Ror = ActiveSheet.Range("H20").Value / 2
Alpha = ActiveSheet.Range("H21").Value
Ri = ActiveSheet.Range("H17").Value
' Properties:
rho_liq = ActiveSheet.Range("H8").Value
visc = ActiveSheet.Range("H9").Value
rho_air = ActiveSheet.Range("H10").Value
surfT_liquid = ActiveSheet.Range("H11").Value
' Preliminaries
Pi = 3.14159265
Ain = ActiveSheet.Range("H23").Value
V1 = ActiveSheet.Range("H24").Value
' Swirl Chamber
Omega = V1 * (Rswrl - Ri)
Re = Omega / visc
R = Rswrl
' Coefficients (constants - variable values must be added.
' Coefficient A1
npts = 100
b = 1#
a = 0#
h = (b - a) / npts

' First and Last Integration Points
First = (2 * a - a ^ 2) - (2 * a - a ^ 2) ^ 2
Last = (2 * b - b ^ 2) - (2 * b - b ^ 2) ^ 2
sum = First
' Odd terms
For k = 1 To npts - 1 Step 2
psi = 2 * k * h - (k * h) * (k * h)
sum = sum + 4 * (psi - psi * psi)
Next k
' Even terms
For k = 2 To npts - 1 Step 2
psi = 2 * k * h - (k * h) * (k * h)
sum = sum + 2 * (psi - psi * psi)
Next k
' Add final value
sum = sum + Last
A1 = h / 3 * sum
ActiveSheet.Range("R5") = A1

' Coefficient A2
npts = 50
b = 1#
a = 0#
h = (b - a) / npts
' First and Last Integration Points
First = (2 * a - a ^ 2)
Last = (2 * b - b ^ 2)
sum = First
' Odd terms
For k = 1 To npts - 1 Step 2
psi = 2 * k * h - (k * h) * (k * h)
sum = sum + 4 * psi
Next k
' Even terms
For k = 2 To npts - 1 Step 2
psi = 2 * k * h - (k * h) * (k * h)
sum = sum + 2 * psi
Next k
' Add final value
sum = sum + Last
A2 = h / 3 * sum
ActiveSheet.Range("R6") = A2

' Coefficient B1 and C2
b = 1#
a = 0#
B1 = (2 - 2 * b) - (2 - 2 * a)
ActiveSheet.Range("R7") = B1
C2 = B1
ActiveSheet.Range("R9") = C2

' Coefficient B2
npts = 100
b = 1#
a = 0#
h = (b - a) / npts
' First and Last Integration Points
First = (2 * a - a ^ 2) ^ 2
Last = (2 * b - b ^ 2) ^ 2
sum = First
' Odd terms
For k = 1 To npts - 1 Step 2
psipsi = (2 * k * h - (k * h) * (k * h)) * (2 * k * h - (k * h) * (k * h))
sum = sum + 4 * psipsi
Next k
' Even terms
For k = 2 To npts - 1 Step 2
psipsi = (2 * k * h - (k * h) * (k * h)) * (2 * k * h - (k * h) * (k * h))
sum = sum + 2 * psipsi
Next k
' Add final value
sum = sum + Last
B2 = h / 3 * sum
ActiveSheet.Range("R8") = B2
'Set initial conditions
W1 = 0
Del1 = 0.05
' Swirl Chamber Coefficient Assembly
C(1, 1) = Del1 * W1
C(1, 2) = Del1 ^ 2
C(2, 1) = ((A2 - B2) * Del1 * W1 ^ 2)
C(2, 2) = ((A2 - 2 * B2 + 1) * Del1 ^ 2 * W1)
D(1, 1) = B1 / A1
D(2, 1) = C2 * W1
ActiveSheet.Range("T5") = C(1, 1)
ActiveSheet.Range("U5") = C(1, 2)
ActiveSheet.Range("T6") = C(2, 1)
ActiveSheet.Range("U6") = C(2, 2)
ActiveSheet.Range("X5") = D(1, 1)
ActiveSheet.Range("X6") = D(2, 1)
' Get F
F = WorksheetFunction.MMult(WorksheetFunction.MInverse(C), D) 'this is the problem statement.

' Spin Chamber
' Coefficients
A3 = A1
B3 = A2
B33 = B2
C3 = B1
B4 = C2
A4 = A1
ActiveSheet.Range("R10") = A3
ActiveSheet.Range("R11") = B3
ActiveSheet.Range("R12") = B33
ActiveSheet.Range("R13") = C3
ActiveSheet.Range("R14") = B4
ActiveSheet.Range("R15") = A4

' Spray Orifice
Omega = ActiveSheet.Range("H25")
Re = ActiveSheet.Range("H26")
R = Ror

' Solve Momentum Equations

' Zone 1

End Sub

**MagPower** · 02-13-2019, 05:40 PM

Just curious, is this a problem from fluid mechanics?

**dennygibson** · 02-13-2019, 06:20 PM

Paul … yes … it's the beginnings of a theoretical, viscous simulation of a pressure swirl atomizer. Let me see if I can get this right this time ... I attached the workbook … hopefully you can access it. I changed the initial conditions to small numbers so that each cell in the matrix had a value and got the inverse to write to the worksheet. Apparently you can't invert a matrix with zeros in all but one array location (guess it makes sense).

Now I'm on to figuring out Runge Kutta 4th order with matrices … not a clue yet.

**dennygibson** · 02-13-2019, 06:25 PM

and let me try to post the code with tags:

Option Explicit
Option Base 1
Sub Coefficients()
Dim Rswrl As Variant
Dim Lswrl As Variant
Dim Ror As Variant
Dim Alpha As Variant
Dim rho_liq As Variant
Dim visc As Variant
Dim rho_air As Variant
Dim surfT_liquid As Variant
Dim Omega As Variant
Dim Re As Variant
Dim Ain As Variant
Dim V1 As Variant
Dim Ri As Variant
Dim R As Variant
Dim psi As Variant
Dim psipsi As Variant
Dim a As Variant
Dim b As Variant
Dim h As Variant
Dim sum As Variant
Dim First As Variant
Dim Last As Variant
Dim A1 As Variant
Dim A2 As Variant
Dim B1 As Variant
Dim B2 As Variant
Dim C2 As Variant
Dim A3 As Variant
Dim B3 As Variant
Dim B33 As Variant
Dim C3 As Variant
Dim B4 As Variant
Dim A4 As Variant
Dim W1 As Variant
Dim Del1 As Variant
Dim C()As Variant
Dim D() As Variant
Dim F As Variant
Dim Pi As Variant
ReDim C(1 To 2, 1 To 2)
ReDim D(1 To 2, 1 To 1)
'ReDim F(1 To 2, 1 To 1)
Dim npts As Long
Dim k As Long
' Input Values:
' Geometry:
    Rswrl = ActiveSheet.Range("H18").Value / 2
    Lswrl = ActiveSheet.Range("H14").Value
    Ror = ActiveSheet.Range("H20").Value / 2
    Alpha = ActiveSheet.Range("H21").Value
    Ri = ActiveSheet.Range("H17").Value
' Properties:
    rho_liq = ActiveSheet.Range("H8").Value
    visc = ActiveSheet.Range("H9").Value
    rho_air = ActiveSheet.Range("H10").Value
    surfT_liquid = ActiveSheet.Range("H11").Value
' Preliminaries
    Pi = 3.14159265
    Ain = ActiveSheet.Range("H23").Value
    V1 = ActiveSheet.Range("H24").Value
' Swirl Chamber
    Omega = V1 * (Rswrl - Ri)
    Re = Omega / visc
    R = Rswrl
'   Coefficients (constants - variable values must be added.
'    Coefficient A1
        npts = 100
        b = 1#
        a = 0#
        h = (b - a) / npts
    
'    First and Last Integration Points
        First = (2 * a - a ^ 2) - (2 * a - a ^ 2) ^ 2
        Last = (2 * b - b ^ 2) - (2 * b - b ^ 2) ^ 2
        sum = First
'    Odd terms
        For k = 1 To npts - 1 Step 2
            psi = 2 * k * h - (k * h) * (k * h)
            sum = sum + 4 * (psi - psi * psi)
        Next k
'    Even terms
        For k = 2 To npts - 1 Step 2
            psi = 2 * k * h - (k * h) * (k * h)
            sum = sum + 2 * (psi - psi * psi)
        Next k
 '   Add final value
        sum = sum + Last
        A1 = h / 3 * sum
        ActiveSheet.Range("R5") = A1
    
'    Coefficient A2
        npts = 50
        b = 1#
        a = 0#
        h = (b - a) / npts
'    First and Last Integration Points
        First = (2 * a - a ^ 2)
        Last = (2 * b - b ^ 2)
        sum = First
'    Odd terms
        For k = 1 To npts - 1 Step 2
            psi = 2 * k * h - (k * h) * (k * h)
            sum = sum + 4 * psi
        Next k
'    Even terms
        For k = 2 To npts - 1 Step 2
            psi = 2 * k * h - (k * h) * (k * h)
            sum = sum + 2 * psi
        Next k
 '   Add final value
        sum = sum + Last
        A2 = h / 3 * sum
        ActiveSheet.Range("R6") = A2
 
 '   Coefficient B1 and C2
        b = 1#
        a = 0#
        B1 = (2 - 2 * b) - (2 - 2 * a)
        ActiveSheet.Range("R7") = B1
        C2 = B1
        ActiveSheet.Range("R9") = C2
 
 '   Coefficient B2
        npts = 100
        b = 1#
        a = 0#
        h = (b - a) / npts
'    First and Last Integration Points
        First = (2 * a - a ^ 2) ^ 2
        Last = (2 * b - b ^ 2) ^ 2
        sum = First
'    Odd terms
        For k = 1 To npts - 1 Step 2
            psipsi = (2 * k * h - (k * h) * (k * h)) * (2 * k * h - (k * h) * (k * h))
            sum = sum + 4 * psipsi
        Next k
'    Even terms
        For k = 2 To npts - 1 Step 2
            psipsi = (2 * k * h - (k * h) * (k * h)) * (2 * k * h - (k * h) * (k * h))
            sum = sum + 2 * psipsi
        Next k
 '   Add final value
        sum = sum + Last
        B2 = h / 3 * sum
        ActiveSheet.Range("R8") = B2
'Set initial conditions
        W1 = 0
        Del1 = 0.05
'    Swirl Chamber Coefficient Assembly
        C(1, 1) = Del1 * W1
        C(1, 2) = Del1 ^ 2
        C(2, 1) = ((A2 - B2) * Del1 * W1 ^ 2)
        C(2, 2) = ((A2 - 2 * B2 + 1) * Del1 ^ 2 * W1)
        D(1, 1) = B1 / A1
        D(2, 1) = C2 * W1
        ActiveSheet.Range("T5") = C(1, 1)
        ActiveSheet.Range("U5") = C(1, 2)
        ActiveSheet.Range("T6") = C(2, 1)
        ActiveSheet.Range("U6") = C(2, 2)
        ActiveSheet.Range("X5") = D(1, 1)
        ActiveSheet.Range("X6") = D(2, 1)
'    Get F
        F = WorksheetFunction.MMult(WorksheetFunction.MInverse(C), D)
        
        
' Spin Chamber
'    Coefficients
        A3 = A1
        B3 = A2
        B33 = B2
        C3 = B1
        B4 = C2
        A4 = A1
        ActiveSheet.Range("R10") = A3
        ActiveSheet.Range("R11") = B3
        ActiveSheet.Range("R12") = B33
        ActiveSheet.Range("R13") = C3
        ActiveSheet.Range("R14") = B4
        ActiveSheet.Range("R15") = A4

' Spray Orifice
        Omega = ActiveSheet.Range("H25")
        Re = ActiveSheet.Range("H26")
        R = Ror
    
' Solve Momentum Equations

        
' Zone 1
        
End Sub

**dennygibson** · 02-13-2019, 06:29 PM

MagPower ... sorry ... I got your reply mixed up with Paul's ... yes, this is a fluid mechnics problem ... trying to simulate swirling flow in a simplex atomizer.

**Paul_Hossler** · 02-13-2019, 07:12 PM

@dennygibson --

+1 on the CODE tags

+1 on attaching WB with 'real' data

. I didn't understand it, but I'll just go with the flow

+1 on running to completion with no debug messages

Few minor suggestions about declaring variables you might consider

a. Since they're numbers, a Double is more efficient that a Variant
b. Sum is a Excel keyword, so a different one might be less confusing
c .You can Dim arrays with dimensions if the size isn't going to change
d. You can use Excel's Pi


Dim Rswrl As Double '   Double is more efficient than Variants if numbers
Dim Lswrl As Double
Dim Ror As Double
Dim Alpha As Double
Dim rho_liq As Double
Dim visc As Double
Dim rho_air As Double
Dim surfT_liquid As Double
Dim Omega As Double
Dim Re As Double
Dim Ain As Double
Dim V1 As Double
Dim Ri As Double
Dim R As Double
Dim psi As Double
Dim psipsi As Double
Dim a As Double
Dim b As Double
Dim h As Double
Dim sum As Double   '   Sum is keyword - different one would be better
Dim First As Double
Dim Last As Double
Dim A1 As Double
Dim A2 As Double
Dim B1 As Double
Dim B2 As Double
Dim C2 As Double
Dim A3 As Double
Dim B3 As Double
Dim B33 As Double
Dim C3 As Double
Dim B4 As Double
Dim A4 As Double
Dim W1 As Double
Dim Del1 As Double
Dim C(1 To 2, 1 To 2) As Double '   no need to ReDim - just size at beginning
Dim D(1 To 2, 1 To 1) As Double
Dim C_inv As Variant            '   these have to be Variant since you assign to them
Dim F As Variant
Dim Pi As Double
Pi = Application.WorksheetFunction.Pi   '   you can use the built in or make a Const
'Const Pi as double = 3.14159265358979

But what the heck … it works

**dennygibson** · 02-13-2019, 07:45 PM

Paul ... thanks! I appreciate the suggestions. Would this be the place to get some help using runge-kutta 4th order to solve simultaneous equations (using matrix operations) or am I barking up the wrong tree?

**MagPower** · 02-13-2019, 09:01 PM

Originally Posted by dennygibson

Paul ... thanks! I appreciate the suggestions. Would this be the place to get some help using runge-kutta 4th order to solve simultaneous equations (using matrix operations) or am I barking up the wrong tree?

Denny,

I don't know if anyone in this forum has a background using RK to solve a system of DEs. A math forum might be a better place.

-- Russ

**dennygibson** · 02-13-2019, 09:18 PM

Russ ... thanks. I’ll see what I can find.

**MagPower** · 02-13-2019, 09:23 PM

Originally Posted by dennygibson

Russ ... thanks. I’ll see what I can find.

Denny - I found this site. Hopefully this helps:

https://www.researchgate.net/publica...e-Kutta_Method

Russ

**dennygibson** · 02-13-2019, 09:40 PM

Russ ... thanks! I’ll check it out

**dennygibson** · 02-16-2019, 01:45 PM

Hey! Things are going well but I have a variable type mismatch that I haven't been able to figure out. The three statements are highlighted in red. I'm not done compiling but I got hung up there. The workbook is attached here New Formulation.xlsm Take a look and save me again?

Option Explicit
Option Base 1
Sub Coefficients()
Dim Rswrl As Double
Dim Lswrl As Double
Dim Ror As Double
Dim Alpha As Double
Dim rho_liq As Double
Dim visc As Double
Dim rho_air As Double
Dim surfT_liquid As Double
Dim Omega As Double
Dim Re As Double
Dim Ain As Double
Dim V1 As Double
Dim Ri As Double
Dim R As Double
Dim psi As Double
Dim psipsi As Double
Dim a As Double
Dim b As Double
Dim h As Double
Dim sum As Double
Dim First As Double
Dim Last As Double
Dim A1 As Double
Dim A2 As Double
Dim B1 As Double
Dim B2 As Double
Dim C2 As Double
Dim A3 As Double
Dim B3 As Double
Dim B33 As Double
Dim C3 As Double
Dim B4 As Double
Dim A4 As Double
Dim C(1 To 2, 1 To 2) As Double
Dim D(1 To 2, 1 To 1) As Variant
Dim C_inv As Variant
Dim F As Variant
Dim Zswrl As Double
Dim dz1 As Double
Dim K1 As Double
Dim K2 As Double
Dim K3 As Double
Dim K4 As Double
Dim L1 As Double
Dim L2 As Double
Dim L3 As Double
Dim L4 As Double
Dim F1 As Double
Dim F2 As Double
Dim Z1 As Double
Dim Z2 As Double
Dim Z3 As Double
Dim Z4 As Double
Dim Del1 As Variant
Dim Del2 As Double
Dim Del3 As Double
Dim Del4 As Double
Dim W1 As Variant
Dim W2 As Double
Dim W3 As Double
Dim W4 As Double
Dim Del1St As Variant
Dim W1St As Variant
Dim Pi As Double
Dim npts As Long
Dim k As Long
Dim Nswrl As Long
Dim IZ As Long
' Input Values:
' Geometry:
    Rswrl = ActiveSheet.Range("H18").Value / 2
    Lswrl = ActiveSheet.Range("H14").Value
    Ror = ActiveSheet.Range("H20").Value / 2
    Alpha = ActiveSheet.Range("H21").Value
    Ri = ActiveSheet.Range("H17").Value
' Properties:
    rho_liq = ActiveSheet.Range("H8").Value
    visc = ActiveSheet.Range("H9").Value
    rho_air = ActiveSheet.Range("H10").Value
    surfT_liquid = ActiveSheet.Range("H11").Value
' Preliminaries
    Pi = 3.14159265
    Ain = ActiveSheet.Range("H23").Value
    V1 = ActiveSheet.Range("H24").Value
' Swirl Chamber
'   Coefficients (constants - variable values must be added.
'    Coefficient A1
        npts = 20
        b = 1#
        a = 0#
        h = (b - a) / npts
    
'    First and Last Integration Points
        First = (2 * a - a ^ 2) - (2 * a - a ^ 2) ^ 2
        Last = (2 * b - b ^ 2) - (2 * b - b ^ 2) ^ 2
        sum = First
'    Odd terms
        For k = 1 To npts - 1 Step 2
            psi = 2 * k * h - (k * h) * (k * h)
            sum = sum + 4 * (psi - psi * psi)
        Next k
'    Even terms
        For k = 2 To npts - 1 Step 2
            psi = 2 * k * h - (k * h) * (k * h)
            sum = sum + 2 * (psi - psi * psi)
        Next k
 '   Add final value
        sum = sum + Last
        A1 = h / 3 * sum
        ActiveSheet.Range("R5") = A1
    
'    Coefficient A2
        npts = 20
        b = 1#
        a = 0#
        h = (b - a) / npts
'    First and Last Integration Points
        First = (2 * a - a ^ 2)
        Last = (2 * b - b ^ 2)
        sum = First
'    Odd terms
        For k = 1 To npts - 1 Step 2
            psi = 2 * k * h - (k * h) * (k * h)
            sum = sum + 4 * psi
        Next k
'    Even terms
        For k = 2 To npts - 1 Step 2
            psi = 2 * k * h - (k * h) * (k * h)
            sum = sum + 2 * psi
        Next k
 '   Add final value
        sum = sum + Last
        A2 = h / 3 * sum
        ActiveSheet.Range("R6") = A2
 
 '   Coefficient B1 and C2
        b = 1#
        a = 0#
        B1 = (2 - 2 * b) - (2 - 2 * a)
        ActiveSheet.Range("R7") = B1
        C2 = B1
        ActiveSheet.Range("R9") = C2
 
 '   Coefficient B2
        npts = 20
        b = 1#
        a = 0#
        h = (b - a) / npts
'    First and Last Integration Points
        First = (2 * a - a ^ 2) ^ 2
        Last = (2 * b - b ^ 2) ^ 2
        sum = First
'    Odd terms
        For k = 1 To npts - 1 Step 2
            psipsi = (2 * k * h - (k * h) * (k * h)) * (2 * k * h - (k * h) * (k * h))
            sum = sum + 4 * psipsi
        Next k
'    Even terms
        For k = 2 To npts - 1 Step 2
            psipsi = (2 * k * h - (k * h) * (k * h)) * (2 * k * h - (k * h) * (k * h))
            sum = sum + 2 * psipsi
        Next k
 '   Add final value
        sum = sum + Last
        B2 = h / 3 * sum
        ActiveSheet.Range("R8") = B2
        
' Zone 1
'       Circulation Constant
            Omega = V1 * (Rswrl - Ri)
'       Reynolds Number
            Re = Omega / visc
'       Initial value
            W1 = 0.001
            Del1 = 0.001
            Z1 = 0
'       Swirl Chamber Height
            Zswrl = ActiveSheet.Range("H18").Value
'       Number of Steps in Swirl Chamber Height
            Nswrl = 20
'           Delta-Z
            dz1 = Zswrl / Nswrl
        For IZ = 1 To Nswrl
'       Swirl Chamber Coefficient Matrix
            C(1, 1) = Del1 * W1
            C(1, 2) = Del1 ^ 2
            C(2, 1) = ((A2 - B2) * Del1 * W1 ^ 2)
            C(2, 2) = ((A2 - 2 * B2 + 1) * Del1 ^ 2 * W1)
            D(1, 1) = B1 / A1
            D(2, 1) = C2 * W1
        
'   Vector F Initial Conditions
            F = WorksheetFunction.MMult(WorksheetFunction.MInverse(C), D)
            ActiveSheet.Range("Y5") = F(1, 1)
            ActiveSheet.Range("Y6") = F(2, 1)
            F1 = F(1, 1)
            F2 = F(2, 1)
            K1 = dz1 * F1
            L1 = dz1 * F2
            Z1 = Z1 + dz1 / 2
            Del1 = Del1 + K1 / 2
            W1 = W1 + L1 / 2
            C(1, 1) = Del1 * W1
            C(1, 2) = Del1 ^ 2
            C(2, 1) = ((A2 - B2) * Del1 * W1 ^ 2)
            C(2, 2) = ((A2 - 2 * B2 + 1) * Del1 ^ 2 * W1)
            D(1, 1) = B1 / A1
            D(2, 1) = C2 * W1
            F = WorksheetFunction.MMult(WorksheetFunction.MInverse(C), D)
            F1 = F(1, 1)
            F2 = F(2, 1)
            K2 = dz1 * F1
            L2 = dz1 * F2
            Z2 = Z1 + dz1 / 2
            Del2 = Del1 + K2 / 2
            W2 = W1 + L2 / 2
            C(1, 1) = Del2 * W2
            C(1, 2) = Del2 ^ 2
            C(2, 1) = ((A2 - B2) * Del2 * W2 ^ 2)
            C(2, 2) = ((A2 - 2 * B2 + 1) * Del2 ^ 2 * W2)
            D(1, 1) = B1 / A1
            D(2, 1) = C2 * W2
            F = WorksheetFunction.MMult(WorksheetFunction.MInverse(C), D)
            F1 = F(1, 1)
            F2 = F(2, 1)
            K3 = dz1 * F1
            L3 = dz1 * F2
            Z3 = Z1 + dz1 / 2
            Del3 = Del1 + K2 / 2
            W3 = W1 + L2 / 2
            C(1, 1) = Del3 * W3
            C(1, 2) = Del3 ^ 2
            C(2, 1) = ((A2 - B2) * Del3 * W3 ^ 2)
            C(2, 2) = ((A2 - 2 * B2 + 1) * Del3 ^ 2 * W3)
            D(1, 1) = B1 / A1
            D(2, 1) = C2 * W3
            F = WorksheetFunction.MMult(WorksheetFunction.MInverse(C), D)
            F1 = F(1, 1)
            F2 = F(2, 1)
            K4 = dz1 * F1
            L4 = dz1 * F2
            Z3 = Z1 + dz1 / 2
            Del4 = Del1 + K3
            W4 = W1 + L3
            C(1, 1) = Del4 * W4
            C(1, 2) = Del4 ^ 2
            C(2, 1) = ((A2 - B2) * Del4 * W4 ^ 2)
            C(2, 2) = ((A2 - 2 * B2 + 1) * Del4 ^ 2 * W4)
            D(1, 1) = B1 / A1
            D(2, 1) = C2 * W4
            
            Del1 = Del1 + (K1 + 2 * K2 + 2 * K3 + K4) / 6
            W1 = W1 + (L1 + 2 * L2 + 2 * L3 + L4) / 6
            Del1St(IZ) = Del1  'Type Mismatch
            W1St(IZ) = W1     'Type Mismatch
            Z1St(IZ) = Z1       'Type Mismatch
'       Adjust Z1
            Z1 = Z1 + dz1
        Next IZ

' Spin Chamber
'    Coefficients
        A3 = A1
        B3 = A2
        B33 = B2
        C3 = B1
        B4 = C2
        A4 = A1
        ActiveSheet.Range("R10") = A3
        ActiveSheet.Range("R11") = B3
        ActiveSheet.Range("R12") = B33
        ActiveSheet.Range("R13") = C3
        ActiveSheet.Range("R14") = B4
        ActiveSheet.Range("R15") = A4

' Spray Orifice
        Omega = ActiveSheet.Range("H25")
        Re = ActiveSheet.Range("H26")
        R = Ror
    
' Solve Momentum Equations
      
        
End Sub

**Paul_Hossler** · 02-16-2019, 03:50 PM

You Dim-ed Del1St and W1St as Variant, and plain variants don't take an Index (IZ)

I did not see where Z1St was Dim-ed at all

            Del1St(IZ) = Del1  'Type Mismatch
            W1St(IZ) = W1     'Type Mismatch
            Z1St(IZ) = Z1       'Type Mismatch

I'm guessing you wanted

           Dim  Del1St(1 to 20) As Double
           Dim  W11St(1 to 20) As Double
           Dim  Z1St(1 to 20) As Double

**dennygibson** · 02-16-2019, 04:42 PM

Paul ... thanks! As usual you were very helpful!

Thread: VBA Matrix Operations Help Needed.

Thread Tools

Display

VBA Matrix Operations Help Needed.

Next Issue

Posting Permissions