Bascom的迷你矩阵代数代码教程 -

标题: Bascom的迷你矩阵代数代码教程 [打印本页]

作者: taotie 时间: 2024-10-29 12:40
标题: Bascom的迷你矩阵代数代码教程
'                                           Bascom的迷你矩阵代数代码
'                                           加法、减法、乘法、倒置
'                                           (Natalius Kiedro，2010年3月）

'我错过了矩阵代数的学习，因为最近我在使用AVR系列（参见我之前关于这个问题的帖子）时遇到了一个需要用到‘矩阵乘法’的问题。

'20年前，我通过编写一个程序来玩转矩阵，这个程序可以在电脑屏幕上绘制并旋转生物分子的3D图像。那个
'时期的另一个项目是关于动力学数据的动态模拟和非线性数据拟合。像3D旋转和参数优化这类事情，在很大
'程度上都依赖于矩阵代数。下面是一个关于矩阵的简要介绍，互联网上有成千上万页的资料可以提供更多信息。

'矩阵入门简介’
'===============
'矩阵是一个数字表。在大多数情况下，数字以类似正方形的形式排列，如下所示：

'          | A(1,1) A(1,2) A(1,3) | | 1 2 3 |
'矩阵A  = | A(2,1) A(2,2) A(2,3) | = | 4 5 6 | 是一个3x3矩阵
'             | A(3,1) A(3,2) A(3,3) | | 7 8 9 |

'矩阵中的数字A(i,j)被称为矩阵的元素。元素A(1,1)的值为1，元素A(2,3)的值为6，以此类推。
'变量i和j是矩阵的索引，i指向行，j指向列，元素就位于该行和该列的交叉点上。

'现在，让我们来看看普通代数。我们有四个每个人都喜欢使用的基本运算：

'加法：C = A + B
'减法：C = A - B
'乘法：C = A * B
'除法：C = A / B

'它们在矩阵代数中的等价运算分别被称为矩阵加法、矩阵减法、矩阵乘法，以及——哎呀——没有所谓的矩阵除法。
'在矩阵代数中，除法是分步进行的。让我们用普通代数来演示一下，看看它在矩阵代数中是如何进行的：

'C = 1 / B：在矩阵代数中，我们称之为矩阵求逆。
'C = C * A：这是矩阵乘法。
'因此，C = 1 / B * A = A / B（在矩阵代数中，这种表示并不严谨，仅用于说明概念）
'所以，在矩阵代数中，除法并不是一个基本运算，因为它可以分解为矩阵求逆和矩阵乘法。

'在普通代数中，数字1是特殊的，因为：

'C = C * 1：（当然，纳塔利乌斯，我们都知道这一点！）
'C = C / 1：（拜托，纳塔利乌斯，你想让我们觉得无聊吗，巴斯科默？）

'不，我不想让任何人觉得无聊。我只是想和你分享矩阵的乐趣！就像尼奥几年前做的那样。
'让我们重新加载这个矩阵，它是……的等价物（此处原文被截断，但根据上下文，可以推测“
'它是……的等价物”可能是指矩阵运算在某种情况下的应用或等价表示）。
'数字1的特殊性在于：

'| 1 0 0 |
'| 0 1 0 | 它被称为单位矩阵，也有人称之为
'| 0 0 1 | 统一矩阵——当然，还有像我这样的人喜欢称之为单位矩阵。

'从矩阵的1,1,1路径穿过的线被称为矩阵的对角线。元素0被称为非对角元素。单位矩阵的特别之处在于，
'其对角线上的元素全部为1，而非对角线上的元素全部为0。

'为什么数字1和单位矩阵如此特殊呢？

'想象一下，如果在Bascom中无法使用乘法和除法，而你需要从零开始，在汇编语言中自己编写这些运算的代码。
'你会如何证明你的算法是有效的呢？

'最简单的方法就是将你的输入数字乘以1或除以1，然后观察数字是否保持不变——正如预期的那样。
'因此，如果我们将这个原理应用到矩阵代数中，我们期望：

'矩阵A乘以单位矩阵 = 矩阵A

'如果矩阵乘法的算法是有效的。如果：

'矩阵A乘以矩阵A的逆矩阵 = 单位矩阵

'那么我们就知道我们的矩阵求逆算法是有效的。很简单，对吧？

