第5章 MATLAB数学计算

数学计算是MATLAB最擅长的领域，MATLAB对于几乎所有数学领域都提出了非常杰出的解决方案，掌握了MATLAB软件工具，对于高等数学的理解和应用必然上升一个新台阶。
5.1初等数学
在初等数学中，很多离散数学和多项式的问题都需要类似MATLAB这样的计算软件参与解决。
5.1.1离散数学
离散数学严格地讲并不属于初等数学的一部分，只是这里讲的离散数学，主要是指关于质数、约数、倍数、有理数这样的计算，是许多高等数学和计算机学的基础，如图51所示。


图51离散数学例程



说明： 
(1) gcd()和lcm()函数分别代表最大公约数(Greatest common dvisor)和最小公倍数(Least common multiple)的计算。
(2) perms()函数用于得到一个所有可能的排列，其实是一个应用场景较广的函数，比如穷举算法以及一些集合论计算中。
(3) rat()函数用于取得输入量的有理分式，后跟参数为计算容差，根据不同的容差得到不同精度的有理分式逼近。
5.1.2多项式

多项式(Polynomial)是中小学阶段就熟悉的概念； 在MATLAB中，多项式由向量来表达，向量中倒数第n个数值对应着多项式的n-1次项，也就是说最后一位代表常数项。多项式常用的函数就是求值(polyval)、求根(roots)、由根反求多项式(poly)、显示多项式(poly2sym)，如图52展示了如何对多项式作图以及如何显示出多项式。


图52多项式例程



说明： 
(1) 注意，多项式向量中，如遇到系数为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)，使用前述“斜线求解法”与“矩阵的逆求解法”在求解原理上大相径庭，后者无论从精度、速度、稳定性上都略逊一筹，从图中的结果也可看出，前者求解精度足够高，可以显示为整数的形式(当然实际还是double类)，而后者显示为short形式，说明求解精度有限，可以输入format long将显示格式改为长型，则更为明显。
(3) 利用符号形式进行求解往往会收到不错的效果，并且可以避免计算精度的问题，simplify()函数是用于对符号表达式进行自动整理与化简的函数，可以看到化简后的解与solve1完全一致。
5.3微积分
微积分(Calculus)是现代科学的核心数学基础之一，对于经典物理世界的建模几乎完全依靠在微积分的港湾里，这其中的底层逻辑在于，人类对于经典物理世界的理解就是“连续的”且“有因果关系的”，这两个特点正中微积分下怀，因而
微积分是对物理世界建模的第一工具。

