OpenCascade Matrix

OpenCascade Matrix

eryar@163.com

摘要Abstract：本文对矩阵作简要介绍，并结合代码说明OpenCascade矩阵计算类的使用方法。 

关键字Key Words：OpenCascade、Matrix、C++ 

一、引言Introduction

矩阵的研究历史悠久，拉丁方阵和幻方在史前年代已有人研究。作为解决线性方程的工具，矩阵也有不短的历史。1693年，微积分的发现者之一莱布尼茨建立了行列式论（theory of determinants）。1750年，克拉默又定下了克拉默法则。1800年代，高斯和威廉若尔当建立了高斯－若尔当消去法。 

从其发展历史可以看出，从行列式论的提出到克拉默法则总共历时50多年之久，而在我们《线性代数》的课本中，也就是一章的内容，学习时间可能只有一个星期。难怪这些概念这么抽象晦涩，也是情有可原的，哈哈。不过，当理解这些概念，也就掌握了一个强大的工具，因为它也使有些计算统一、简单。 

矩阵是高等代数中的常见工具，也常见于统计分析等应用数学学科中。在物理学中，矩阵用于电路学、力学、光学和量子物理中都有应用；计算机科学中，三维变换就是使用了矩阵的线性变换。矩阵的运算是数值分析领域的重要问题。 

在计算机图形相关的软件中，矩阵是必不可少的重要工具。如：OpenGL、DirectX中的变换都是用矩阵实现的；OpenCascade中就有矩阵计算类（math_Matrix）；OpenSceneGraph中也有矩阵的计算类（osg::Matrix）等等。矩阵在计算中主要用来表达线性变换，统一三维空间中物体的位置和方向变换。 

本文主要介绍OpenCascade中矩阵类的使用方法。 

二、矩阵的运算 Matrix Calculation

矩阵的运算主要有矩阵的加法、数与矩阵相乘、矩阵与矩阵相乘、矩阵的转置、方阵的行列式、共轭矩阵、逆矩阵。还有通过矩阵来求解线性方程组的解问题。 

1. 矩阵的加法 Matrix Addition

应该注意，只有当两个矩阵是同型矩阵时，才能进行加法运算。矩阵加法满足下列运算规律：设A、B、C都是m×n矩阵： 

2. 数与矩阵相乘 Matrix Multiplication

数乘矩阵满足下列运算规律（设A、B为m×n矩阵，λ、μ为数）： 

矩阵相加与矩阵相乘合起来，统称为矩阵的线性运算。 

3. 矩阵与矩阵相乘 Matrix Multiplication

并把此乘积记作：C = AB 

按此定义，一个1×s行矩阵与一个s×1列矩阵的乘积是一个1阶方阵，也就是一个数： 

由此表明乘积矩阵AB＝C的（i,j）元cij就是A的第i行与B的第j列的乘积。必须注意：只有当每一个矩阵（左矩阵）的列数等于第二个矩阵（右矩阵）的行数时，两个矩阵才能相乘。

矩阵的乘法虽不满足交换律，但仍满足下列结合律和分配律（假设运算都是可行的）： 

4. 矩阵的转置 Matrix Transpose

矩阵转置也是一种运算，满足下列运算规律（假设运算都是可行的）：

5. 方阵的行列式 Determinant of the Square Matrix

应该注意，方阵与行列式是两个不同的概念，n阶方阵是n×n个数按一定方式排列成的数表，而n阶行列式是这些数（也就是数表A）按一定的运算法则所确定的一个数。

由A确定｜A｜的这个运算满足下述运算规律（设A、B为n阶方阵，λ为数）： 

6. 逆矩阵 Inverse of Matrix

定义：对于n阶矩阵A，如果一个n阶矩阵B，使AB＝BA＝E，则说矩阵A是可逆的，并把矩阵B称为A的逆矩阵，简称逆阵。 

方阵的逆矩阵满足下列运算规律： 

7. 线性方程组的数值解法

在自然科学和工程中有很多问题的解决都归结为求解线性方程组或者非线性方程组的数学问题。例如，电学中的网络问题，用最小二乘法求实验数据的网线拟合问题，三次样条的插值问题，用有限元素法计算结构力学中的一些问题或用差分法解椭圆型偏微分方程的边值问题等等。 

求解线性方程组的直接法主要有选主元高斯消去法、平方根法、追赶法等。如果计算过程中没有舍入误差，则此种方法经过有限步算术运算可求得方程组的精确解，但实际计算中由于舍入误差的存在和影响，这种方法也只能求得方程组的近似解。目前这种方法是计算机上解低阶稠密矩阵的有效方法。 

