浅谈对矩阵理论、数值分析及其matlab应用的理解-4

《浅谈对矩阵理论、数值分析及其matlab应用的理解-4》由会员分享，可在线阅读，更多相关《浅谈对矩阵理论、数值分析及其matlab应用的理解-4（36页珍藏版）》请在装配图网上搜索。

1、3数值分析数值分析是研究用计算机求解各种数学问题的数值方法及其理论的一门学科。数值分析也称为数值计算方法、科学计算。随着计算科学与技术的进步和发展，科学与工程计算的应用范围已扩大到许多的学科领域，新的、有效的数值方法不断出现，形成了许多新型交叉学科。如：计算金融学、计算物理学、计算力学、计算化学等。继理论方法和实验方法之后，数值计算方法已成为科学研究的第三种基本手段，科学计算已经成为了人们进行科学活动必不可少的科学方法和工具。所以，数值分析既是一个基础性的，同时也是一个应用性的数学学科，与其他学科的联系十分紧密。 在计算机问世前, 由于受到计算工具的限制, 许多很好的数值方法难以应用。计算机出

2、现后, 其在通用性、标准化、大容量、高精度等方面都取得了惊人的成就, 特别是近20 年高速巨型计算机的成功研制与发展, 为数值分析方法提供了极好的机会, 大大提高了作为计算方法理论基础的数值分析方法课程的重要性与实用性。利用数值计算方法和计算机来解决科学与工程中的问题，通常称为科学与工程计算，或简称科学计算。科学计算是人们在社会实践中，特别是在科学和工程实践中普遍面临的问题：原问题 抽象为数学问题（建立原问题的数学模型） 求解该数学问题（必然涉及数值计算） 回答原问题。由于科学计算本身的迅速发展及不断取得成效，使得科学计算与传统的理论研究和实验研究并列称为当今研究的三大方法。而科学计算与具体科

3、学的交叉发展，又形成了诸如计算力学、计算物理、计算化学、计算生物学等新的边缘科学。这些学科与“数值分析”息息相关。误差来源误差是实验科学术语。指测量结果偏离真值的程度。对任何一个物理量进行的测量都不可能得出一个绝对准确的数值，即用测量技术所能达到的最完善的方法，测出的数值也和真实值存在差异，这种测量值和真实值的差异称为误差。分为绝对误差和相对误差。也可以根据误差的来源分为系统误差（又称偏性）和随机误差（又称机会误差）计量或测定中的误差是指测定结果与真实结果之间的差值。 数值计算普遍存在于科学研究和工程应用中，由于误差的存在，一般难以获得精确的计算结果，产生误差的原因主要有以下几个方面：1.模型

4、误差：数学模型对实际问题的仅是刻画；基于对实际问题近似描述的数学模型进行数值计算，例如利用函数的n阶Taylor展式计算函数值；2.观测误差：数学模型或计算公式中通常包含若干参数，这些参数往往是通过观测或实验得到的，这样得到的参数与其真值之间有一定的差异即所谓的观测误差，例如描述弹簧受迫振动的二阶线性常系数微分方程中的质量、阻尼系数和弹性系数等。更一般地：对物体的长宽高、电压、温度、速度的量测等。3.截断误差：许多数学运算是通过极限过程定义的，如微分、积分以及无穷级数求和等，由于计算机只能完成有限的算术运算和逻辑运算，所以在利用计算机进行计算是需要把无限的计算过程用有限的计算过程代替，由此产生

5、的误差成为截断误差；4.舍入误差：实际计算时只能按有限位进行，特别是里用计算机计算，由于计算机的有限位的限制，对参与运算的数据以及运算结果往往要进行舍入，例如利用公式计算圆的面积时， 需用有限的小数代替，由此产生的误差成为舍入误差。数值分析与计算主要研究截断误差和舍入误差，研究误差产生的原因、分析算法的误差以及控制计算过程中误差的扩散，由此把握计算结果的精度。a) 尽量减少计算步骤，简化计算；减少运算次数，既可以提高解题速度，又有可能使计算中的舍入误差积累减少。例：计算多项式的值。若采用逐项计算然后相加的算法所需的乘法次数为加法次数为n,若采用秦九韶算法则只需n次乘法和n次加法即可算出。b)