'Bascom中的矩阵代数
'========================
'矩阵代数通常使用二维数组进行计算，但Bascom并不了解二维数组。而且它也不需要了解，
'因为汇编语言总是将计算机内存视为线性（一维）排列的。让我们回顾一下本教程开头给出的二维矩阵排列的例子。
'在Bascom中，我们将二维矩阵当作一维矩阵来处理。因此，矩阵A变为：

'矩阵A = |A(1,1) A(1,2) A(1,3) A(2,1) A(2,2) A(2,3) A(3,1) A(3,2) A(3,3)|

'= | 1 2 3 4 5 6 7 8 9 |

'= | A(1) A(2) A(3) A(4) A(5) A(6) A(7) A(8) A(9)|

'为了将针对二维矩阵编写的矩阵代数运算转换为适用于“Bascom化”的一维矩阵的运算，
'我们只需将元素A(i,j)的索引i,j转换为一个指向第ij个元素A(ij)（即A(1)至A(9)）的单索引ij。

'这非常简单：以下是一个只执行此操作的“空”矩阵循环：

'Dim I As Byte, I0 As Byte, J As Byte, Ij As Byte, N As Byte
' N是矩阵的大小。对于3x3矩阵，N为3；对于NxN矩阵，N为相应的值
'N = 3
'For I = 1 TO N          '1, 2, 3
'I0 = I - 1             '0, 1, 2
'I0 = I0 * N             '0, 3, 6（计算当前行之前所有行的元素总数）
'For J = 1 TO N          '例如，当I=2时，I0=3，则J=1,2,3
'Ij = I0 + J             '计算一维数组中的索引位置，即4, 5, 6
'Next J
'Next I

'断线内为AI给出的释义：
'////////////////////////////////////////////////////////////////////////////////////////////////////////////
'在矩阵程序中，变量通常用于表示矩阵的维度、索引位置以及辅助计算。以下是对您提供的代码中各个变量的具体含义的解释：

' N As Byte：
'含义：N代表矩阵的大小，即矩阵的行数和列数。对于NxN矩阵，N的值是相同的。
'示例：对于3x3矩阵，N为3。
'I As Byte：
'含义：I用作矩阵行的索引，通常从1开始，到N结束。
'示例：对于3x3矩阵，I的取值为1, 2, 3。
'I0 As Byte：
'含义：I0是一个辅助变量，用于计算当前行之前所有行的元素总数。它首先通过I - 1得到当前行的相对位置（从0开始计数），
'然后乘以N（矩阵的列数/大小）来得到之前所有行的元素总数。
'示例：对于3x3矩阵，当I=1时，I0=0；当I=2时，I0=3（因为第一行有3个元素）；当I=3时，I0=6（因为前两行共有6个元素）。
'J As Byte：
'含义：J用作矩阵列的索引，通常从1开始，到N结束。
'示例：对于3x3矩阵，当I=2（即处理第二行）时，J的取值为1, 2, 3。
'Ij As Byte：
'含义：Ij是一个辅助变量，用于计算当前元素在一维数组中的索引位置。它通过I0 + J得到，
'其中I0表示当前元素之前所有元素的数量，J表示当前元素在当前行中的位置。
'示例：对于3x3矩阵，当I=2（第二行）且J=1时，Ij=I0+J=3+1=4，表示第二行第一个元素在一维数组中的索引为4。
'这段代码的目的是遍历一个NxN矩阵的所有元素，并通过计算每个元素在一维数组中的索引位置，以便
'可以对这些元素进行操作（如读取、修改等）。这种方法在处理需要将二维数组（如矩阵）转换为一维数组进行处理的场景时非常有用。
'////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

'请保持这种“直观”的变量表示法，如Ij、Ji、Ik等，因为当N=10、I=6、J=3时，I0=50。很容易看出Ij=I0+j，因为53=50+3。

'下面的示例和子程序可以在模拟器中运行。我没有配置串行输出，因此您无法直接看到输出结果。

'祝您在矩阵运算中玩得开心！

'作者：纳塔利乌斯

$regfile = "m328pdef.dat"
$crystal = 8000000
$baud = 19200

Const Nbyc = 4
Const Mbyc = 16                                           'Nbyc * Nbyc

Declare Sub Matadd()
Declare Sub Matsub()
Declare Sub Matmul()
Declare Sub Matinv()
Declare Sub Matprint(byval Ksi As Single)

