方程左右两侧差值的“最小二乘法”来表示。 MATLAB矩阵计算
通过inv(矩阵)命令可以很容易获得逆矩阵。矩阵方程的解以矩阵除法运算符“\\”表示：
>> A=[ 3 -1 0; -1 6 -2; 0 -2 10]; >> B=[1; 5; 26]; >> X=A\\B X = 1.0000 2.0000 3.0000
行列式
行列式多用在稳定理论和评估矩阵奇异性上。为什么需要知道行列式？大多数时候你是想知道行列式是否为零。但是，当行列式为零、或者数值上非常接近于零时，由于前面提到的“舍入误差”原因很难进行数值计算。这是奇异值分解的另一个实例。
MATLAB用det(A)命令计算奇异值。在MATLAB命令行中手工输入以下矩阵，或者从matrix2.dat文件中复制以下矩阵：
>> A=[0.45, -0.244111, -0.0193373, 0.323972, -0.118829; -0.244111, 0.684036, -0.103427, 0.205569, 0.00292382; -0.0193373,-0.103427,0.8295, 0.0189674, -0.011169; 0.323972, 0.205569,0.0189674, 0.659479, 0.197388; -0.118829,0.00292382,-0.011169, 0.197388, 0.776985]
用det命令可以计算出行列式：
>> det(A) ans = -1.9682e-008
主轴理论：特征值和特征向量
MATLAB有预置的计算矩阵特征值和特征向量的函数：
>> eig(A) ans = -0.0000 0.7000 0.8000 0.9000 1.0000
当使用以下形式调用eig()函数时，它会以矩阵V的列向量形式返回特征向量：
>> [V,D]=eig(A) V =
-0.6836 0.0000 -0.5469 -0.4785 -0.0684 -0.4181 0.6162 0.1831 0.4530 -0.4547 -0.0837 0.4003 0.6189 -0.6232 0.2479 0.5409 0.2582 -0.2415 -0.4042 -0.6474 -0.2416 -0.6272 0.4755 -0.1190 -0.5550
D =
- 0.0000 0 0 0 0 0 0.7000 0 0 0 0 0 0.8000 0 0 0 0 0 0.9000 0 0 0 0 0 1.0000
eigs()函数是eig()函数的变形，它能够计算稀疏矩阵的特征值和特征向量。在下一节中将结合COMSOL Multiphysics介绍它的应用。
如果矩阵A的行列式 非常接近于零，那么一个特征值也非常接近于零。特征向量是矩阵A的化零空间——给出映射到零的方向：
>> A*V(:,1) ans = 1.0e-007 * 0.2669 0.1633 0.0327 -0.2112 0.0943
所有特征向量可以通过乘以特征值等于矩阵A的特性来证明这一点，例如：
>> A*V(:,2)./ V(:,2) ans = 0.7000 0.7000 0.7000 0.7000 0.7000
在MATLAB中，“./”符号表示元素对元素的除法，上面的冒号表示矩阵所有的列。
由于系统接近奇异，所有包含该矩阵的解都很差，对这一点我们不应该感到
奇怪。例如：
>> B=[0; 1; 0; 1; 0]; >> A\\B ans = 1.0e+006 * 2.1487 1.3142 0.2631 -1.7001 0.7593
由于A的所有元素都是个位数，被除数B矢量元素也是个位数，方程(27)的解也应该由个位数组成，而不是百万量级。对于化学工程师，这意味着对于一个质量守恒系统，入口流速为1kg/hr，由于质量守恒限制出口也为1kg/hr，但是反应器中流速为1000000kg/hr。这是不可能的。但这就是由近奇异矩阵给出的计算结果。
奇异值分解（SVD）
奇异值分解在很多情况下可以给出一个较好的解。类似于特征值和特征向量的主轴理论，所有矩阵都只有唯一的分解：
A?U?diag?VT (28)
这里的U和V都是正交方阵，diag是包含奇异值的对角矩阵。在解出U，V和diag的基础上，(27)式可以轻易解得
A?1??UT?V??1/diagj?? (29)
U和V为正交矩阵意味着其转置矩阵也是其逆矩阵。对角阵的逆矩阵为对角元素
取倒数。所以求解的唯一问题就是何时一个或多个奇异值(diagj)接近于零。这时(1/diagj)是一个非常大的数，导致求解结果错误。一个比较好的近似方法是将这些奇异值的(1/diagj)设为零！
??UTb x?V??1/diagj?? (30)
几乎为零元素的替代向量应该在数量级上最小，近似满足方程。
在我们矩阵A的例子中，如果只有一个输出，MATLAB命令svd()给出奇异值，如果有三个输出，命令svd()给出：
>>[U,D,V]=svd(A) U =
-0.0684 -0.4785 0.5469 0.0000 -0.6836 -0.4547 0.4530 -0.1831 -0.6162 -0.4181 0.2479 -0.6232 -0.6189 -0.4003 -0.0837 -0.6474 -0.4042 0.2415 -0.2582 0.5409 -0.5550 -0.1190 -0.4755 0.6272 -0.2416 D =
1.0000 0 0 0 0 0 0.9000 0 0 0 0 0 0.8000 0 0 0 0 0 0.7000 0 0 0 0 -0.0000 1.0000 V =
-0.0684 -0.4785 0.5469 0.0000 -0.6836 -0.4547 0.4530 -0.1831 -0.6162 -0.4181 0.2479 -0.6232 -0.6189 -0.4003 -0.0837 -0.6474 -0.4042 0.2415 -0.2582 0.5409 -0.5550 -0.1190 -0.4755 0.6272 -0.2416
由于A为对称矩阵，必然有U＝V。
最小数量级的奇异值分解命令如下：
>> ss=[1. 1./0.9 1./0.8 1./0.7 0]; >> dinv=diag(ss); >> V*dinv*U’*B ans = 0.0893 1.2820 0.1479 1.0317 -0.2130
这是一个物理上更容易接受的答案，例如，前面讨论的质量守恒系统中的内部气体流速。
由于有限元方法基于矩阵，所以线性系统理论对于COMSOL Multiphysics建模非常重要。当产生的刚度矩阵变得接近奇异时，COMSOL Multiphysics可能不能给出满足要求的结果。MATLAB中的矩阵计算和稀疏化方法可轻易地用于检查COMSOL Multiphysics结果的合理性，也可以用于对系统自然动力学的特征分析。这些想法将在下一节中作为COMSOL Multiphysics的一个模型例子加以介绍。
6.1 非均匀介质中的热传导
经常遇到的稳态热传导方程，是可以分析求解的最简单的偏微分方程：泊松方程。但是对于复杂几何情况，会比较难处理。作者认为某些通解很难推导和收敛[5]。所以，数值解比通解更好，进一步讲，传热过程的任何改变都可能会破坏分析结果。本节我们考虑非均匀厚板中的一维传热问题，在板中有分布热源：
?Td?dT??k??f(x)dx?dx?
x?0 (31)
?1Tx?1?0打开耦合MATLAB/COMSOL Script的COMSOL Multiphysics模型导航栏：
按照表9中的步骤建立有分布热源的介质传热问题。根据计算结果可以画出一条几乎线性的曲线，斜率为－1。由于这是个线性问题，求解器日志显示两步迭代收敛。计算结果得T|x=0.5=0.474097，而解析解为T|x=0.5～0.475：
T?1?1312x?x36?x412 (32)
下面尝试k＝1－x/2。这种情况下也有解析解，但是由于数值复杂性，需要采用对数和分支切割。解析解为T|x=0.5～0.550，你的计算结果呢？
表9 非均匀介质中的热传导——分布热源情况 Model Navigator 选择1D空间维数，COMSOL Multiphysics: Heat Transfer: Conduction: Steady state analysis，设定因变量：u 选择Element：Lagrange-Linear，完成 Draw菜单 Specify objects: Line, Coordinates弹出菜单，x：0 1 名称：interval，完成 Physics菜单： 选择边界1 Boundary settings 设定边界条件为：温度；T0＝1 选择边界2 设定边界：T1＝0，完成 Physics菜单： 选择域1 Subdomain settings 设定k=1; Q=-x*(1-x) 选择Init选项卡；设定T(t0)=1-x，完成 Meshing 点击三角形按钮，划分网格 Solver 点击求解按钮（＝） Post-processing 设定x＝0.5 Data display 值：0.474097[K]，表达式：T，位置：(0.5) 下面来看线性系统理论。打开“File”菜单，选择“Export FEM structure as fem”。这是我们第二次导出为FEM结构，有必要简单介绍一下。对于任何MATLAB/COMSOL Script变量，在命令行中输入变量名字，MATLAB/COMSOL Script将会给出变量值或其数据结构。试一下：
>> fem fem =
version: [1x1 struct] appl: {[1x1 struct]} geom: [1x1 geom1] mesh: [1x1 femmesh] border: 1 outform: ’general’ form: ’general’ units: ’SI’ equ: [1x1 struct] bnd: [1x1 struct]