6、尽量避免两个相近的数减；在数值计算中，两个相近的数作减法时有效数字会损失。例：原始数据有9位有效数字，运算结果只有1位有效数字。c)尽量避免较大的数和较小的数相加；大数“吃掉”小数是指计算过程中，较小的数加不到较大的数中。这种现象有时会产生严重的后果。例：较小的数被较大的数“吃”掉了。d) 尽量避免绝对值较小的数做除数；在用计算机作运算时，绝对值很小的数做除数会溢出停机。而且当绝对值很小的除数稍微有一点误差时，对计算结果影响很大。例：设 ，计算结果起了很大的变化，因此，在计算中必须避免绝对值很小的数作除数。e) 选择数值稳定的算法定义：如果一个算法的计算结果对原始数据的误差（扰动）以及计算过程

7、中的舍入误差不敏感，则称该算法是数值稳定的，否则称为数值不稳定的。基本数值计算方法及其（结合Matlab相关函数）应用示例线性代数方程组的解法 线性代数方程组的解法在数值分析中占有极其重要的位置，一方面，在工程技术领域，经常以线性代数方程组作为其基本模型，例如，电学网络问题，热传导问题、质谱仪数据分析等；另一方面，在许多有效的数值方法中，求解线性代数方程组是其中关键的一步，比如样条差值法、矩阵特征值问题、微分方程数值解法等，都离不开线性方程组的求解。 求解线性代数方程组的数值方法可分为两大类：直接法和迭代法。直接法是指若没有舍入误差的影响，经过有限次算数运算可求得方程组的精确解的方法。由于实际

8、计算时舍入误差不可避免，故直接法得到的解并不精确。迭代法是利用迭代公式产生的向量序列去逼近精确解，当迭代收敛时，可得到满足精确要求的近似解。一般的，对于低阶的稠密线性方程组以及大型带形方程组的解，采用直接法比较有效，而对于大型稀疏（非常形）方程组则用迭代法求解比较有利。.1高斯消元法高斯消去法的求解思路：把一般的线性方程组化成（上或下）梯形的形式。例：使用Matlab软件，采用高斯消元法解方程Ax=b程序见附录1。.2高斯列主元消元法 高斯消元法是按方程组和未知元的自然顺序（即方程组对应增广矩阵的行和列的自然顺序）逐步消元的。高斯消元法中若主元的绝对值很小，将使相应的乘数的绝对值很大，从而导致

9、舍入误差的严重扩散（放大），严重影响所得解的精确度。改进措施是，调整方程组或未知元的次序（即对增广矩阵的行或列进行交换），使得系数矩阵中绝对值较大的元素作主元，这就是所谓的“选主元”技巧。列主元高斯消元法只是调整方程组（增广矩阵的行）的次序，不改变未知元（增广矩阵）的次序，与顺序高斯消元法相比，只是在每一步消元前增加一步按列选主元的工作。例：使用Matlab软件，采用列主元高斯消元法解方程Ax=b程序见附录1。.3直接三角分解法 高斯消去法的消元过程，从代数运算的角度看就是用一个下三角矩阵左乘方程组的系数矩阵，且乘积的结果为上三角矩阵，即例：使用Matlab软件，采用分解法解方程Ax=b程序见

10、附录1。.4迭代法迭代法是解线性方程组的另一类方法，特别适用与解大型稀疏线性方程组（阶数高、零元素多的方程组）。作为一种求解数值问题的通用方法，迭代法的基本思路是统一的，即利用设计好的迭代公式所产生的迭代序列去逐步逼近精确解。迭代法求解线性方程组的基本思想是a.不追求“一下子”得到方程组的解，而是在逐步逼近方程组的精确界的迭代过程中获得满足精度要求的近似值解，这一点与直接法不同；b.通过对问题的转化，避免（困难的）矩阵求逆运算。由迭代格式确定如下的迭代算法：对于给定的线性方程组，可以写成不同的（无穷多）迭代格式，有意义的（可用的）迭代格式应具有收敛性生成的序列收敛于方程组的解；好的迭代法应具有