Dim Iby As Byte , I0by As Byte , Ijby As Byte , Ikby As Byte , Irby As Byte
Dim Jby As Byte , J0by As Byte
Dim Kby As Byte , K0by As Byte , Kjby As Byte
Dim Rby As Byte , R0by As Byte , Riby As Byte , Rjby As Byte , Rrby As Byte , Rsby As Byte
Dim Sby As Byte , S0by As Byte , Siby As Byte , Srby As Byte , Ssby As Byte

Dim Iin As Integer
Dim Isi As Single , Jsi As Single , Ksi As Single , Msi As Single
Dim A(mbyc) As Single , B(mbyc) As Single , C(mbyc) As Single
'Temporary:
Dim M1(mbyc) As Single , M2(mbyc) As Single , M3(mbyc) As Single , M4(mbyc) As Single , M5(mbyc ) As Single

'用（伪）随机数填充矩阵A
'
For Iby = 1 To Mbyc
Iin = Rnd(1000)
Isi = Iin
A(iby) = Isi / 100
Next Iby

Print "输入矩阵A"
Matprint A(1)
Matinv
Print "倒置矩阵B"
Matprint B(1)
Matmul
Print "将矩阵A与其逆相乘"
Matprint C(1)
Print "如果算法是正确，那么它就是单位矩阵。"
'请注意，当从矩阵及其逆矩阵重新建立单位矩阵时，浮点运算的有限精度永远不会导致真正的零出现。
'零会表现为极小的正数或负数。这就是为什么打印出的是-0.00和0.00的原因。
End 'end program

'*******************************************************************************
'******************************* 矩阵加法 *******************************
'*******************************************************************************
Sub Matadd                                                 'C = A + B
For Iby = 1 To Mbyc
      C(iby) = A(iby) + B(iby)
Next Iby
End Sub

'*******************************************************************************
'******************************* 矩阵减法 ****************************
'*******************************************************************************
Sub Matsub                                                 'C = A - B
For Iby = 1 To Mbyc
      C(iby) = A(iby) - B(iby)
Next Iby
End Sub

'*******************************************************************************
'*************************** 矩阵乘法*****************************
'*******************************************************************************
Sub Matmul                                                 'C = A * B
For Iby = 1 To Nbyc
      I0by = Iby - 1
      I0by = I0by * Nbyc
   For Jby = 1 To Nbyc
         Ijby = I0by + Jby
         C(ijby) = 0
         For Kby = 1 To Nbyc
            K0by = Kby - 1
            K0by = K0by * Nbyc
            Kjby = K0by + Jby
            Ikby = I0by + Kby
            Jsi = A(ikby) * B(kjby)
            C(ijby) = C(ijby) + Jsi
         Next Kby
   Next Jby
Next Iby
End Sub

'*******************************************************************************
'******************************* 矩阵求逆 ******************************
'*******************************************************************************
Sub Matinv                                                 'B = 1/A
'
'将原始矩阵A复制到M1中
'
For Iby = 1 To Mbyc
      M1(iby) = A(iby)
Next Iby
'
'填充两个单位矩阵：对于i等于j的情况，M2(i,j) = 1，M4(i,j) = 1；
'对于i不等于j的情况，M2(i,j) = 0，M4(i,j) = 0。
For Iby = 1 To Nbyc
      I0by = Iby - 1
      I0by = I0by * Nbyc
   For Jby = 1 To Nbyc
         Ijby = I0by + Jby
         If Iby = Jby Then
            M2(ijby) = 1
            M4(ijby) = 1
         Else
            M2(ijby) = 0
            M4(ijby) = 0
         End If
   Next Jby
Next Iby

For Rby = 1 To Nbyc                   'Row byte  RBy
      R0by = Rby - 1                      '0, 1, 2, 3..
      R0by = R0by * Nbyc                '0,Nbyc,2Nbyc,3Nbyc
      Isi = 0
      '------------------------------------'Pivotation中心点或轴
   For Iby = Rby To Nbyc             '矩阵或表格的对角线开始提取或构建行
         I0by = Iby - 1
         I0by = I0by * Nbyc
         Irby = I0by + Rby             'index of element in 1D: Ir一维中元素的索引
         Jsi = M1(irby)
         Jsi = Abs(jsi)
         If Jsi >= Isi Then             '找到最大的元素
            Isi = Jsi                   '
            Sby = Iby                   '记住用于交换行的（操作或步骤）
         End If
   Next Iby
      S0by = Sby - 1
      S0by = S0by * Nbyc
   For Iby = 1 To Nbyc
         Riby = R0by + Iby             'Row R , column I “R”代表行号，而“I”代表列号。
         Siby = S0by + Iby             'Row to swap,       “要交换的行”
         Swap M1(riby) , M1(siby)
   Next Iby
      Rrby = R0by + Rby
      Jsi = M1(rrby)
      Jsi = Abs(jsi)
   If Jsi < 10e-30 Then
      Print "奇异矩阵：无法进行求逆"