微积分的英文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()函数第三输入项表示要求的是几阶导数，此项为空即为一阶导数。
(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点的泰勒展开公式为： 



f(x)=f(x0)+f′(x0)(x-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)就是对周期性函数的简谐波分解，也可以称其为“简谐波仿真系统”。前面讲的泰勒展开相当于对函数在时域上的多项式仿真，而傅里叶展开则相当于对函数在频域上的简谐波仿真，两者均是人类解析物理世界的重要模型。

对于一个周期为T的函数f(x)，在其一个周期范围［-L,L］内，可以展开为如下级数形式： 


f(x)=a02+∑∞n=1ancosnπLx+ bnsinnπLx


其中，






an=1L∫L-Lf(x)cosnπxLdx,n=0,1,2,…

bn=1L∫L-Lf(x)sinnπ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'邻点的线性插值(默认方法)C02
'nearest'距样本网格点最近的值不连续2
'next'下一个抽样网格点的值不连续2
'previous'上一个抽样网格点的值不连续2
'pchip'邻点网格点处数值的分段三次插值C14
'makima'基于阶数最大为 3 的多项式的分段函数C12
'spline'邻点网格点处数值的三次插值

(比'pchip'需要更多内存和计算时间)C24

5.4.2二维网格数据插值
对于二维网格数据插值采用interp2()函数，同样也有几种插值方法，如图512所示，与一维稍有不同，函数格式为： 

interp2(x0, y0, z0, x1, y1, 'method')



图512二维网格数据插值例程


注意interp2()函数只应用于“网格数据”，即要求在平面某范围内，所有网格点上均有数据。该函数常用方法及说明如表52所示。


表52interp2()函数常用设置及说明



方法说明连续性每维度最少已知点

'linear'邻点的线性插值(默认方法)C02
'nearest'距样本网格点最近的值不连续2
'cubic'邻点网格点处数值的三次卷积插值C14
'makima'基于阶数最大为 3 的多项式的分段函数C12
'spline'邻点网格点处数值的三次插值

(比'cubic'需要更多内存和计算时间)C24

5.4.3二维一般数据插值
对于非网格数据，MATLAB提供了一个面向更为一般数据的插值函数griddata()函数，如图513所示，格式为： 

griddata(x0, y0, z0, x, y, 'method')



图513二维一般数据插值例程


说明： 

(1) 4种插值方法如图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.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()，也有优化工具箱(Optimization Toolbox)提供的等效函数fminunc()，函数名意为find minimum of unconstrained multivariable function，下面用著名的Rosenbrock香蕉函数来举例，如图517所示，该函数是一个常用来测试最优化算法性能的非凸函数。


图517无约束优化例程



说明： 
(1) 设置函数时，注意要将多元自变量写为x(1)、x(2)这样的形式。
(2) 从图中可见，fminsearch()函数不如优化工具箱提供的fminunc()函数收敛的速度快，对于绝大部分其他类型函数也是类似的效果，因此后者拥有更优异的计算性能，是优化问题的首选函数。
5.5.3线性规划
线性规划问题(Linear programming)，是“有约束优化问题”中最简单的一类问题，在线性规划中，目标函数和约束函数都是线性的，约束函数可能是不等式、等式、自变量上下界，数学描述为： 


min fTxs.t.
Ax≤b

Aeqx=beq

lb≤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非线性规划
非线性规划问题，是“有约束优化问题”中较为复杂的一类，约束中不仅包含线性不等式、等式、自变量上下界，还可能包括非线性的不等式及等式，数学表达形式为： 


minf(x)s.t.
Ax≤b

Aeqx=beq

lb≤x≤ub

C(x)≤0

Ceq(x)=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)最小的，就是最佳方案。也就是说，伤员在运输过程中所消耗的每一分钟都是损失，耽误时间越久就越危险，如何让最大危险降至最小，就是最大值最小化问题。该问题的数学描述与前述非线性规划非常类似： 


minx 
maxif(x)s.t.
Ax≤b

Aeqx=beq

lb≤x≤ub

C(x)≤0

Ceq(x)=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常微分方程解析解
常微分方程(Ordinary Differential Equation，ODE)相比于偏微分方程更简单易于求解，对于“线性常微分方程”和一些“特殊的低阶非线性常微分方程”一般可以求得解析解，而绝大多数的非线性常微分方程是没有解析解的，这时就退而求其次，可以求得方程的数值解，得到了求解区域的数值解，也就可以作出该解函数的图像，无论是理解函数性质还是分析函数数值都可以胜任，甚至可以再利用数据进行函数拟合(5.4节)的方法，得到一个近似的解析解。

对常微分方程求解析解首先需要定义符号变量与符号函数，进而将微分方程定义为符号方程，如果方程有初始条件也需要定义，然后直接使用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()。更为精微的函数区分可进入doc文档中搜索“选择ODE求解器”查看。


图522常微分方程数值解例程



说明： 
(1) 设定的求解范围tspan是一个向量，向量的步长即为求解步长，步长越小则精度越高向量越大。
(2) 注意设定微分方程时，必须将求解函数(y)和自变量(t)一起列入参数列表中。
5.6.3微分方程Simulink求解
Simulink既是MATLAB软件的一个工具箱，同时也是几乎可以与MATLAB相并列的一款独立软件，
Simulink的学习不是本书的主要内容，但由于Simulink也可以用于解一些微分方程，并且还十分方便，因此在这里简要介绍。
Simulink将所有数据看成是可以传递的信号，从一个模块单元输出并输入另一个模块单元中，Simulink打包了各行各业可能会用到的如此多的模块，以至于几乎没有可以新增的模块了，用户所要做的，就是将模块单元拖入到模型中，连线并进行设置，即可运行计算。Simulink计算的对象称为“模型”，其意为对于现实物理系统的“数字孪生”，擅长处理所有可以称为系统的问题，微分方程(组)就是最简单且典型的一类系统，下面举例，解微分方程组： 