11、较高的收敛速度。例：使用Matlab软件，采用Jacobi迭代法验证生成线性方程组是否收敛，若收敛求出谱半径，并改变谱半径的大小观察其与收敛速度的关系。保证Jacobi迭代法，Gauss-Seidel迭代法都收敛的情况下，选择最佳松弛因子计算线性方程组的解。经过实验发现运行时谱半径反别为0.5984、0.2256、0.1215时，收敛速度分别为37、26、20由上面三次运算结果可以看出谱半径越小收敛速度越快。SOR引入松弛因子后收敛速度明显加快。程序见附录2。3.4.2非线性代数方程组的解法一般非线性代数方程可有如下形式：设 在区间上连续。确定方程在该区间内是否有实根，如果有实根，求某个或全部

12、的实根（近似值）或方程的近似解。.1 二分法设在区间上连续、单调并且 。如方程首先要确定适当的包含根的区间，这可以依据闭区间上连续函数的介值定理来确定。注：二分法适用于一个方程的场合，收敛速度是线性的。迭代次数的估计：.2 黄金分割法在区间内取对称的两点： 使得求根（方法）程序如下： 按这种方法选取点和 ，每次去掉的区间长度至少是原区间长度的0.618倍，上述适用于一个方程的场合，收敛速度是线性的。迭代次数的估计：例：使用Matlab软件，采用黄金分割法解初始区间-3, 3， 收敛精度为的解。程序见附录3。.3 牛顿迭代法牛顿迭代法构造说明：设方程右端的函数 在所确定的包含唯一根的区间内具有一

13、阶连续导数 ，那么由一阶微分的性质，有 由上式可得从而得到牛顿迭代格式 以及迭代法: 。Newton迭代法的几何含义曲线在点 处的切线方程为 。 该切线与 轴的交点横坐标就是Newton迭代公式确定的：例：使用Matlab软件，利用函数Newton 迭代法求非线性方程组其中。的所有实根，误差不超过，并在途中画出所求实根对应的各点。所得图形如下，程序见附录4。差值方法建立函数表达式、计算函数值等问题，都是数学中最基础同时也很有实际需要的问题。差值问题就是这种问题中的一类近似处理方法。这些问题的产生，可能是已知函数太复杂，或者是已知函数只是一张表格，而这两者都不便使用。这时，我们就希望通过已知函数

14、的一些离散的值，按照一定的要求，构造出此函数的另外一个近似解表达式，一种常用的办法是从某个性质优良、便于结算的函数类中选出一个函数，如代数多项式，寻找的方法就是插值法。定理：给定个曲线上点 ，如果 ，互不相同，那么，在所有次数不超过次的多项式函数中，存在唯一的多项式函数，满足条件。.1拉格朗日插值法拉格朗日插值法是多项式插值中最基础的一种。构造插值多项式的基函数 ， 拉格朗日插值多项式：拉格朗日插值法的不足：在实际问题中，观测的数据可能会不断增加，如果用拉格朗日插值公式构造插值多项式，那么，每当增加数据就要重新计算多项式的系数，由此增加许多不必要的计算工作量。例：使用Matlab软件利用Lag

15、range差值方法求温度随水深（近似）变化函数表达式。在夏季，较大湖泊的水体按深度被跃变层分为上部的变温层和下部的均温层。水体的分层化对环境工程中污染问题的研究具有重要的意义，例如，有机物的分解将导致被跃变层隔离的底部水体中氧急剧减少。按温度随水深的变化曲线，跃变层位于水深处：。 现有美国普拉特湖（Platte Lake） 的一组数据：深度(): 0温度(): 所得图形如下图所示，程序见附录5。所得方程为：.2三次样条插值法三次样条插值（简称Spline插值）是通过一系列形值点的一条光滑曲线，数学上通过求解三弯矩方程组得出曲线函数组的过程。早期工程师制图时，把富有弹性的细长木条（所谓样条）用压