End If
      'End Sub

      Ssby = S0by + Sby
      Rsby = R0by + Sby
      Srby = S0by + Rby
      M4(rrby) = 0
      M4(ssby) = 0
      M4(rsby) = 1
      M4(srby) = 1
      '-----------------------------------' 高斯-约旦消元法（Gauss-Jordan Elimination）
   For Iby = 1 To Nbyc
         I0by = Iby - 1
         I0by = I0by * Nbyc
         Irby = I0by + Rby
         For Jby = 1 To Nbyc
            J0by = Iby - 1
            J0by = J0by * Nbyc
            Ijby = I0by + Jby
            Rjby = R0by + Jby
            If Iby = Rby Then
               If Iby = Jby Then
                  M3(rrby) = 1 / M1(rrby)
               Else
                  M3(ijby) = M1(ijby) / M1(rrby)
                  M3(ijby) = 0 - M3(ijby)
               End If
            Elseif Jby = Rby Then
               M3(irby) = M1(irby) / M1(rrby)
            Else
               Jsi = M1(rjby) * M1(irby)
               Jsi = Jsi / M1(rrby)
               M3(ijby) = M1(ijby) - Jsi
            End If
         Next Jby
   Next Iby
      '-----------------------------' 矩阵乘法 M5 = M4 × M2
   For Iby = 1 To Nbyc
         I0by = Iby - 1
         I0by = I0by * Nbyc
         For Jby = 1 To Nbyc
            Ijby = I0by + Jby
            M5(ijby) = 0
            For Kby = 1 To Nbyc
                  K0by = Kby - 1
                  K0by = K0by * Nbyc
                  Kjby = K0by + Jby
                  Ikby = I0by + Kby
                  Jsi = M4(ikby) * M2(kjby)
                  M5(ijby) = M5(ijby) + Jsi
            Next Kby
         Next Jby
   Next Iby
      '---------------------------' 用M3,M5覆盖矩阵M1,M2
   For Iby = 1 To Nbyc
         I0by = Iby - 1
         I0by = I0by * Nbyc
         For Jby = 1 To Nbyc
            Ijby = I0by + Jby
            M2(ijby) = M5(ijby)
            M1(ijby) = M3(ijby)
            If Iby = Jby Then  ' 重建单位矩阵 M4
               M4(ijby) = 1
            Else
               M4(ijby) = 0
            End If
         Next Jby
   Next Iby
Next Rby
'-------------------------------' 矩阵乘法 B = M1×M2
For Iby = 1 To Nbyc
      I0by = Iby - 1
      I0by = I0by * Nbyc
   For Jby = 1 To Nbyc
         J0by = Jby - 1
         J0by = J0by * Nbyc
         Ijby = I0by + Jby
         B(ijby) = 0
         For Kby = 1 To Nbyc
            K0by = Kby - 1
            K0by = K0by * Nbyc
            Ikby = I0by + Kby
            Kjby = K0by + Jby
            Jsi = M1(ikby) * M2(kjby)
            B(ijby) = B(ijby) + Jsi
      Next Kby
   Next Jby
Next Iby
End Sub

'*******************************************************************************
'******************************* 打印矩阵 **************************************
'*******************************************************************************
Sub Matprint(byval Ksi As Single)
Jby = 0 : Kby = 0
If Ksi = A(1) Then
   Kby = 1
Elseif Ksi = B(1) Then
   Kby = 2
Elseif Ksi = C(1) Then
   Kby = 3
End If
For Iby = 1 To Mbyc
      Select Case Kby
            Case 1
                  Msi = A(iby)
            Case 2
                  Msi = B(iby)
            Case 3
                  Msi = C(iby)
            Case Else
               Exit Sub
      End Select
      If Msi >= 0 Then
         Print " ";
      End If
      Print Fusing(msi , "#.##") ; "  ";
      Incr Jby
      If Jby = Nbyc Then
         Print
         Jby = 0
      End If
Next Iby
End Sub

Proteus 仿真输出结果：

欢迎光临 (http://www.51hei.com/bbs/)