x·1=-x1+x2

x·2=-2x2+1

x·3=-x1x2



首先进入Simulink模块，方法是在命令行中输入Simulink即可，选择新建空模型(Blank model)，单击窗口上方
“模块浏览器”按钮(Library browser)，选择需要的模块向模型中拖动，如图523所示，左侧即为搭建的模型，三个积分器的初值均设为0，使用Mux混路器模块将三个输入组合为一个向量输出，再使用Fcn函数模块计算方程等式的右侧，将计算输出回路给到对应积分器上，并且可以在x(t)输出端插入一个示波器Scope，这样可以直观看到计算结果。



图523微分方程Simulink求解



从示波器图中可以看出，x1与x2均稳定在0.5位置，而x3稳定为一条直线，这样即相当于求得了微分方程的数值解。Simulink的建模求解方式极为强大，虽然对于简单的小型问题，使用MATLAB解方程函数会更为快捷，然而对于求解较大规模的问题、模块化系统问题时，则更显示出Simulink的四两拨千斤，尤其对于“时间延迟性微分方程”来说，更是只有用Simulink才能实现求解。
5.6.4抛物椭圆形偏微分方程
常微分方程(ODE)的解析解已经比较难求，到了“偏微分方程”
(Partial Differential Equation，PDE)，连数值解都困难了，在MATLAB的核心组件中，只对一类偏微分方程提供了数值解求解函数，这一类称为“一阶抛物椭圆形PDE”，数学形式如下： 


c
x,t,u,
ux
ut=
x-m
x
xmf
x,t,u,ux
+s
x,t,u,ux



实际求解例程如图524所示。


图524抛物椭圆形偏微分方程例程


一阶抛物椭圆形PDE只是诸多偏微分方程中的一类，因此MATLAB核心组件对于偏微分方程的求解还远远不够； 实际上，MATLAB将其强大的偏微分方程求解功能集成在了优秀的“偏微分方程工具箱”(Partial Differential Equation Toolbox,PDET)中，而且还拥有图形用户界面(GUI)，可以让用户非常方便地求解偏微分方程，因此建议初学者或不以此为研究的用户首选PDE工具箱。
5.6.5偏微分方程工具箱
一个标量物理量在空间中的表示为： 

u=u(x,y,z,t)

这也是目前人类对于世界的认知，即三维空间加一维时间(其中时间单方向流动)，如果把该物理量的值按颜色绘制在空间中，则显示出一幅“云图”，这就是“三维空间中的物理场”，而场中每点的量值还会随着时间变化，这就是所谓“瞬态物理场”。物理场是关于时间与空间的函数，导数代表着物理量的原因，积分代表物理量的结果，该函数所要满足的客观规律即为偏微分方程，这就是为什么物理学几乎全部是由偏微分方程所支撑，也是偏微分方程为什么会成为科学工程领域至关重要的环节的原因了。

物理学中常见的偏微分方程最多包含二阶偏导数，而MATLAB提供的偏微分方程工具箱可以非常简易有效地求解“二阶偏微分方程”，还包括多元PDE以及PDE组。二阶PDE的数学形式如下： 


m2ut2+dut-(cu)+au=f


式中为Nabla算子，表示梯度，·表示散度，如果c为常数，则梯度的散度其实是代表“所有非混合二阶偏导数”之意，即： 


·cu=c2x21+2
x22+…+2x2nu=cΔu


式中Δ为Laplace算子。
二阶PDE与二阶代数方程有所类似，借用按圆锥曲线的分类来划分，包含如下三大类方程形式： 
(1) 椭圆形方程(Elliptic)，此时m=0，d=0，数学形式为： 


-·(cu)+au=f(x,t)


椭圆形方程没有对于时间的一阶及二阶偏导，因此定解问题中只有边界条件而没有初值条件，主要用来描述物理中的定常、平衡、稳定状态，最典型的就是泊松方程及二维拉普拉斯方程(梯度的散度恒为零)，物理中常用的有定常状态的电磁场、引力场和反应扩散现象等。
(2) 抛物线方程(Parabolic)，此时m=0，数学形式为： 



dut-·(cu)+au=f