16、铁固定在样点上，在其他地方让它自由弯曲，然后沿木条画下曲线。成为样条曲线。从数学上抽象就得到三次样条函数这一概念。实际计算是还需要引入边界条件才能完成计算。边界通常有自然边界（边界点的导数为0），夹持边界（边界点导数给定），非扭结边界（使两端点的三阶导与这两端点的邻近点的三阶导相等）。a、确定三次样条函数的条件根据三次样条插值的要求，样条插值函数在每个小区间 ， ，上都是三次多项式，每个三次多项式有4个系数，总共需要确定个系数，因此，需要个独立条件才能保证唯一地确定满足要求的三次样条插值函数。已知：条件b给出了条件：， 条件c要求 在区间上具有二阶连续导数，所以，插值函数在中间插值节点 ， ，

17、处，必须满足条件：， ； ， ； ， ； 这样已经有个条件：还需要2个条件，才能保证对插值函数的唯一性要求；为此，通常在插值区间的端点处附加2个条件：第一类边界条件： 固定端点的斜率：固定边界条件： ，. ， 第二类边界条件： 给定端点的的二阶导数：自由边界条件： 。 第三类边界条件：周期性条件： ， 。 b、 三次样条插值函数的构造1. 设待求的三次样条插值函数 在各插值节点处的一阶、二阶导数，由于插值函数在各小区间，上都是三次多项式，所以二阶导数 是一次多项式，记 ， 利用Lagrange插值公式， 在各个区间上可表示为：， 此式的两端积分两次，得到样条函数每个区间上的表达式： ， 在每个