高斯消去法是一个古老的求解线性方程组的方法，但由于它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。 

用高斯消去法解方程组的基本思想是用矩阵行的初等变换将系数矩阵化为具有简单形式的矩阵，即三角形矩阵，再通过回代过程，求出方程组的解。 

三、OpenCascade矩阵计算

在OpenCascade中，矩阵的计算类是math_Matrix，它可以有任意的行数和列数，实现了矩阵大部分的运算。以下通过程序实例来说明math_Matrix的使用方法： 

/*
*    Copyright (c) 2013 eryar All Rights Reserved.
*
*        File    : Main.cpp
*        Author  : eryar@163.com
*        Date    : 2013-07-17 21:46
*        Version : 1.0v
*
*    Description : Demonstrate how to use math_Matrix class.
*                  题目来自《线性代数》同济 第四版
*                  
*/

#include <math_Matrix.hxx>

#pragma comment(lib, "TKernel.lib")
#pragma comment(lib, "TKMath.lib")

int main(int argc, char* argv[])
{
    math_Matrix A(1, 3, 1, 3);
    math_Matrix B(1, 3, 1, 3);
    math_Matrix C(1, 3, 1, 3);

    // 1. Test Matrix Transpose and multiply.
    // 题目来自P53习题二 第3题
    // Set value of Matrix A
    A.Init(1);
    A(2, 3) = -1;
    A(3, 2) = -1;

    // Set value of Matrix B
    B(1, 1) = 1;    B(1, 2) = 2;    B(1, 3) = 3;
    B(2, 1) = -1;   B(2, 2) = -2;   B(2, 3) = 4;
    B(3, 1) = 0;    B(3, 2) = 5;    B(3, 3) = 1;

    C.TMultiply(A, B);

    cout<<"Matrix A: "<<endl;
    cout<<A<<endl;
    cout<<"Matrix B: "<<endl;
    cout<<B<<endl;
    cout<<"Matrix C: "<<endl;
    cout<<C<<endl;

    // 2. Test the Determinant of the Matrix.
    // 题目来自P12 例7
    math_Matrix D(1, 4, 1, 4);
    D(1, 1) = 3;    D(1, 2) = 1;    D(1, 3) = -1;    D(1, 4) = 2;
    D(2, 1) = -5;   D(2, 2) = 1;    D(2, 3) = 3;     D(2, 4) = -4;
    D(3, 1) = 2;    D(3, 2) = 0;    D(3, 3) = 1;     D(3, 4) = -1;
    D(4, 1) = 1;    D(4, 2) = -5;   D(4, 3) = 3;     D(4, 4) = -3;

    cout<<"Matrix D: "<<endl;
    cout<<D<<endl;
    cout<<"Determinant of Matrix D: "<<D.Determinant()<<endl;

    // 3. Calculate Inverse of the Matrix.
    // 题目来自P54 11题(3) 和P56 30题(1)
    math_Matrix E(1, 3, 1, 3);

    E(1, 1) = 1;    E(1, 2) = 2;    E(1, 3) = -1;
    E(2, 1) = 3;    E(2, 2) = 4;    E(2, 3) = -2;
    E(3, 1) = 5;    E(3, 2) = -4;   E(3, 3) = 1;

    cout<<"Inverse of the Matrix E: "<<endl;
    cout<<E.Inverse()<<endl;

    // P56 30题(1)
    math_Matrix F(1, 4, 1, 4);
    
    F(1, 1) = 5;    F(1, 2) = 2;
    F(2, 1) = 2;    F(2, 2) = 1;
    F(3, 3) = 8;    F(3, 4) = 3;
    F(4, 3) = 5;    F(4, 4) = 2;

    cout<<"Inverse of the Matrix F: "<<endl;
    cout<<F.Inverse()<<endl;

    return 0;
}

计算结果为：

 1 Matrix A:
 2 math_Matrix of RowNumber = 3 and ColNumber = 3
 3 math_Matrix ( 1, 1 ) = 1
 4 math_Matrix ( 1, 2 ) = 1
 5 math_Matrix ( 1, 3 ) = 1
 6 math_Matrix ( 2, 1 ) = 1
 7 math_Matrix ( 2, 2 ) = 1
 8 math_Matrix ( 2, 3 ) = -1
 9 math_Matrix ( 3, 1 ) = 1