包含场对于时间的一阶偏导，因此不仅需要边界条件，还需要一个初值条件，一般用于描述能量耗散系统，物理中常见的抛物线方程有一维热传导方程。
(3) 双曲线方程(Hyperbolic)，此时d=0，数学形式为： 


m2ut2-·(cu)+au=f


包含场对于时间的二阶偏导，因此需要边界条件和两个初值条件，没有一阶时间偏导，可以理解为能量不耗散，一般描述能量守恒系统，常见的比如一维弦振动(波动)方程，波动方程的扰动是以有限速度传播，因而其影响区和依赖区是锥体状的。
PDE工具箱对于上述这些类型的二阶PDE可以很方便地在二维或三维几何空间上完成求解，基本原理是采用三角形及四面体网格划分，采用有限单元法求解并对结果后处理得到物理场云图。
PDE工具箱提供两种使用方法： 一是打开GUI(相当于一个软件界面)进行单击与输入操作； 二是使用命令代码输入。后者的功能涵盖了前者，对于初学者可以使用操作界面实现初步功能，再使用“代码自动生成功能”得到与操作对应的代码并在其基础上修改。下面按着操作步骤举例。
1. 打开PDE工具箱
方法是输入命令pdetool。该界面的菜单栏就代表着操作顺序： 选项(Options)、几何(Draw)、边界(Boundary)、方程(PDE)、网格
(Mesh)、求解(Solve)、作图(Plot)，如图525所示。


图525PDE工具箱界面