18、小区间的左端点 处 ， 和右端点 处 ， ，利用插值条件：得到一阶导数的表达式： ， (2由此确定了个系数： 。再利用一阶导数在中间个点处的连续性：，；得到个等式（方程）： 经进一步整理：， 最终归结为求解 阶的三对角线性方程组 的问题，其中： 对于自然边界条件，有例：使用Matlab软件利用三次样条差值方法求温度随水深（近似）变化函数表达式。（题目如上题）所得图形如下图所示，程序见附录6。所得方程为：。分析：根据常识我们可以知道，随着向湖面或湖底延伸，水温应该呈现自然的延伸状态，lagrange插值法没有展现这种趋势，而且在各个点之间变化比较大，而三次样条插值法加上自己设置的边界条件，所给出

19、的曲线较好地展示了这些特性，因此，相比较而言，三次样条插值法更适用于此场合，更具有效性。.3曲线拟合的最小二乘法给出一组离散点，确定一个函数逼近原函数，插值是这样的一种手段。在实际中，数据不可避免的会有误差，插值函数会将这些误差也包括在内。因此，我们需要一种新的逼近原函数的手段：不要求过所有的点（可以消除误差影响）；尽可能表现数据的趋势，靠近这些点。设近似函数为： 函数值 与观测值 之差称为残差 可以用残差来衡量近似函数的好坏。拟合问题的几何背景是寻求一条近似通过给定离散点的曲线，故称曲线拟合问题。最小二乘法一般形式 为线性无关的基函数，求驻点，令 即令 则法方程组可写成以下形式 函数空间的基

20、然后列出法方程 函数空间的基然后列出法方程 例：使用Matlab软件利用最小二乘法拟合法按要求做下题。炼钢是通过氧化降低含碳量的过程。钢液含碳量的多少直接影响冶炼时间的长短，因此，必需了解它们之间的关系。下表给出的是一组熔毕碳 （炉料熔化完毕时）与精炼时间（从炉料熔化完毕到出炉）的数据。含碳量:104121134147150177180191190204时间: 分钟 100125135155170185200205210235选择不同次数的多项式，试用该组数据，画出数据散点图和拟合曲线图形，比较拟合的效果，选择你认为适当的次数，并说明理由.所得图形如下图所示，程序见附录6当n=1时 当n=2时

21、 当n=3时 当n=4时 结果分析：因为所给数据近似线性分布当n=1或2时数据点均匀分布在拟合线的两侧，即n=1或2就可满足结果需求。.4最佳平方逼近多项式 所谓函数最佳逼近就是从指定的一类简单的函数中寻找一个和给定的函数“最贴近”的函数。从几何（空间）的角度看，函数最佳逼近就是从指定的函数集合（函数空间）中寻找一个函数（点），该函数和给定的函数（定点）最贴近（该点到定点的距离最小）。为此，要给出函数之间距离的概念和定义，这时要把函数看做点（函数空间中的点）。区间 上连续函数组全体 构成一个函数空间，可以看做该空间中的一个点；所有次数不超过次的多项式集合 也构成函数空间。显然有 例：使用Mat

22、lab软件利用最佳平方逼近多项式求函数在区间上距离不超过 的最佳平方逼近多项式。在同一窗口画出该函数和最佳逼近多项式曲线，写出该多项式的表达式。所得图形如下图所示，程序见附录7N =8整理得：数值积分和数值微分公式.1数值积分对于给定的函数，如果1. 的函数关系式比较复杂；2. 未知，而仅仅知道该函数在个点，处的函数值 .则希望能用相对简单的计算方法，求得在上的定积分的近似值。a、 两点公式（梯形公式）b、 三点公式（Simpson公式）c、复化的梯形求积公式d、复化的Simpson求积公式复化的求积公式是受“利用低次多项式分段插值”能够提高逼近精度且计算简单的事实启发，将上述的“在整个积分区

23、间上用一个插值多项式代替被积函数”的做法变为“用分段插值多项式代替被积函数”。f、Monte Carlo 求积法 函数在区间上的平均值为当那么有由此，利用在上服从均匀分布的随机数生成函数unifrnd(a,b,m,n)生成，代入。例：使用Matlab软件利用复化的梯形求积公式、复化的Simpson求积公式、Monte Carlo 求积法求解下题。下面是一处地质岩层断面上部边缘的深度测量数据。水平距离():0深度(): 1.1 试利用复化的梯形求积法求该组数据所在曲线与基准线（ 轴）在范围 内所围成图形面积.画出数据散点图和图形的示意图.所得图形如下图所示，程序见附录81.2 试利用复化的Sim

24、pson求积法求该组数据所在曲线与基准线（ 轴）在范围 内所围成图形面积.画出数据散点图和图形的示意图.所得图形如下图所示，程序见附录8利用Monte Carlo法求面积（注：各做 次试验，取平均值作为面积值）。程序见附录8结果分析：法一、选择相邻两点为一组按照梯形的面积公式即可求出相邻两点的面积，然后相加即可。法二、选择相邻三点为一组，中间的点位插入点，然后按照Simpson求积法求出相邻三点的面积，然后相加即可。.2数值微分对于给定的函数，如果1. 的函数关系式比较复杂；2. 未知，而仅仅知道该函数在个点，处的函数值 .则希望能用相对简单的计算方法，求得导数的近似值。基于上述考虑，选择的方

25、法之一是利用函数的插值多项式的导数作为函数 导数的近似值，例如Lagrange插值多项式，由因而有这里需要说明一点的是，尽管 和的函数值可能相差不多，但是它们的导数有可能相差很大。常微分方程的数值解法 微分方程（组）是科学研究和工程应用中最常用的数学模型之一。如揭示质点运动规律的Newton第二定律：和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓解析解或公式解，如一阶线性微分方程的初值问题：的解为：但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实：I. 绝大多数的常微分方

26、程和常微分方程组得不到解析解；II. 实际应用中往往只需要知道常微分方程（组）的解在（人们所关心的,感兴趣的）某些点处的函数值（可以是满足一定精度要求的近似值）；如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得解的表达式，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。一般的一阶常微分方程初值问题是指如下的一阶常微分方程的定解问题： 微分方程（组）的初值问题通常是对一动态过程演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。解的存在且唯一性条件（定理）：如果函数在区域上连续，

27、且存在（Lipschitz）常数,对自变量 满足不等式（Lipschitz条件）：那么初值问题（Cauchy 问题）存在唯一的连续可微解 。.1 Euler 方法假设要求在点（时刻） ， ， 处初值问题的解 的近似值。首先对式的两端积分，得对其右边，如果被积函数用积分下限处的函数值代替被积函数作积分（从几何上的角度看，是用矩形面积代替曲边梯形面积），则有进而得到下式给出的递推算法Euler 方法改进的Euler 方法对于的右边，如果被积函数用积分限和处的函数值的算术平均值代替（几何上，是用梯形面积代替曲边梯形面积），则有进而得到下式给出的递推算法：通常算法改进的上式方法比Euler 方法的精度

28、高，但是，按算法Euler方法求 时要解（非线性）方程，这是上式方法不如Euler上式算法比Euler方法精度高的优点；II 尽可能地利用Euler方法计算简单的长处；人们采取了如下的称之为改进的Euler 方法的折衷方案：预测 修正 .2 RungeKutta法RungeKutta法是按选取区间 上函数 变化率的个数的多少和截断误差的阶数 来区分的一系列方法，如I. 二阶的RungeKutta法（改进的Euler方法） II. 三阶的RungeKutta法 III. 四阶的RungeKutta法III.1 古典形式 III.2 Gill 公式(具有减小舍入例：使用Matlab软件利用Eule

29、r 法，4阶R-K公式解下题。给定如下常微分方程Cauchy 问题：1.分别取参数，初值和步长：，1.1 取 ，在区间 上分别利用Euler 法，4阶R-K公式和matlab的ode45求问题（1）的数值解；的曲线；2.1 分别取 和 ，重复1.1的实验和1.2的做法，给出你根据观察，对微分方程（1）的解 变化趋势的判断。注：上述的实验2.1和2.2 初值可按如下方法随机确定：实验结果及分析如下，程序见附录9运行结果：只分别去了他们最后一个结果 结果分析：由（）分别取 和，重复1.1的实验和1.2的做法得到的图可以看出当区间变大时，得到的解趋于常数1.有图像也可以看出由4阶R-K公式得到的曲线

30、比Euler 法得到的图形曲线等贴近Matlab得到的图象。.3 Matlab法 ode是Matlab专门用于解微分方程的功能函数。该求解器有变步长（variable-step）和定步长（fixed-step）两种类型。不同类型有着不同的求解器，其中ode45求解器属于变步长的一种，采用Runge-Kutta算法；和他采用相同算法的变步长求解器还有ode23。ode45表示采用四阶，五阶Runge-Kutta单步算法,截断误差为(x)3。解决的Nonstiff(非刚性)常微分方程。解决数值解问题的首选方法是ode45，若长时间没结果，应该就是刚性的，可换用ode23试试。Matlab中求微分方

31、程数值解的函数有七个：ode45，ode23，ode113，ode15s，ode23s，ode23t，ode23tb 。ode45 此指令以 Runge-Kutta (4,5) 演算法为基础，运算上是属单步骤计算運y(tn)，过仅需时间点前之解y(tn-1)。 ode45函数指令可作为第一阶段评估的数据，然后再思考是否采取下一步骤。ode23与ode45函数同以 Runge-Kutta 方程式为基础，并以Bogacki与Shampine配对。在温和刚性的情況，容许度较粗略时比ode45 更有效率。此指令也是单步骤解题。ode15s以数值微分方程(NDFs).为基础之可变顺序指令。可选择后向微分

32、方程BDFs, (或称为Gear 法)。ode15s是一個多步骤指令。若执行前认为问题特性属于刚性，亦或使用过ode45指令但无法达成目的时，可以试用f ode15s指令。ode23s以Rosenbrock修正方程式二阶为基础。由于它是单步骤解题指令，故在粗略容许范围下可能比ode15s 更有效率。它可能解一些ode15s效率低或无法解的刚性问题。ode23t采用自由间机的梯形法。在问题属于中等刚性或不需含括数值阻尼的特性时可以使用此指令解题。ode23tb执行TR-BDF2初期利用梯形法含Runge-Kutta 方程式而第二阶段采用后向微分法二阶。如同ode23s一样，在粗略容许度下此法比o

33、de15s有效率。语法T,Y = ode45(odefun,tspan,y0) T,Y = ode45(odefun,tspan,y0,options)T,Y,TE,YE,IE = ode45(odefun,tspan,y0,options) sol = ode45(odefun,t0tf,y0.)T,Y = ode45(odefun,tspan,y0)odefun 是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名tspan 是区间 t0 tf 或者一系列散点t0,t1,.,tfy0 是初始值向量 Simulink与信号处理T 返回列向量的时间点 Y 返回对应T的求解列向量T,Y = o

34、de45(odefun,tspan,y0,options)options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等T,Y,TE,YE,IE =ode45(odefun,tspan,y0,options)在设置了事件参数后的对应输出TE 事件发生时间 YE 事件解决时间IE The index i of the event functionthat vanishes. sol =ode45(odefun,t0 tf,y0.) sol 结构体输出结果例：使用Matlab软件利用matlab法解下列题目二阶常微分方程是一类典型的、有十分广泛应用背景的数学模型。令例1、考虑

35、如下二阶常微分方程Cauchy 问题：取求解（仿真）的终止时刻 中的某个值。令是如下特征方程的根: 。 参数分别取如下五组：， （4-1）， （4-2）， （4-3）， （4-4） ， （4-5） 试分别用Matlab的ode45求上述问题的数值解。并按每组参数（4-1）-（4-5），分别在在同一个图形窗口画出所求数值解在区间 的图像。1.2 针对（4-1）-（4-5）中的每组参数，观察所求数值解曲线随时间的演化趋势，说明演化趋势的特征；参数（4-4）对应的演化特征具有典型的意义，工程中用专门名词(?)刻画这一现象。 （4-1）-（4-5）中的每组参数，求初值问题（2）的公式解，并结合公式解再

36、讨论问题1.2. 实验结果及分析如下，程序见附录10 取不同值运行结果时 时时 时根据图形类型的趋势可将图形分为3中类型：即图形（1）、（5），图形（2）、（4）和图形（3）。图形的演化趋势与各参数密切相关。由图形（1）、（5）可知当C=1时，曲线很快趋于稳定。而C=0时，曲线将很不稳定。对比图（3）、（4），仅仅d的取值不同，由图（2）、（4）可知当d=时，则形成了周期式的震动。图（4）形成的曲线图形在工程上我们称为共振。（附程序见附录10）当时运行结果是例2给定Van der Pol微分方程组的初值问题： 。 其中是参数。令得到一个平衡点（状态）： 。令在平衡点处特征多项式：1. 依次取参

37、数：求解（仿真）的终止时刻：初值： 利用Matlab 函数ode45（适用于非刚性问题） 或ode15s（适用于刚性问题）求数值解。1.1在不同的图形窗口分别画出：I. 随时间的演化曲线；II. 分别以 为横、纵坐标，在平面（相平面）画出变化曲线（轨迹）。1.2 依据对I.和II.的观察，阐述曲线演化特征与参数 之间的关系。1.3 针对 的情形，取不同的初值 进行试验，并重复1.1的做法，阐述曲线演化特征与初值之间的关系。2. 依次取参数：求解（仿真）的终止时刻：初值： 。 2.1 利用Matlab 函数ode15s进行求数值解；2.2 重复1.1的做法；2.3 依据对随时间演化曲线的观察，说

38、明两者的变化速度有无明显差异，与参数 的增大是否有关系；2.4 按照参数（9）的取值顺序，取相同的初值和终止时刻，分别利用ode45 和ode15s求数值解，并利用函数tic 和toc记录求解用时。然后最小二乘法分别拟合上述解法用时toc与参数 的关系曲线。 在同一图形窗口画出这两条曲线，说明两者之间的主要差异。运行结果如图所示，程序见附录11 结果：当0时，曲线振幅在时间达到一定数值后趋于0，而在 平面内是不闭合的；当0时，曲线是周期震荡，在平面内是闭合的。在=0时，是共振，在平面内成椭圆形。 运行结果如图所示，程序见附录11当=0时，曲线形状与初值无关，除了没有振幅外，其他都是周期震荡，只是振幅不同。而在平面内曲线都是封闭图形。运行结果如图所示，程序见附录11运行结果如图所示，程序见附录11由图知：当 值变大时用ode45 和ode15s求解所需的时间明显不同，且ode15s求解所需的时间远小于ode45，所以当求解刚性问题而且 值较大时，用ode15s要明显好于ode45。两者出现差异，也主要是 值的原因。附录表3-5材料清单综合布线系统设备配置清单及报价7.2平面布置图