第5章　数学： MATLAB数学计算 MATLAB的真正力量蕴藏在其对数学运算的精湛处理之中。无论是微积分、线性代数，还是更为复杂的数学分支，MATLAB都呈现了卓越的解决方案。精通MATLAB的工具和功能，读者的高等数学技能将提升至新的高度。 5.1初等数学 在初等数学的学习中，我们经常遇到离散数学和多项式问题，这些问题往往需要借助像MATLAB这样的计算软件来解决。掌握MATLAB，就意味着读者能以更高效、更精确的方式处理这些数学挑战。 5.1.1离散数学 虽然离散数学并非严格意义上的初等数学范畴，但我们这里提到的离散数学内容，主要涉及质数、约数、倍数、有理数等计算问题，这些都是高等数学和计算机科学的基础知识。掌握这些基础概念，对深入理解更高级的数学和计算机科学领域至关重要。离散数学例程如图51所示。 图51离散数学例程 说明： （1） gcd()函数和lcm()函数分别为求最大公约数（Greatest Common Divisor）和求最小公倍数（Least Common Multiple）的函数。 （2） perms()函数用于得到一个向量中所有元素的所有可能的排列，其实是一个应用场景较广的函数，比如常被应用于穷举算法以及一些集合论计算中。 （3） rat()函数用于取得输入量的有理分式，其参数为计算容差，根据不同的容差得到不同精度的有理分式逼近。 5.1.2多项式 在MATLAB中，多项式的表达采用向量形式，向量中的每个元素代表多项式相应次幂的系数，其中，向量的最后一个元素对应常数项。MATLAB内置了一些处理多项式的常用函数，例如求值函数polyval（）、求多项式根的函数roots（）、由根构建多项式的函数poly（），以及将多项式转换为符号表达式的函数poly2sym（）。多项式例程如图52所示，读者还可以利用MATLAB绘制多项式的图形，来丰富多项式的表达形式。 图52多项式例程 说明： （1） 在MATLAB中表示多项式的向量里，即使某些次幂的系数为0，也需要用0值来正确占位，以保持多项式结构的完整性。 （2） 绘制多项式图形的核心方法实际上是计算多项式在各个点上的值，并据此进行作图。 （3） 多项式显示函数poly2sym（）的本质是将多项式向量转换为符号矩阵。 5.2线性代数 线性代数，被称为“第二代数学模型”，它是构筑了现代科学的基石之一，并以其矩阵理论为核心，成为了高等数学思维中至关重要的一环。线性代数是研究“线性问题”的代数理论，线性是线性代数的第一特性。我们可以这样初步地理解线性——“线性就是简单性”，就是“可以任意叠加而不会产生意外”。线性代数无论在数学、自然科学，还是高等数学学习里都具有中心地位。线性代数还是高等数学思维的敲门砖，所谓高等数学，就是将事物抽象到本质结构，再通过结构来反向推知事物规律。因此，线性代数是最值得所有人深度体会与思考的一门数学。我们坚信，线性代数的思想是高等数学中最典型的范例。因此，学习线性代数无疑是进入高等数学世界的最佳途径。 5.2.1矩阵基础运算 在线性代数的学习中，矩阵的基本运算诸如提取对角线、计算迹（即对角线元素之和）、求秩，以及计算行列式等，是不可或缺的操作。在MATLAB中，这些操作的函数可以同时应用于符号矩阵和数值矩阵，如图53所示。 图53矩阵基本运算例程 说明： （1） diag()函数对于不同的输入会有不同的输出，如果输入的是矩阵，则输出对角线向量； 如果输入的是一个向量，则会依照该向量输出一个对角矩阵。 （2） rank()函数对于符号阵的计算，是按照最大秩获得结果的，即把所有符号都当作非零值来计算，实际意义微弱，因此一般并不常用。数值计算必然是存在精度的，精度大小与算法关系密切，比如MATLAB中求秩的算法是基于矩阵的奇异值分解的，在rank()函数里还可以后跟一个给定误差参数，即为机器精度。 5.2.2矩阵分解 MATLAB为用户提供了一系列全面的线性代数函数，涵盖了几乎所有的矩阵分解操作，这些函数同样适用于符号矩阵，如图54所示。 图54三角分解与特征分解例程 说明： （1） 在MATLAB中，对于符号阵的LU分解，可以得出与线性代数教材中一致的结果； 然而，由于考虑到算法稳定性的问题，在数值计算LU分解中，得到的所谓下三角阵（L）与定义中并不一致，而是需要进行PLU分解才能得到真正的单位下三角阵，这里需要多得到一个置换矩阵P，即A=P′LU。 （2） 特征分解在矩阵分解中至关重要，eig()函数提供分解算法得到特征对角阵D和特征向量矩阵V，实现满足AV = VD的特征阵分解计算。 5.2.3线性方程及矩阵的逆 线性方程组是许多问题的基础抽象模型，例如，坐标系变换、运动分析、解析几何中直线交点问题等，它们的应用范围极其广泛。MATLAB，这款根植于线性代数学科的科学计算软件，在解决线性方程问题上展现了其领先的技术实力。它以简洁的符号封装了适用于各种特殊情况的算法，并且在求解速度上做到了极致的优化，如图55所示。 图55解线性方程例程 说明： （1） 对于Ax=b形式的线性方程求解非常简洁，解即为“A＼b”，注意斜线的方向，记忆方法是，斜线方向指示A为分母项，是从等号左侧移动到右侧的分母项。 （2） inv()函数用于求矩阵的逆（Inverse），使用前述“斜线求解法”与“矩阵的逆求解法”在求解原理上大相径庭，后者无论从精度、速度、稳定性上都略逊一筹，从图55中的结果也可看出，前者求解精度足够高，可以显示为整数的形式（当然实际还是double类），而后者显示出了short形式，说明求解精度有限，可以使用“format long”命令将显示格式改为长型，则更为明显。 （3） 利用符号形式进行求解往往会得到不错的效果，并且可以避免计算精度的问题，simplify()函数是用于对符号表达式进行自动整理与化简的函数，可以看到化简后的解simpSlolve2与solve1完全一致。 5.3微积分 微积分是现代科学的核心数学基础之一，对于经典物理世界的建模几乎完全依靠在微积分的港湾里，这其中的底层逻辑在于，人类对于经典物理世界的理解就是“连续的”且“有因果关系的”，这两个特点正与微积分的特点相符，因而，微积分是对物理世界建模的第一工具。 微积分的英文Calculus其实就是“计算方法”之意，在没有计算机的时代，科学家与工程师把微积分看作很有效的计算工具，因此也诞生了许多计算微积分的手算算法，当然这些手算算法也不是万能的，而计算机软件工具促使这一切向前飞跃，MATLAB强大的符号计算引擎，已经可以瞬间得到许多人工无法得到的解析解，让微积分的计算从此不再是科学家与工程师需要考虑的难点。 5.3.1极限的艺术 微积分的基础是连续与极限，在MATLAB中有limit()函数可以用于求极限，且兼容符号计算与数值计算，如图56所示。 图56极限例程 说明： （1） 极限函数limit()的第三个参数为极限点的位置，可以是正无穷（inf）或负无穷（-inf），第四个参数可以定义所求极限为右单边极限（right）还是左单边极限（left）。 （2） 在本例中取a为1后使用fplot()函数作图时，会提示警告函数处理数组输入时行为异常，原因是在0点附近的虚线范围内，函数值的绝对值过大，导致无法正常作图。 （3） 对于多元函数的极限求解，仍然是使用limit()函数，原理是先对一元求极限，再对另一元求极限，下面两种代码形式等价： limit(limit(f,x,x0), y, y0) 和 limit(limit(f,y,y0), x, x0) 5.3.2导数： 原函数的“因” “导数”是微积分的核心概念，“导”的含义从组词中隐约可见，比如“导致”“引导”“导向”等，《史记·孙子吴起列传》中讲“善战者，因其势而利导之”，导的含义其实是“引领方向改变”，所以导数大则原函数方向改变多，导数小则原函数方向改变少，导数为零则方向不变。因而可以说，“导数是原函数的原因”。 MATLAB中求解导数使用diff()函数，极为简洁实用，如图57所示，这里与导数的英文Derivative并不一致，原因是diff()同时也是差分（Differences）函数，在3.4节的矩阵运算中提及过。 图57导数例程 说明： （1） diff()函数的第三个参数表示要求的是几阶导数，此项为空即为1阶导数。 （2） 本例展示了单引号与双引号在文本中的显示方法，单引号文本使用双引号引用，而双引号的文本使用单引号引用，这样不会引发代码歧义。 （3） 对于多元函数的偏导数，与前述极限的思路一致，代码形式如下： diff(diff(f,x,m), y, n) 或 diff(diff(f,y,n), x, m) 5.3.3积分： 原函数的“果” 积分（Integral）与导数是相反的运算，导函数是原函数的原因，而积分是原函数的结果。积分包括不定积分与定积分，不定积分是求导函数的原函数，求得结果是由一个函数代表的一族函数，而定积分是求导函数在一个区间内的面积，求得结果是一个值，代表着导函数在该区间内引发的变化，如图58所示。 图58积分例程 说明： （1） 积分函数int()对于不定积分与定积分通用，且并没有多次积分的输入参数，只能使用函数嵌套多层来实现。 （2） 不定积分的结果函数省略了一个常数C。 （3） 对于不可积的函数，即使是MATLAB也无能为力。 （4） 当定积分区间的一边是无穷值时，称为无穷积分，只需要将区间输入值设置为-inf或inf即可。 5.3.4泰勒展开： 多项式仿真工具 泰勒展开（Taylor Expansion）是微积分体系提供的又一伟大工具，它可以将任意包裹着外壳的函数都展开为多项式函数，可以说泰勒展开是用于模拟任意函数的一个多项式仿真工具，与此同时，由于展开的结果中各项对应着导数次数，因此相当于将原函数的变化原因依照主次进行了分解，对于了解原函数的本质有拨云见日的功效。对函数f(x)在x0点的泰勒展开的公式为 fx=fx0+f′x0x-x0+f″(x-x0)22!+…+f(n)(x0)n!(x-x0)n 等号右侧即为“泰勒多项式”，其最后一项为“多项式仿真的误差项”，是(x-x0) n的高阶无穷小项，如图59所示。 图59泰勒展开例程 说明： （1） 泰勒展开函数taylor()的截断阶数（order）项如果不予输入，则默认截断阶数为6，该阶数为绝对阶数，从示例中可见。 （2） 如绘制的图形所示，阶数越多仿真结果越准确，离展开点越近仿真结果越准确，这也是为什么sin(x)在x接近0时可以用x来近似表示的原因了。 5.3.5傅里叶展开： 频域上的简谐波仿真 “周期”是一个显然而又隐秘的概念，周期性的变化称之为“波”，波是物理世界的基本组成，而所有的波都以“简谐波”（正余弦）为宗，也就是说所有的周期函数，都可以分解为简谐波的组合。傅里叶展开（Fourier Expansion）就是对周期性函数的简谐波分解，也可以称其为简谐波动仿真系统。5.3.4节讲的泰勒展开相当于对函数在时域上的多项式仿真，而傅里叶展开则相当于对函数在频域上的简谐波仿真，两者均是解析物理世界的重要模型。 对于一个周期为T的函数f(x)，在其一个周期范围 ［-L,L］ 内，可以展开为如下级数形式： fx=a02+∑∞n=1ancosnπLx+bnsinnπLx 其中， an=1L∫L-LfxcosnπxLdx,n=0,1,2,…bn=1L∫L-LfxsinnπxLdx,n=1,2,3,… 对于非周期函数也一样可以进行傅里叶展开，展开的效果其实就是把 ［-L,L］ 范围认为是周期函数的一个周期，其余部分都默认自动进行了周期性拓展。在MATLAB中并没有专门的傅里叶展开函数，不过完全可以根据上述公式以及前述积分函数写出自己的傅里叶展开算法，如图510所示。 图510傅里叶展开例程 说明： （1） 例程中p表示展开的项数，在算法中对应于循环的次数； L为函数的半周期，即函数周期为2L。 （2） 如图形中所示，在一个周期范围内，傅里叶展开可以将函数值很好地逼近，级数的项数越多仿真精度就越高。 5.4插值与拟合 在科学及工程研究的过程中，有时会得到一些“成组数据”，这些数据可能是一维的或者多维的，并且代表着某几个变量之间的逻辑关系，这些成组数据被称为“样本点数据”。样本点数据可能不够密集，或者其中恰好没有想要得到的位点的数据，用户可以通过算法，向已知数据中“插入新的位点以得到新值”，称为“插值”。插值算法不改变输入数据，而是尽可能地使插入的新值“和谐”，让得到的新值从逻辑上“最可能”符合数据的内存逻辑关系。 从已知数据中得到新数据，往往并不足够，用户还希望得到数据与数据之间内存逻辑关系——函数关系，这样就需要一套算法“模拟出一个函数使之与样本数据相合”，因此称为“拟合算法”，拟合算法首先需要假设一个可能的“原型函数”，再根据样本点数据求出一套“可行参数”，保证与所有样本点都尽可能地“接近”，而且一般来说，都只是接近而不能完全通过样本点，但这样拟合的结果函数已经完全足够用于代表数据之间的真实规律了。 5.4.1一维插值 在已知点范围内进行插值称为“内插”，在范围外进行插值则称为“外插”； 从时间概念上讲，如果要插值的位点处于已知点之后，则称为“预报”。一维插值面向的是一维的已知数据的插值方法，如图511所示，输入的数据为向量形式。 图511一维插值例程 说明： （1） 一维插值函数interp1()对于输入已知点的向量并不要求单调，只要长度一致即可。 （2） 不同的插值方法对于计算结果与计算速度有很大的不同，如图511所示，一般来说，三次样条插值（spline）插值效果最好，计算速度最慢。主要的插值方法及说明如表51所示。 表51一维插值主要的插值方法及说明 方法 说明 连续性 最少已知点 'linear' 邻点的线性插值（默认方法） C0 2 'nearest' 距样本网格点最近的值 不连续 2 'next' 下一个抽样网格点的值 不连续 2 'previous' 上一个抽样网格点的值 不连续 2 'pchip' 邻点网格点处数值的分段三次插值 C1 4 'makima' 基于阶数最大为 3 的多项式的分段函数 C1 2 'spline' 邻点网格点处数值的三次插值 （比'pchip'需要更多内存和计算时间） C2 4 5.4.2二维网格数据插值 对于二维网格数据插值采用interp2()函数，同样也有几种插值方法，如图512所示，与一维稍有不同，代码格式为 interp2(x0, y0, z0, x1, y1, 'method') 图512二维网格数据插值例程 说明： interp2()函数只应用于“网格数据”，即要求在平面某范围内，所有网格点上均有数据。该函数常用方法及说明如表52所示。 表52二维网格数据插值常用设置及说明 方法 说明 连续性 每维度最少已知点 'linear' 邻点的线性插值（默认方法） C0 2 'nearest' 距样本网格点最近的值 不连续 2 'cubic' 邻点网格点处数值的三次卷积插值 C1 4 'makima' 基于阶数最大为 3 的多项式的分段函数 C1 2 'spline' 邻点网格点处数值的三次插值 （比'cubic'需要更多内存和计算时间） C2 4 5.4.3二维一般数据插值 对于非网格数据，MATLAB提供了一个面向更为一般数据的插值函数griddata(),如图513所示，代码格式为 griddata(x0, y0, z0, x, y, 'method') 图513二维一般数据插值例程 说明： （1） 四种插值方法如图513所示，其中，'natural'方法是基于三角剖分的自然邻点插值，支持二维和三维插值，该方法在线性与立方之间达到有效的平衡。另外还有一种'v4'方法，采用双调和样条插值方法，仅支持二维插值，且不是基于三角剖分，目前公认效果较好。 （2） 对于三维的网格样本点，可以使用函数interp3()或者更一般的n维插值函数interpn()； 而对于多维的非网格样本点，则使用与之对应的griddata3()和griddatan()函数。 5.4.4多项式拟合 多项式拟合是一种最常用和实用的拟合手法，前述的泰勒展开本质上也属于一种多项式拟合，只不过要求有已知函数表达式，而拟合行为只需要样本点数据作为输入即可，如图514所示。 图514多项式拟合例程 说明： （1） 多项式拟合函数polyfit()的第三个参数代表多项式拟合的阶数，阶数越多，精度越高。 （2） vpa()函数用于将输入数据转换为变精度数据（Variableprecision Arithmetic），第二个参数表示保留小数点后的数字位数，如若此处不采用vpa()函数，则会输出无必要的长分式数字，影响对于多项式的快速感观。 （3） 从图形中可见，在样本数据点范围内拟合的效果相当好，但是在此范围之外逐渐远离了原始的函数图像，这是拟合动作所无能为力的，如果想得到足够贴近真相的模型，就必须获得足够范围的样本点数据作为输入。 5.4.5最小二乘拟合： 拟合的优化之路 许多场景下，用户希望拟合的函数形式并不一定总是多项式，对于更一般的原型函数，MATLAB的优化工具箱（Optimization Toolbox）提供了函数lsqcurvefit()，基于最小二乘（Leastsquares, LSQ）原理实现曲线拟合（Curvefitting），该函数的输入原型可以是M函数名及匿名函数，可以高效地求解出要拟合的参数，如图515所示。 图515最小二乘拟合例程 说明： （1） 定义原型函数时，注意要将待定参数定义到函数的自变量集中。而在使用fplot()函数对拟合函数作图时，参数a已经不是未知量，需要将其从自变量集中取出。 （2） lsqcurvefit()函数的第二个参数为用户给定的参数初值，一般来说可以随意给定，如果希望提高运算速度，则可以尽量给出预测参数值。 （3） 从图515可见，在原型函数的结构选择准确的前提下，拟合算法可以很好地得到准确的拟合函数，且在样本数据范围外的部分也能很好地进行预测。 5.5代数方程与优化 由已知数与未知数通过代数运算组成的方程称为代数方程，线性方程作为简单的一次代数方程（组）已在5.2节线性代数中展示过解法，本节重点介绍非线性方程的求解问题。 优化（Optimization）问题是运筹学（Operations Research）的主要组成部分，以数学分析和组合学为主要工具，提到数学分析就离不开数学建模，而优化问题的数学模型往往就是代数方程（包括等式与不等式方程），因此优化计算问题最终大多转化为解代数方程问题。 规划（Programming）与优化向来是运筹学中概念区别不清晰的两个词，根据中国运筹学会引用《数学大辞典》的说法来看，优化与规划所包含的内容基本一致，因此也不加以区分； 习惯上宏观概念多称为优化，而具体的类别问题多称为规划，如线性规划、非线性规划等。 5.5.1代数方程的求解 解析求解函数为solve()，数值求解函数为vpasolve()，具体应用如图516所示。 图516代数方程例程 说明： （1） 将方程整理为等号右侧为零的形式后，只需要输入等式左侧即可，等号和等号右侧零的部分会自动补充。 （2） 对于解析解函数solve()，如果所求方程无法求出解析解，系统会提示，并且自动转为使用vpasolve()函数求得数值解。 5.5.2无约束优化 无约束优化是最简单的一类优化问题，其目标就是求某一函数的最小值（最大值只需要乘一个负号即化为最小值问题）； 在计算之前，对函数进行作图有一个直观的印象往往是有必要的，在MATLAB中，既有软件主体自带的函数fminsearch()，也有优化工具箱提供的等效函数fminunc()，函数名意为“find minimum of unconstrained multivariable function”，下面用著名的Rosenbrock香蕉函数来举例，如图517所示，该函数是一个常用来测试最优化算法性能的非凸函数。 图517无约束优化例程 说明： （1） 设置函数时，注意将多元自变量写为x(1)、x(2)这样的形式。 （2） 从图517中可见，fminsearch()函数不如优化工具箱提供的fminunc()函数收敛的速度快，对于绝大部分其他类型函数也是类似的效果，因此后者拥有更优异的计算性能，是优化问题的首选函数。 5.5.3线性规划： 高效决策工具 线性规划问题（Linear Programming），是“有约束优化问题”中最简单的一类问题，在线性规划中，目标函数和约束函数都是线性的，约束函数可能是不等式或等式及自变量上下界，数学描述为 min fTxs.t.Ax≤bAeqx=beqlb≤x≤ub 其中，记号“s.t.”翻译为“使得”，是“subject to”的简写，是极为常用的数学符号。MATLAB为线性规划问题提供了基于单纯形法的函数linprog()，如图518所示，常用代码格式为 ［x, fval］ = linprog(f,A,b,Aeq,beq,lb,ub) 图518线性规划例程 说明： （1） 计算结果显示，当x(1)为4且x(2)为1.25时，函数将在满足约束条件的前提下达到最小值fval为-9.25。 （2） 计算优化问题的前提是问题合理，如果在已知的约束条件下，函数本身并没有最值，那么即使是MATLAB也无能为力，一般会提示“问题是欠约束的”（Problem is unbounded）。 （3） 对于函数的参数中并有的参数，应使用空矩阵符号“［］”来占位。 5.5.4非线性规划 非线性规划问题，是“有约束优化问题”中较为复杂的一类，约束中不仅包含线性不等式、等式、自变量上下界，还可能包括非线性的不等式及等式，数学表达形式为 min f(x)s.t.Ax≤bAeqx=beqlb≤x≤ubCx≤0Ceqx=0 这一类问题都可以用MATLAB优化工具箱中提供的边界约束优化函数fmincon()来解决，函数名意为“find minimum of constrained nonlinear multivariable function”，例程如图519所示，函数常用代码形式为 ［x, fval］ = fmincon(fun,x0,A,b,Aeq,beq,lb,ub) 图519非线性规划例程 说明： （1） fmincon()函数的第二个参数为计算初值； 第三个参数A和第四个参数b为一组，表示约束边界为Ax≤b，第五个参数Aeq与第六个参数beq为一组，表示约束边界为Aeq x=beq。当没有第三个和第四个参数时，使用空矩阵符号“［］”占位。 （2） 注意该函数要求目标函数与约束函数必须都是连续的，否则可能会给出局部最优解而不是全局最优解。 5.5.5最大值最小化问题 最大值最小化问题（Minimax Constraint Problem）其实是一个实际应用中比较常见的问题，却由于问题类型不易被理解而很容易被忽略。举一个具体的例子来理解该问题，比如在一个战场中对于急救中心的选址，要求它到各阵地的运输伤员的时间fi(x)不可以太长，对于所有选址方案中能保证到达所有阵地中所需最长时间max fi(x)最小的，就是最佳方案。也就是说，伤员在运输过程中所消耗的每分钟都是损失，耽误时间越久就越危险，如何让最大危险降至最小，就是最大值最小化问题。该问题的数学描述与前述非线性规划非常类似： minxmaxi if(x)s.t.Ax≤bAeqx=beqlb≤x≤ubCx≤0Ceqx=0 而且在MATLAB中求解计算的形式也非常类似，如图520所示，函数使用fminimax()，其代码格式为 ［x,fval］ = fminimax(fun,x0,A,b,Aeq,beq,lb,ub) 图520最大值最小化例程 说明： fminimax()函数与前述fmincon()函数的参数几乎完全一致，不同的是，fminimax()函数要求输入目标函数组，这样才能对一组函数取最大值。 5.6微分方程 微分方程（Differential Equation,DE）分为常微分方程（Ordinary Differential Equation,ODE）和偏微分方程（Partial Differential Equation,PDE）。常微分方程中的“常”并不是“常系数”的意思，而是“正常”（Ordinary，通常的、普通的）的常，相对于“偏”（Partial，局部的、片面的），理解记忆两种微分方程的区别可以将关注点放在“偏”字上，即包含偏导数的微分方程为偏微分方程，不包含偏导数的微分方程即为常微分方程。目前，人类对于经典物理世界的认知，基本都是建立在微分方程的基础上的，这其中的底层逻辑在于，物理量的导数代表着该量的原因，积分代表该量的结果，而微分方程正是物理量的因果关系，这也是为什么许多学科的数学模型都是微分方程，而对于这些学科而言，解决变量之间的影响关系，就等同于求解微分方程。 5.6.1常微分方程解析解 常微分方程相比于偏微分方程更简单易于求解，对于线性常微分方程和一些特殊的低阶非线性常微分方程一般可以求得解析解，而绝大多数的非线性常微分方程是没有解析解的，这时需要退而求其次，求得方程的数值解，得到了求解区域的数值解，也就可以做出该解函数的图像，这无论是对理解函数性质还是分析函数数值都可以提供极大帮助，甚至可以再利用数据进行函数拟合，得到一个近似的解析解。 对常微分方程求解析解首先需要定义符号变量与符号函数，进而将微分方程定义为符号方程，如果方程有初始条件则也需要定义，然后直接使用dsolve()函数即可完成求解，如图521所示。 图521常微分方程解析解例程 说明： （1） 符号函数的定义即为函数字母加括号带自变量，如y(t)等，与数学中的书写完全一致； 注意符号方程中的等号是双等号（==），而不是赋值符号（=）； 仍然遵从一切都是矩阵的原则，方程组本质是也是方程的矩阵，需要用方括号括起来。 （2） 对于方程组的解，dsolve()函数将输出一个结构体，结构体的字段即为求解的函数。 （3） 常数C的角标完全取决于该程序之前的使用情况，角标值仅区分意义，同一个角标代表同一个常数。 （4） 对于无法求出解析解的方程，软件会提示“无法找到解析解”（Unable to find explicit solution），这时就要考虑使用数值解法了。 5.6.2常微分方程数值解 科学家很早就意识到解析解不是大多数微分方程的解决方案，于是纷纷研究探索数值解法，比如欧拉法、龙格库塔法、亚当斯莫尔顿法等，MATLAB内置的求解ODE的函数就是基于这些方法的，比如其中最常用而强大的ode45()函数就是基于四五阶龙格库塔法的算法。 MATLAB共内置8个ODE函数，对于初学者以及不专门研究它的用户，建议不必深究其中的区别，实际使用时选取方法是： 首选ode45()函数，如图522所示； 如果无法解算或速度极慢，则可尝试改用ode15s()函数； 如果ode45()函数可以求解，但是速度略慢而精度并不要求那么高，则可尝试改用ode23()函数。更为精微的函数区分可进入帮助文档中搜索“选择ODE求解器”查看。 图522常微分方程数值解例程 说明： （1） 设定的求解范围tspan是一个向量，向量的步长即为求解步长，步长越小则精度越高、算量越大。 （2） 注意设定微分方程时，必须将求解函数（如y）和自变量（如t）均列入参数列表中。 5.6.3微分方程Simulink求解 Simulink既是MATLAB的一个工具箱，同时也是几乎可以与MATLAB相并列的一款独立软件，就是因为其实在过于强劲，在某些方面甚至掩盖了MATLAB的光芒。由于Simulink也可以用于解一些微分方程，并且还十分方便，因此在这里作简要介绍。 Simulink将所有数据看成可以传递的信号，从一个模块单元中输出并输入另一个模块单元中，Simulink打包了各行各业可能会用到的大量的模块，以至于几乎没有可以新增的模块了，用户所要做的，就是将模块单元拖入模型中，连线并进行设置，即可运行计算。Simulink计算的对象称为“模型”，其意为对于现实物理系统的“数字孪生”，擅长所有可以称为系统的问题，微分方程（组）就是最简单且典型的一类系统，下面举一实例，解微分方程组： x·1=-x1+x2x·2=-2x2+1x·3=-x1x2 首先进入Simulink模块，方法是在命令行中输入simulink再按下Enter键即可，选择新建空模型（blank model），单击窗口上方模块浏览器按钮（library browser），选择需要的模块向模型中拖动，如图523所示，左侧即为搭建的模型，三个积分器的初值均设为0，使用Mux混路器模块将三个输入组合为一个向量输出，再使用Fcn函数模块计算方程等式的右侧，再将计算输出回路连接到对应积分器上，并且可以在x(t)输出端插入一个示波器Scope，这样便可以看到计算的输出结果。在本书配套的代码文件中，模型文件为pdeModel.slx，读者可以直接打开运行。 图523微分方程Simulink求解 从示波器图中可以看出，x1与x2均稳定在0.5的位置，而x3稳定为一条直线，这样即相当于求得了微分方程的数值解。Simulink的建模求解方式极为强大，虽然对于简单的小型问题，使用MATLAB解方程函数会更为快捷，然而对于求解较大规模的问题、模块化系统问题时，则更显示出Simulink的四两拨千斤，尤其对于时间延迟性微分方程来说，更是只有用Simulink才能实现求解。 5.6.4抛物椭圆型偏微分方程 常微分方程的解析解已经比较难求，到了偏微分方程，连求解数值解都困难了，在MATLAB的核心组件中，只对一类偏微分方程提供了数值解求解函数，这一类称为一阶抛物椭圆型PDE，数学形式如下： cx,t,uuxut=x-mxxmfx,t,uux+sx,t,uux 实际求解例程如图524所示。 图524抛物椭圆型偏微分方程例程 一阶抛物椭圆型PDE只是诸多偏微分方程中的一类，因此，MATLAB核心组件对于偏微分方程的求解还远远不够，实际上，MATLAB将其强大的偏微分方程求解功能集成在了优秀的“偏微分方程工具箱”（Partial Differential Equation Toolbox,PDET）中，而且还拥有图形用户界面，可以让用户非常方便地求解偏微分方程，因此建议初学者或不以此为研究的用户首选PDET。 5.6.5偏微分方程工具箱 一个标量物理量在空间中的表示为 u=u(x,y,z,t) 这也是目前人类对于世界的认知，即三维空间加一维时间（其中，时间单方向流动），如果把该物理量的值按颜色绘制在空间中，则会显示出一幅“云图”，这就是“三维空间中的物理场”，而场中每点的量值还会随着时间变化，这就是所谓的“瞬态物理场”。物理场是关于时间与空间的函数，导数表示物理量的原因，积分表示物理量的结果，因而该函数所要满足的客观规律即为偏微分方程，这就是为什么物理学几乎全部由偏微分方程所支撑起来，也是偏微分方程为什么会成为科学工程领域至关重要的环节的原因。 物理学中常见的偏微分方程最多包含二阶偏导数，而MATLAB提供的偏微分方程工具箱可以非常简易有效地求解“二阶偏微分方程”，还包括多元PDE以及PDE组，二阶PDE的数学形式如下： m2ut2+dut-·cu+au=f 其中为Nabla算子，表示梯度，·表示散度，如果c为常数，则梯度的散度其实是代表“所有非混合二阶偏导数”之意，即 ·cu=c2x21+2x22+…+2x2nu=cΔu 其中，Δ为Laplace算子。 二阶PDE与二阶代数方程有所类似，借用按圆锥曲线的分类来划分，包含如下三大类方程形式。 （1） 椭圆型方程（Elliptic），此时m=0，d=0，数学形式为 -·cu+au=f(x,t) 椭圆型方程没有对于时间的一阶及二阶偏导，因此定解问题中只有边界条件而没有初值条件，主要用来描述物理中的定常、平衡、稳定状态，最典型的就是泊松方程及二维拉普拉斯方程（梯度的散度恒为零），物理中常用的有定常状态的电磁场、引力场和反应扩散现象等。 （2） 抛物线方程（Parabolic），此时m=0，数学形式为 dut-·cu+au=f 包含场对于时间的一阶偏导，因此不仅需要边界条件，还需要一个初值条件，一般用于描述能量耗散系统，物理中常见的抛物线方程有一维热传导方程。 （3） 双曲线方程（Hyperbolic），此时d=0，数学形式为 m2ut2-·cu+au=f 包含场对于时间的二阶偏导，因此需要边界条件和两个初值条件，没有一阶时间偏导，可以理解为能量不耗散，一般描述能量守恒系统，常见的比如一维弦振动（波动）方程，波动方程的扰动是以有限速度传播的，因而其影响区和依赖区是锥体状的。 PDET对于上述这些类型的二阶PDE可以很方便地在二维或三维几何空间上完成求解，基本原理是采用三角形及四面体网格划分，采用有限单元法求解并对结果后处理得到物理场云图。 PDET提供两种使用方法，一是打开GUI（相当于一个软件界面）进行单击与输入操作，二是使用命令代码输入。后者的功能涵盖了前者，对于初学者可以使用GUI实现初步功能，再使用代码自动生成功能得到与操作对应的代码并在其基础上修改。下面按照操作步骤举一实际例子。 1. 打开PDET界面 方法是输入命令pdetool并按下Enter键，或在App界面中找到PDE Modeler。界面的菜单栏就代表着操作顺序： 选项（Options）、几何（Draw）、边界（Boudary）、方程（PDE）、网格（Mesh）、求解（Solve）、作图（Plot），如图525所示。 图525PDET初始界面 2. 选择PDE类型 位于选项（Option）菜单栏下的应用（Application）选项，也位于界面按钮行右侧，10个选项分别意为： 通用标量方程（Generic scalar）、通用系统方程组（Generic system）、机械结构平面应力场（Structural mechanics, plane stress）、机械结构平面应变场（Structural mechanics, plane strain）、静态电场（Electrostatics）、静态磁场（Magnetostatics）、交流电磁场（AC power electromagnetics）、直流电场（Conductive media DC）、热传导温度场（Heat transfer）、扩散（Diffusion）。 正确选择这些物理场，可以在界面上简化参数的输入，获得正确的引导提示。 3. 确定几何域 在GUI中指画出平面形状，有5项工具： 长方形、中心长方形、（椭）圆、中心（椭）圆和多边形，还有一项旋转工具，每创建一个图形，界面上Set formula后的文本框内都会出现该图形的“字母+序号”，字母有R（矩形）、Q（正方形）、E（椭圆形）、C（圆形）、P（多边形），加减号表示图形之间的布尔运算，可以自行修改顺序与符号以得到目标求解区域。 双击图形可以打开图形位置尺寸的准确设置窗口，如图526所示设置了中心为零点半径为1的正圆。 图526确定几何域 对应代码为 pdeellip(0,0,1,1,0,'E1'); 画图函数只有三个： 画（椭）圆函数pdeellip（）、画矩形函数pderect（）、画多边形函数pdepoly（）。还可以载入三维的标准模板库（Standard Template Library，STL）模型，作为三维几何域，此功能并不包含在界面上，需要使用代码来载入，函数为importGeometry()，代码形式如下： importGeometry(model,'BracketWithHole.stl'); 4. 设置边界条件 进入边界（Boundary）菜单栏，选择边界模式（Boundary mode）进入状态，这时，边界上会显示颜色与箭头，再进入边界条件设置（Specify boundary conditions），对于单个偏微分方程的边界条件有如下两类。 (1) 狄利克雷（Dirichlet）边界条件，也称第一类边界条件，该边界条件定义了边界处的值，数学形式为 hu=r 其中，h与r是关于x,t,u,u/x的函数，更为常见的是h与r都比较简单，可以将方程向右整理，令h为1即可。 (2) 诺依曼（Neumann）边界条件，也称第二类边界条件，定义了场在边界处的导数（变化率），数学形式为 n(cu)+qu=g 其中，q与g是关于x,t,u,u/x的函数，u/n表示x向量法向的偏导数。 如果是求解偏微分方程组的话，当然也不排除其中有不同的方程使用不同类型的边界条件，这时称为“混合边界条件”。 本例设置最简单的狄利克雷边界条件，让场量在边界处均为0,如图527所示。 图527边界条件对话框 不同的边界条件由不同的颜色区分，狄利克雷边界条件为红色，诺依曼边界条件为蓝色。 对于本例的设置，对应代码为 applyBoundaryCondition(model,'dirichlet','Edge',1:model.Geometry.NumEdges,'u',0) 5. 设置偏微分方程 选择菜单栏中的偏微分方程（PDE），选择设置偏微分方程（PDE Specification），进入PDE的设置界面，本例选择椭圆型方程（Elliptic），设置c=1, a=0, f=10，这也是软件的默认参数，如图528所示。 图528偏微分方程设置 如图528所示的界面上方已将方程的形式列出，可以参考对应位置的字母，方便设置。对应代码为 specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',1); 6. 设置网格 单击网格（Mesh）菜单栏，进入网格模式（Mesh mode）可以看到软件自动为区域划分了三角形网格，可以对网格进行一些更详细的设置和处理，比如细化网格（Refine mesh），细化前后的网格如图529所示。 图529网格细化前后对比 设置最大网格尺寸为0.1，对应代码为 generateMesh(model,'Hmax', 0.1); 如需显示网格划分的情况，可以使用如下代码对网格作图： pdemesh(model); 7. 求解方程 单击求解（Solve）菜单栏下的求解偏微分方程（Solve PDE），默认求解结果如图530所示。 图530方程解的云图 如图530所示，边缘上确实如边界条件设置的那样均为零。对应的代码如下，注意求解的直接结果为一个结构体： results = solvepde(model); u = results.NodalSolution; 如需绘图，使用pdeplot()函数，代码如下： pdeplot(model,'XYData',u) 8. 后处理作图 自动生成的图像往往不满足要求，这时可以单击作图选项（Plot selection），在界面上有多种功能可选，而且可以对作图的对象进行选择甚至输入，如图531所示为选择使用箭头标注场量的梯度方向，并使用jet作为配色方案，同时勾选了三维绘图选项。 图531多维数组/矩阵1 PDET既然可以解偏微分方程，就可以解决许多物理问题，比如热传递问题，如图532所示。 图532多维数组/矩阵2 还可以导入三维的STL模型，实现三维空间分析，如图533所示为一个支架结构的变形分析。 图533多维数组/矩阵3 PDET实质上是一个专为有限元分析而设计的仿真软件。在本质上，单一物理场问题可以归结为偏微分方程的求解，而多物理场耦合问题则涉及偏微分方程组的处理。作为多物理场有限元仿真的先驱，COMSOL Multiphysics实际上是MATLAB PDET的一个高度发展和演化形式，因此COMSOL不仅继承了求解偏微分方程的数学功能，还扩展了这些能力。如果遇到一些MATLAB难以攻克的偏微分方程挑战，可以转向使用COMSOL软件进行求解。如图534所示，COMSOL软件提供了一个数学接口，该接口能解决的偏微分方程类型比MATLAB更加丰富多样。在物理数学方程的数值求解方面，特别是在二维平面和三维空间问题上，更推荐使用像COMSOL这类专业的有限元软件。 图534COMSOL软件中提供的数学工具模块 5.7概率统计 统计学（Statistics）主要有两大方面的内容： 统计描述和统计推断。统计描述，用数字、函数、图像等描述一个概率分布，这个分布可能是总体的，也可能是样本的，用数字来描述分布的，称为分布度量。统计推断，是用来解释“样本分布”与“总体分布”之间的关系的，它分为两大问题： 参数估计（Parameter Estimation）和假设检验（Hypothesis Test）。 5.7.1概率分布 连续随机变量的概率分布函数F(x)，表示随机变量ξ满足ξ