10 math_Matrix ( 3, 2 ) = -1
11 math_Matrix ( 3, 3 ) = 1
12 
13 Matrix B:
14 math_Matrix of RowNumber = 3 and ColNumber = 3
15 math_Matrix ( 1, 1 ) = 1
16 math_Matrix ( 1, 2 ) = 2
17 math_Matrix ( 1, 3 ) = 3
18 math_Matrix ( 2, 1 ) = -1
19 math_Matrix ( 2, 2 ) = -2
20 math_Matrix ( 2, 3 ) = 4
21 math_Matrix ( 3, 1 ) = 0
22 math_Matrix ( 3, 2 ) = 5
23 math_Matrix ( 3, 3 ) = 1
24 
25 Matrix C:
26 math_Matrix of RowNumber = 3 and ColNumber = 3
27 math_Matrix ( 1, 1 ) = 0
28 math_Matrix ( 1, 2 ) = 5
29 math_Matrix ( 1, 3 ) = 8
30 math_Matrix ( 2, 1 ) = 0
31 math_Matrix ( 2, 2 ) = -5
32 math_Matrix ( 2, 3 ) = 6
33 math_Matrix ( 3, 1 ) = 2
34 math_Matrix ( 3, 2 ) = 9
35 math_Matrix ( 3, 3 ) = 0
36 
37 Matrix D:
38 math_Matrix of RowNumber = 4 and ColNumber = 4
39 math_Matrix ( 1, 1 ) = 3
40 math_Matrix ( 1, 2 ) = 1
41 math_Matrix ( 1, 3 ) = -1
42 math_Matrix ( 1, 4 ) = 2
43 math_Matrix ( 2, 1 ) = -5
44 math_Matrix ( 2, 2 ) = 1
45 math_Matrix ( 2, 3 ) = 3
46 math_Matrix ( 2, 4 ) = -4
47 math_Matrix ( 3, 1 ) = 2
48 math_Matrix ( 3, 2 ) = 0
49 math_Matrix ( 3, 3 ) = 1
50 math_Matrix ( 3, 4 ) = -1
51 math_Matrix ( 4, 1 ) = 1
52 math_Matrix ( 4, 2 ) = -5
53 math_Matrix ( 4, 3 ) = 3
54 math_Matrix ( 4, 4 ) = -3
55 
56 Determinant of Matrix D: 40
57 Inverse of the Matrix E:
58 math_Matrix of RowNumber = 3 and ColNumber = 3
59 math_Matrix ( 1, 1 ) = -2
60 math_Matrix ( 1, 2 ) = 1
61 math_Matrix ( 1, 3 ) = -8.88178e-017
62 math_Matrix ( 2, 1 ) = -6.5
63 math_Matrix ( 2, 2 ) = 3
64 math_Matrix ( 2, 3 ) = -0.5
65 math_Matrix ( 3, 1 ) = -16
66 math_Matrix ( 3, 2 ) = 7
67 math_Matrix ( 3, 3 ) = -1
68 
69 Inverse of the Matrix F:
70 math_Matrix of RowNumber = 4 and ColNumber = 4
71 math_Matrix ( 1, 1 ) = 1
72 math_Matrix ( 1, 2 ) = -2
73 math_Matrix ( 1, 3 ) = 0
74 math_Matrix ( 1, 4 ) = 0
75 math_Matrix ( 2, 1 ) = -2
76 math_Matrix ( 2, 2 ) = 5
77 math_Matrix ( 2, 3 ) = -0
78 math_Matrix ( 2, 4 ) = -0
79 math_Matrix ( 3, 1 ) = 0
80 math_Matrix ( 3, 2 ) = 0
81 math_Matrix ( 3, 3 ) = 2
82 math_Matrix ( 3, 4 ) = -3
83 math_Matrix ( 4, 1 ) = -0
84 math_Matrix ( 4, 2 ) = -0
85 math_Matrix ( 4, 3 ) = -5
86 math_Matrix ( 4, 4 ) = 8
87 
88 Press any key to continue . . .

线性方程组的计算方法请参考另一篇blog： 

使用OpenCASCADE的Math功能解线性方程组

http://www.cppblog.com/eryar/archive/2012/06/21/179629.html 

四、结论 Conclusion

OpenCascade的math_Matrix实现了矩阵的大部分运算，也可以作为一个小型的C++矩阵计算库来使用，在没有安装Matlab、Mathematica等软件的情况下，可以使用这个库进行一些计算。 

另外，对《计算方法》中算法感兴趣的读者可以参考OpenCascade其中代码的实现，加深对相关算法的理解。 

五、参考文献 Bibliography

1. 同济大学应用数学系 线性代数 高等教育出版社 2003 

2. 易大义，沈云宝，李有法 计算方法 浙江大学出版社 2002 

PDF Version: OpenCascade Matrix