其中“matlab”告诉femlab打开一个matlab命令窗口，“path”根据COMSOL命令库建立matlab命令窗口。
首先运行COMSOL Multiphysics和模型导航栏窗口，然后根据表8中的步骤建立模型。该实例给了一个独立变量u和一个一维空间坐标x。h和r是系数模式中狄利克雷边界条件的两个句柄。如果你希望边界速度u为常数U0，可以通过设定h＝1，r＝U0来实现。点击三角形按钮生成默认网格（15个），然后点击两个嵌套三角形按钮加密网格（30个）。
表8 反应－扩散系统中两点边界值问题的常微分方程实例 选择1D空间维数，COMSOL Multiphysics: PDE modes: PDE, coefficient形式 Model Navigator 设置因变量：u 选择Element：Lagrange-Linear，完成 Specify objects: Line，Coordinates弹出菜单，x：0 1 Draw菜单 名称：interval，完成 选择域1 Physics菜单： 选择Dirichlet边界条件，设定h=1; r=1 Boundary settings 选择域2 选择Neumann边界条件，完成 选择域1 Physics菜单： Subdomain settings Meshing Solve菜单： Solver parameters 设定C=-1; f=0.833*u; da =0 选择Init选项卡；设定u(t0)=1-x 完成 点击三角形按钮进行网格划分 Stationary linear求解器默认设置，完成 General选项卡，设定解的形式为“Coefficient” 求解（＝） 选择边界2，输入表达式u。Reports窗口显示： 值：0.861167，表达式：u，边界：2 Post-processing： Point evalution 根据计算结果可以画出图5，显然边界条件需要满足：x＝0处u＝1，在x＝1处斜率减为零。但是COMSOL Multiphysics有没有按照我们设想的求解呢？
图5 稳态下的浓度曲线
现在用稳态非线性求解器再计算一次。首先查看求解器菜单中的日志，迭代13次收敛。你注意到最终值从0.86降到0.69了吗？为什么会这样呢？线性（系数型）求解器只在初始时刻u(t0)=1-x时计算一次R(u)，所以方程(24)只需要一次迭代。而非线性求解器在每次迭代时计算一次R(u)，计算时使用上次迭代的u值。所以非线性求解器在向收敛靠近时，经过几次迭代可能完全“忘记”了初时猜测。
检验一下这种解释对不对。改变初值将会改变稳定的线性解。返回子域设置对话框，设定初值为u(t0)=1，最终得到的u(x=1)是多少？再试一下稳态非线性求解器，能得到和其它初值一样的结果吗？
该例说明为方程选择正确求解器的重要性。如果函数f依赖于任何非独立变量，就应该选择非线性求解器。线性求解器更快，但是它假设偏微分方程的系数不依赖于非独立变量u（否则为非线性问题）。如果不确定，还是应该用非线性求解器。别忘了，(21)式满足R(u)=ku，是线性问题，但是COMSOL Multiphysics使用非线性求解器才得到正确解！收敛比较慢也是模型形式导致的结果——选中雅克比求解器选项，非线性求解器只需要两次迭代就可以得出正确结果。
我们一开始认为(22)式是该问题有限微分矩阵方程，但是实际上最后用(24)式描述COMSOL Multiphysics的有限元问题。因为我们这里使用拉格朗日线性微元，在这种特殊情况下有限元和有限微分矩阵一致。根据这一点，我们看一下MATLAB/COMSOL Script如何计算该COMSOL Multiphysics问题。
打开“File”菜单，选择“Export FEM Structure as fem”。这就把现在的计算结果转化为MATLAB/COMSOL Script工作空间的数据结构，然后用MATLAB/COMSOL Script预置的函数和命令来计算它。
在MATLAB/COMSOL Script工作空间输入以下命令：
>>x=fem.mesh.p; >>u=fem.sol.u; >>plot(x,u)
将会弹出一个包含网格节点u值的MATLAB/COMSOL Script图形窗口。不要怀疑你的图形失真。这是因为COMSOL Multiphysics存储网格点和对应结果的方式是为了使矩阵变得稀疏和紧凑。可以用MATLAB排序命令对网格点进行排序，进而使图形有意义。
在MATLAB工作空间输入以下命令：
>>[xx, idx]=sort(x); >>plot (xx, u(idx))
最后考虑MATLAB对有限元矩阵的访问。Fem文件结构没有包含有限元刚度矩阵，但是它包含重建这个矩阵的FEMLAB函数的必要信息。这对有限元方法非常重要，这个FEMLAB函数是assemble，输入以下命令：
>>[K,L,M,N]= assemble(fem); >>K’/15
这张图包含上次COMSOL Multiphysics的计算结果，类似图5。实际上，我们之所以能够快速搞清楚fem求解结果的结构，是因为这是只有一个独立变量的一维问题。多变量和多维产生的网格和结果的数据结构只能用COMSOL Multiphysics的工具/函数来读取。图6给出了由MATLAB命令spy(K’)产生的松散结构。很有必要将它与(23)式相互比较，因为能够清楚地看到具有统一网格划分的一维拉格朗日线性基元和生成矩阵方程的有限微分法是类似的。
图6 由MATLAB命令spy(K')产生的K'松散结构
下面看一下MATLAB中矩阵的松散结构，所有的元素只有1，-2和1，位于三条对角线上。这是有限元法的刚度矩阵，对于未知类型，它等于(23)式。如果返回“Subdomain settings”对话框的“element”选项卡，选择拉格朗日二次元素，再次求解并导出FEM，生成上面的K，你会发现虽然矩阵也很疏松，但是和拉格朗日线性元素完全不同。
练习
1.6 偏微分方程的系数形式有一项αu，令α＝0.833，f＝0，再次计算以上“反应
－扩散”的例子，比较稳态线性和非线性求解器得出的结果。你能解释为什么这样的表达式会导致这样的结果吗？这样的表达式对刚度矩阵K又有什么影响？如果把Da选为某特定值，例如Da△x2＝2，使得对角元素几乎为零，你认为求解时会出现什么困难？
6．方法4：线性系统分析
MATLAB和COMSOL Multiphysics的核心是线性系统分析。本节将简要回顾线性算法理论——本科生工程数学中的“矩阵方程”就是这类典型问题。好消息是不需要你自己动手编程处理矩阵，这就是MATLAB的用途：提供一个工程矩阵计算子程序库的用户界面。应当注意到，对于本章的例子，COMSOL Script也能很好地实现这个功能。以前科学计算主要集中在矩阵计算的高效和松散方法上。Golub和Van Loan[3]的书是一本很好的专家级矩阵计算指南，但是对于MATLAB入门水平，Hanselman和Littlefield[4]的书是个不错的选择。
简单来讲，标准的矩阵方程如下：
a11x1?a12x2?a13x3???a1NxN?b1a21x1?a22x2?a23x3???a2NxN?b2a31x1?a32x2?a33x3???a3NxN?b3?aM1x1?aM2x2?aM3x3???aMNxN?bM (26)
这里有N个与M方程有关的未知变量xj。系数aij和右侧常数项bi都是未知数。在工程应用中，通常将模型转化为这样的线性方程系统。例如，质量和能量守恒通常都会生成这样的线性方程系统。
求解可能性
当N＝M时，约束条件和未知数相同，所以可能求出唯一解xj。但是如果一个或者多个方程彼此线性相关（行退化），或者所有的方程都只含由某些特定变量组成的未知数（列退化），就不能得出唯一解。对于稀疏矩阵，行退化和列退化相同，退化方程会导致奇异。但是从数值求解讲，至少还有两点会导致错误：
? 如果方程彼此不是严格线性相关，一些方程也可能在计算过程中由于舍
入误差而导致变得非常接近线性相关。
? 求解过程中累积的舍入误差可能会“淹没”真实解，这在N比较大的时
候容易发生。这种情况下，计算过程没有出错，但是计算结果会不满足原始方程。 线性系统指南
没有“典型”的线性方程系统，但是一个粗略的想法就是舍入误差变得不可忽略：
? 当N在20－50之间时，可以用单精度常规方法求解，不需要对以上两
种错误进行修正。
? 当N在几百数量级时可以用双精度求解。
? 当N在几千数量级时，如果系数矩阵稀疏（几乎都是零），由于系数特
性可以求解。对于稀疏矩阵，MATLAB有特殊的数据形式和一套针对性的函数。
但是在工程实际和物理过程中，有很多问题本身就是奇异或者非常接近奇异的。你可能会发现N＝10也很难求解。用奇异值分解法有时可以解决奇异值问题，将其分解为非奇异值。
数值线性代数的常规任务
方程(26)可以写为紧凑矩阵方程形式：
A?x?b (27)
? 求解未知矢量x，这里A是方阵，b是已知矢量。 ? A为常数，求解多个b矢量时的解。 ? 计算方阵A的逆矩阵A－1。 ? 计算方阵A的行列式。
? 如果M 欠定系统。这种情况下可能没有解或者有多个解。解空间由特殊解xp乘以任何线性相关的N－M维矢量构成，称为矩阵A的化零空间。我们的任务是通过奇异值分解找出解空间。 ? 如果M>N，通常没有符合(26)的矢量解x。经常会出现这种过定系统， 我们的任务主要是寻找最接近于满足方程的折衷解。通常接近程度以(26)