2. 选择PDE类型
单击位于选项(Option)菜单栏下的应用(Application)选项，或者单击位于界面按钮行右侧的下拉组件，可以看到10个
PDE类型备选项： 通用标量方程(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');

画图函数只有3个： 画(椭)圆函数(pdeellip)、画矩形函数(pderect)和画多边形函数(pdepoly)。还可以载入三维的STL模型，作为三维几何域，此功能并不包含在界面上，需要使用代码来载入，函数为importGeometry()，代码形式如下： 

importGeometry(model,'BracketWithHole.stl');

4. 设置边界条件
选择菜单栏中的“边界”(Boundary)→“边界模式”(Boundary mode)进入状态，这时边界上会显示颜色与箭头，再进入边界条件设置(Specify boundary conditions)，对于单个偏微分方程的边界条件有如下两类： 
① 
狄利克雷边界条件(Dirichlet)，也称“第一类边界条件”，该边界条件定义了边界处的值，数学形式为： 


hu=r



式中h与r是关于x,t,u,u/x的函数，更为常见的是h与r都比较简单，可以将方程向右整理，令h为1即可。

② 诺依曼边界条件(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偏微分方程设置


该界面上方已将方程的形式列出，可以参考对应位置的字母，方便设置。对应代码如下： 

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方程解的云图


可见边缘上确实如边界条件设置的那样均为零。对应的代码如下，注意求解的直接结果为一个结构体： 

results = solvepde(model);

u = results.NodalSolution;

如需绘图，使用pdeplot()函数： 

pdeplot(model,'XYData',u)

8. 后处理作图
自动生成的图像往往不满足要求，这时可以选择“作图选项”(Plot Selection)，在界面上有多种功能可选，而且可以对作图的对象进行选择甚至输入，如图531所示为选择使用箭头标注场量的梯度方向，并使用jet作为配色方案，同时勾选了
“三维绘图”复选框。


图531多维数组/矩阵1



PDE工具箱既然可以解偏微分方程，那就可以解决许多物理问题，比如热传递问题，如图532所示。


图532多维数组/矩阵2


还可以导入三维STL模型，实现三维空间里的分析，如图533所示为一个支架结构的变形分析。


图533多维数组/矩阵3



可以看出，PDE工具箱其实就是一个有限元仿真软件，实际上，单物理场本质上就是偏微分方程，而多物理场耦合本质上就是偏微分方程组； 多物理场有限元仿真领域的领军软件COMSOL Multiphysics其实就是MATLAB的PDE工具箱独立发展进化的优秀产品，因此COMSOL软件里还包含解各种偏微分方程的数学功能，所以一些MATLAB难以求解的偏微分方程问题，可以使用COMSOL软件尝试求解，图534所示即为COMSOL软件中的数学接口，从中看出可求解的偏微分方程种类比MATLAB中更为全面。


图534COMSOL软件中提供的数学工具模块


5.7概率统计

“统计学”(Statistics)主要有两大方面的内容： “统计描述”和“统计推断”。统计描述，用数字、函数、图像等描述一个概率分布，这个分布可能是总体的，也可能是样本的，用数字来描述分布的，称为“分布度量”。统计推断，是用来解释“样本分布”与“总体分布”之间关系的，它分为两大问题： 参数估计(Parameter Estimation)和假设检验(Hypothesis Test)。
5.7.1概率分布

连续随机变量的概率分布函数F(x)，表示随机变量ξ满足ξ<x的概率，即F(x)=P(ξ<x) ，因此也被更形象地称为累积分布函数
(Cumulative Distribution Function,CDF)； 由于它可以通过简单减法得到随机变量落入任何范围内的概率，所以CDF是最为易用的概率描述函数。而概率密度函数(Probability Density Function,PDF)用于描述随机变量的输出值在某个取值点附近的可能性，一般记为p(x)，两者之间存在简单的积分关系： 


F(x)=∫x-∞p(t)dt


MATLAB提供了强大的统计和机器学习工具箱(Statistics and Machine Learning Toolbox，SMLT)，其中关于概率统计的相关函数均有涉及，如图535所示绘制了最典型的离散和连续概率分布——泊松和正态概率分布。


图535概率分布与概率密度例程


说明： 
(1) 概率密度函数pdf()与概率分布函数cdf()的输入格式非常清晰： 

pdf('name',x,A,B,C,D)

cdf('name',x,A,B,C,D)

其中name表示分布的名称， A~D即为对应分布的参数设置，输出向量与输入的向量x长度一致，因此可以直接作图。

(2) MATLAB共提供33种概率分布选项，基本上兼容了所有的概率分布，最常用的3种离散型及8种连续型分布如表53所示，如需查看所有分布选项，可以doc搜索pdf或cdf。
其中均匀分布也有离散型，即“等概率模型”(古典概型)，由于比较简单因此并没有对应函数。


表53最常用的3种离散型及8种连续型分布



类 别分布字段参数A参数B

离散型

二项分布'Binomial'n(试验次数)p(单次成功率)
泊松分布'Poisson'λ(均值)—
几何分布'Geometric'p(单次成功率)—

连续型

均匀分布'Uniform'a(下界)b(上界)
指数分布'Exponential'μ(均值)—

正态分布'Normal'μ(均值)σ(标准差)
瑞利分布'Rayleigh'b(规模参数)—
卡方分布'Chisquare'ν(自由度)—
伽马分布'Gamma'a(形状参数)λ(规模参数)
学生分布'T'ν(自由度)—
贝塔分布'Beta'a(形状参数)b(形状参数)


如果需要快速地可视化概率分布与概率密度曲线，MATLAB提供了一个概率分布函数App(Probability Distribution Function App)，相当于一个方便快捷的小软件，打开方式在命令行中输入disttool即可。如图536所示即为正态分布的CDF与PDF图像，界面上可以任意选择分布并输入参数，还可以使用鼠标交互式地捕捉坐标位置。


图536概率分布函数App


5.7.2伪随机数
在科学研究和统计分析活动中经常要用到随机数据，然而由计算机生成的随机数据，是通过算法并按照给定的分布规律计算出来的，称为“伪随机数”。

在MATLAB核心函数集中，提供了3个最常用的随机数生成函数，分别为： 产生均匀分布随机数的rand()、产生均匀分布随机整数的randi()、产生标准正态分布随机数的randn()。在统计和机器学习工具箱(SMLT)中还提供了更为通用和强大的随机数产生函数random()，它可以按任意概率分布规律产生随机数，其代码形式如下： 

R = random('name',A,B,C,D)

R = random('name',A,B,C,D,［size］)

可见与前述pdf()及cdf()函数所能处理的概率分布类型及分布参数输入形式完全一致，后跟向量为输出矩阵的尺寸规模。利用随机数产生函数绘制一些直方图如图537所示。


图537伪随机数例程


说明： 1e4为科学计数法，等价于10000，这里设置目的是产生一个由随机数组成的特定尺寸的向量。
5.7.3统计量分析
对于一组数据或一个分布，用一些数字量来作为某些整体性质的度量，这些量称为统计量，比如均值(mean)、方差(var)、标准差(std)、原点矩(无对应函数)、中心矩(moment)等，如图538所示为一个实例，其中构造了一个概率分布对象，并使用圆点表示法提取了参数，可以理解为解析解，是精确的统计量； 而通过随机量计算的统计量一定会存在偏差，当随机量数目越多时，偏差会有变小的趋势。


图538统计量分析例程



说明： 
(1) 均值函数mean()、方差函数var()、标准差函数std()的输入参数不仅可以是向量，还可以是矩阵，这时相当于对矩阵中所有的列向量进行计算，因而输出的是一个行向量，如果需要对矩阵中所有元素进行计算，则可以采用简单的代码形式： 

mean(x(:))

(2) 函数makedist()用于创建概率分布对象，该分布对象的属性中包含特定的分布参数，如mu和sigma等，还包含一些方法(函数)，也可以直接使用圆点表示法调用，比如mean、var及std等，与前述对于数据的分析函数完全一致。

(3) MATLAB没有提供原点矩计算函数，但根据定义易知k阶原点矩的代码如下： 

Ak = sum(x.^k)/length(x)

5.7.4参数估计

参数估计的应用场景是这样的，有时需要根据一组数据(样本数据)的基本形态做出推测，它应该是满足某分布规律的，比如正态分布规律，因此有理由推定总体数据也是满足正态分布规律的，那么这个规律表示成函数形式具体是什么样的，或者说分布函数里的参数应该是多少？这个计算过程就是“参数估计”。

本书观点认为“参数估计”(Parameter Estimation)这个名称略有晦涩，其实涵义比较接近于概率分布函数的拟合(Probability Distribution Fit)，只不过有一点在概念上需要注意，参数估计是用样本数据来估计总体的函数，因此样本数据的个数也是决定算法结果的关键。
对参数的估计，并不是只得到一个值，而是一个估计值加一个估计值的变化范围，称为“置信区间”(Confidence Interval)，与这个区间相对应的可信度称为“置信水平”(Confidence Level)，置信水平越接近1，则置信区间就会越小，一般最常用的置信水平是95%，因此MATLAB工具箱中提供的通用参数估计函数fitdist()即默认为95%的置信水平，工具箱还提供一些对应于特定分布的参数估计函数，比如正态分布估计函数normfit()、泊松分布估计函数poissfit()、均匀分布估计函数unifit()等，如图539所示，这些函数可以对置信水平进行设置，得到不同的置信区间。


图539参数估计例程


说明： 
(1) 通用参数估计函数fitdist()的常用代码形式如下： 

pd = fitdist(x,distname)

(2) normfit()函数的第二输入项意为显著性水平，显著性水平与置信水平的和为1； CI表示置信区间。

对于参数估计计算，MATLAB还提供了专用App——参数估计器(Distribution Fitter)，打开方式是在命令行中输入distributionFitter即可，界面友好、功能强大，如图540所示为本节例程中数据x的估计结果PDF及CDF。



图540参数估计器App


5.7.5假设检验
假设检验的应用场景是这样的，比如一个射击运动员说： “我射击的平均成绩是8环”，但是对于分析者来说，他说的话不知真假，只能称之为“假设”，如何证明这个假设，可以根据“大数定律”让他射击很多次取均值来判定，但实际情况很可能是样本数量并不足够多，这时就需要特定的检验算法，这也就是假设检验的过程。
假设检验(Hypothesis Test)先假设总体分布的参数，再用样本来检验这个参数的可信度(置信水平)，是一种持怀疑态度的基于反证法思想的验证，实例如图541所示。


图541假设检验例程


说明： 
对于疑似符合正态分布的数据，当已知标准差时，可以采用Z检验，函数为ztest()，若不已知标准差时，可以使用T检验。h为检验结论，为0时表明不拒绝假设，也就是接受假设； p为该检验的显著性水平； ci为置信区间。
本章小结
本章从初等数学、线性代数、微积分、插值与拟合、代数方程与优化、微分方程、概率统计7个方面深入介绍了MATLAB在数学领域的应用，本章内容的掌握有利于读者对于数学领域的理解更加深入，将数学的学习从记忆升级为理解与应用。