数学建模稳定状态模型及数学建模算法大全现代优化算法简介

第十四章稳定状态模型虽然动态过程的变化规律一般要用微分方程建立的动态模型来描述，但是对于某些实际问题，建模的主要目的并不是要寻求动态过程每个瞬时的性态，而是研究某种意义下稳定状态的特征，特别是当时间充分长以后动态过程的变化趋势。譬如在什么情况下描述过程的变量会越来越接近某些确定的数值，在什么情况下又会越来越远离这些数值而导致过程不稳定。为了分析这种稳定与不稳定的规律常常不需要求解微分方程，而可以利用微分方程稳定性理论，直接研究平衡状态的稳定性就行了。本章先介绍平衡状态与稳定性的概念，然后列举几个这方面的建模例子。§1微分方程稳定性理论简介定义1称一个常微分方程（组）是自治的，如果方程（组）（1）中的，即在中不含时间变量。事实上，如果增补一个方程，一个非自治系统可以转化自治系统，就是说，如果定义，且引入另一个变量，则方程（1）与下述方程是等价的。这就是说自治系统的概念是相对的。下面仅考虑自治系统，这样的系统也称为动力系统。定义2系统（2）的相空间是以为坐标的空间，特别，当时，称相空间为相平面。空间中的点集称为系统（2）的轨线，所有轨线在相空间中的分布图称为相图。定义3相空间中满足的点称为系统（2）的奇点（或平衡点）。奇点可以是孤立的，也可以是连续的点集。例如，系统（3）当时，有一个连续的奇点的集合。当时，是这个系统的唯一的奇点。下面仅考虑孤立奇点。为了知道何时有孤立奇点，给出下述定理：定理1设是实解析函数，且系统（2）的奇点。若在点处的Jacobian矩阵是非奇异的，则是该系统的孤立奇点。定义4设是（2）的奇点，称（=1\*romani）是稳定的，如果对于任意给定的，存在一个，使得如果，则对所有的都成立。（=2\*romanii）是渐近稳定的，如果它是稳定的，且。这样，如果当系统的初始状态靠近于奇点，其轨线对所有的时间仍然接近它，于是说是稳定的。另一方面，如果当时这些轨线趋于，则是渐近稳定的。定义5一个奇点不是稳定的，则称这个奇点是不稳定的。对于常系数齐次线性系统（3）有下述定理。定理2设是系统（3）的通解。则（=1\*romani）如果系统（3）的系数矩阵的一切特征根的实部都是负的，则系统（3）的零解是渐近稳定的。（=2\*romanii）如果的特征根中至少有一个根的实部是正的，则系统（3）的零解是不稳定的。（=3\*romaniii）如果的一切特征根的实部都不是正的，但有零实部，则系统（3）的零解可能是稳定的，也可能是不稳定的，但总不会是渐近稳定的。定理2告诉我们：系统（3）的零解渐近稳定的充分必要条件是的一切特征根的实部都是负的。对于非线性系统，一般不可能找出其积分曲线或轨迹，也就不可能直接导出奇点的稳定性。为克服这一困难，在奇点附近用一个线性系统来近似这个非线性系统，用这个近似系统的解来给出这个奇点的稳定解。定义6设是系统（2）的一个孤立奇点。称系统在点几乎是线性的，如果在的Jacobian矩阵是非奇异的，即。设在的某邻域内连续，并有直到二阶连续偏导数，则由多元函数的Taylor公式，可将展开成，其中是一个常数矩阵，这样得到的线性系统（4）称为系统（2）的线性近似。一开始，人们以为总可以用线性近似系统来代替所研究的原系统。但后来人们发现，这种看法是不对的，或至少说是不全面的，非线性系统中的许多性质，在它的线性近似中不再保留。即使象零解稳定性这样一个问题，也要在一定条件下，才可用它的线性近似系统代替原系统来研究。关于这个问题，我们有下述定理：定理3如果系统（4）的零解是渐近稳定的，或不稳定的，则原系统的零解也是渐近稳定的或不稳定的。然而，如果系统（4）的零解是稳定的，则原系统的零解是不定的，即此时不能从线性化的系统来导出原系统的稳定性。系统（3）在其系数矩阵的行列式的条件下，可知是系统（3）的唯一的平衡点，它的稳定性由特征方程：的根（特征根）决定。定理4设线性系统（3）所对应的特征方程是其中，。设和是它的根，则当时关于奇点有下述结论：（=1\*romani），是稳定结点；（=2\*romanii），是稳定退化结点；（=3\*romaniii），是不稳定结点；（=4\*romaniv），是不稳定退化结点；（=5\*romanv），是不稳定鞍点；（=6\*romanvi），是稳定焦点；（=7\*romanvii），是不稳定焦点；（=8\*romanviii），是不稳定中心。定理5设非线性系统（5）中的和满足条件：（=1\*romani）在点的某邻域内存在连续的一阶偏导数。（=2\*romanii）存在常数，使得，（）又设系统（5）的一次近似系统（3）的特征方程的根没有零实部，则（5）式与（3）式的奇点的类型相同，并有相同的稳定性或不稳定性。§2再生资源的管理和开发渔业资源是一种再生资源，再生资源要注意适度开发，不能为了一时的高产“竭泽而渔”，应该在持续稳产的前提下追求最高产量或最优的经济效益。这是一类可再生资源管理与开发的模型，这类模型的建立一般先考虑在没有收获的情况下资源自然增长模型，然后再考虑收获策略对资源增长情况的影响。2.1资源增长模型考虑某种鱼的种群的动态。在建立模型之前，做如下的基本假设：（=1\*romani）鱼群的数量本身是离散变量，谈不到可微性。但是，由于突然增加或减少的只是单一个体或少数几个个体，与全体数量相比，这种增长率是微小的。所以，可以近似地假设鱼群的数量随时间连续地，甚至是可微地变化。（=2\*romanii）假设鱼群生活在一个稳定的环境中，即其增长率与时间无关。（=3\*romaniii）种群的增长是种群个体死亡与繁殖共同作用的结果。（=4\*romaniv）资源有限的生存环境对种群的繁衍，生长有抑制作用，而且这一作用与鱼群的数量是成正比的。记时刻渔场中鱼量为，我们可以得到所满足的Logistic模型：（6）其中是固有增长率，是环境容许的最大鱼量。由分离变量法求得方程（6）解为（6）式有两个平衡点，即，，其中是不稳定的，在正半轴内全局稳定。2.2资源开发模型建立一个在捕捞情况下渔场鱼量遵从的方程，分析鱼量稳定的条件，并且在稳定的前提下，讨论如何控制捕捞使持续产量或经济效益达到最大。设单位时间的捕捞量与渔场鱼量成正比，比例系数表示单位时间捕捞率，可以进一步分解分解为，称为捕捞强度，用可以控制的参数如出海渔船数来度量；称为捕捞系数，表示单位强度下的捕捞率。为方便取，于是单位时间的捕捞量为。常数，表示一个特定的捕捞策略，即要求捕鱼者每天只能捕捞一定的数量。这样，捕捞情况下渔场鱼量满足方程（7）这是一个一阶非线性方程，且是黎卡提型的。也称为Scheafer模型。希望知道渔场的稳定鱼量和保持稳定的条件，即时间足够长以后渔场鱼量的趋向，并且由此确定最大持续产量。在平衡点处有，方程（7）有两个平衡点，。显然，它们均是方程的解。在的情况下，是一正平衡点。（7）式可改写为（8）易知，当时，；时，，即平衡解是不稳定的，而是稳定平衡解。即在捕捞强度的情况下，渔场鱼量将稳定在的水平，因此产量（单位时间的捕捞量）也将稳定在的水平，即此时可获得持续收获量。当然，当时，，渔场鱼量将逐渐减少至，这时的捕捞其实是“竭泽而渔”，当然谈不上获得持续产量了。如何才能做到渔资源在持续捕捞的条件下为我们提供最大的收益？从数学上说，就是在或的条件下极大化所期望的“收益”，这里的“收益”可理解为产量，则问题就可以数学地叙述为下述优化问题：约束条件为。这里它可以归结为的二次函数的最大值问题。简单的推导不难得到最大持续捕捞强度为，最大持续产量为。捕捞强度是得到最大持续捕鱼量的策略。2.3经济效益模型当今，对鱼类资源的开发和利用已经成为人类经济活动的一部分。其目的不是追求最大的渔产量而是最大的经济收益。因而一个自然的想法就是进一步分析经济学行为对鱼类资源开发利用的影响。如果经济效益用从捕捞所得的收入中扣除开支后的利润来衡量，并且简单地设鱼的销售单价为常数，单位捕捞强度（如每条出海渔船）的费用为常数，那么单位时间的收入和支出分别为，单位时间的利润为利润是渔民所关注的焦点。因此在制定管理策略时所期望极大化的“收益”，这时就应理解为经济利润或净收入而不是鱼的产量。因而所讨论的问题就变成了在使鱼量稳定在的约束条件下的。即求的最大值。容易求出使达到最大的捕捞强度为最大利润下的渔场稳定鱼量最大利润下渔场单位时间的持续产量为最大可持续净收益与前一模型相比较可以看出，在最大效益原则下捕捞强度和持续产量均有减少，而渔场的鱼量有所增加。并且，减少或增加的比例随着捕捞成本的增长而变大，随着销售价格的增长而变小，这显然是符合实际情况的。2.4种群的相互竞争模型有甲乙两个种群，当它们独自在一个自然环境中生存时，数量的演变均遵从Logistic规律。记是两个种群的数量，是它们的固有增长率，是它们的最大容量。于是，对于种群甲有其中，因子反映由于甲对有限资源的消耗导致的对它本身增长的阻滞作用，可解释为相对而言单位数量的甲消耗的食物量（设食物总量为1）。当两个种群在同一自然环境中生存时，考察由于乙消耗同一种有限资源对甲的增长产生的影响，可以合理地在因子中再减去一项，该项与种群乙的数量（相对于而言）成正比，于是，种群甲增长的方程为（9）这里的意义是，单位数量乙（相对而言）消耗的供养甲的食物量为单位数量甲（相对）消耗的供养乙的食物量的倍，类似地，甲的存在也影响了乙的增长，种群乙的方程应该是（10）对可作相应的解释。在两个种群的相互竞争中，是两个关键的指标。从上面对它们的解释可知，表示在消耗供养甲的资源中，乙的消耗多于甲，对可作相应的理解。一般来说，之间没有确定的关系，在此我们仅讨论相互独立的情形。目的是研究两个种群相互竞争的结局，即时，的趋向，不必要解方程组（9）和（10），只需对它的平衡点进行稳定性分析。为此我们解代数方程(11)得到四个平衡点分别为，，，。为分析这些点的稳定性，需使用相空间的技巧。首先找出在平面上使或的区域。注意到，当时，但要使和，当且仅当类似地，当且仅当这样我们得到在平面中,直线(12)把平面划分为和两个区域。类似地，对进行分析得到（=1\*romani），当且仅当（=2\*romanii），当且仅当（=3\*romaniii）直线（13）将平面划分为和两个区域。两直线（12）和（13）之间的位置关系可以由下图的四种情况来说明。每种可能性对应于平衡点的稳定性说明如下：（=1\*romani），由图（b）知，两直线将平面划分为三个区域：（14）（15）（16）可以说明不论轨线从哪个区域出发，时都将趋向。若轨线从出发，由（14）式可知随着的增加轨线向右上方运动，必然进入；若轨线从出发，由（15）式可知轨线向下方运动，那么它或者趋向点，或者进入，但进入是不可能的。因为，如果设轨线在某时刻经直线（12）式进入，则，由式（9）、（10）不难看出。由式（15）、（16）知，故，表明在达到极小值，而这是不可能的，因为在中，即一直是增加的。若轨线从出发，由（16）式可知轨线向左下方运动，那么它或者趋向点，或者进入，而进入后，根据上面的分析最终也将趋向。综合上述分析可以画出轨线示意图。因为直线（12）式上，所以在（12）式上轨线方向垂直于轴；在（13）式上，轨线方向平行于轴。（=2\*romanii），类似的分析可知稳定。（=3\*romaniii），稳定。（=4\*romaniv），不稳定（鞍点）。因为轨线的初始位置不同，其走向也不同或趋于或趋于。根据建模过程中的含义，可以说明点稳定在生态学上的意义：①：意味着在对供养甲的资源的竞争中乙弱于甲，意味着在对供养乙的资源的竞争中甲强于乙，于是种群乙终将灭绝，种群甲趋向最大容量，即趋向平衡点。②：情况与=1\*GB3①正好相反。=3\*GB3③：因为在竞争甲的资源中乙较弱，而在竞争乙的资源中甲较弱，于是可以达到一个双方共存的稳定的平衡状态，这是种群竞争中很少出现的情况。=4\*GB3④：留作习题。§3Volterra模型意大利生物学家D'Ancona曾致力于鱼类种群相互制约关系的研究，在研究过程中他无意中发现了第一次世界大战期间地中海各港口捕获的几种鱼类占捕获总量百分比的资料，从这些资料中他发现各种软骨掠肉鱼，如鲨鱼、鲢鱼等我们称之为捕食者的一些不是很理想的鱼占总渔获量的百分比，在1914～1923年期间，意大利阜姆港收购的捕食者所占的比例有明显的增加：年代1914 1915191619171918百分比11.921.422.121.236.4年代1919192019211922 1923百分比27.316.015.914.810.7他知道，捕获的各种鱼的比例基本上代表了地中海渔场中各种鱼类的比例。战争中捕获量大幅度下降，当然使渔场中食用鱼（食饵）增加，以此为生的鲨鱼也随之增加。但是捕获量的下降为什么会使鲨鱼的比例增加，即对捕食者而不是对食饵有利呢？他无法解释这个现象，于是求助于著名的意大利数学家V.Volterra，希望建立一个食饵—捕食者系统的数学模型，定量地回答这个问题。3.l形成模型为建立这样的模型，我们分别用和记食饵和捕食者在时刻的数量。因为大海中鱼类的资源丰富，可以假设如果食饵独立生存则食饵将以增长率按指数规律增长，即有。捕食者的存在使食饵的增长率降低，设降低的程度与捕食者数量成正比，于是满足方程（17）比例系数反映捕食者掠取食饵的能力。捕食者离开食饵无法生存，若设它独自存在时死亡率为，即，而食饵为它提供食物的作用相当于使死亡率降低，或使之增长。设这个作用与食饵数量成正比，于是满足（18）比例系数反映食饵对捕食者的供养能力。方程（17）和（18）是在没有人工捕获情况下自然环境中食饵与捕食者之间的制约关系，是Volterra提出的最简单的模型。这个模型没有引入竞争项。3.2模型分析这是一个非线性模型，不能求出其解析解，所以我们还是通过平衡点的稳定性分析，研究的变化规律。容易得到方程（17）和（18）的平衡点为,(19)当然，平衡解对我们来说是没有意义的。这个方程组还有一族解，和，。因此，轴和轴都是方程组（17），（18）的轨线。这意味着：方程（17）、（18）在由第一象限出发的每一个解在以后一切时间都保持在第一象限内。当时，方程（17）、（18）的轨线是一阶方程的解曲线。用分离变量方法解得(20)是任意常数。因此，方程（17），（18）的轨线是由式（20）定义的曲线族，我们来证明这些曲线是封闭的。引理1当时，方程（20）定义了一组封闭曲线。证明我们首先来确定当时函数和的性状。利用微积分方法可以作出和的图形。如下图所示。若它们的极大值分别记作和，则不难确定满足，（21），（22）显然，仅当（20）式右端常数时相轨线才有定义。当时，，，将式（21）和（22）与（19）式比较可知正是平衡点，所以是相轨线的退化点。为了考察时轨线的形状，我们只需考虑的情况，其中。首先注意到：方程具有一个解和另一个解。因此，当或时，方程没有解。当或时，这个方程具有唯一的解，而对于，则具有两个解和。较小的解总是小于，较大的解总是大于。当趋于或时，和都趋向于。因此，当和是正数时，由（20）所定义的曲线都是封闭的。而且，这些封闭曲线中的每一条（除和以外），都不含（17）和（18）的任何平衡点。所以（17）和（18）的具有初始条件，的所有的解，都是时间的周期函数。也就是说，（17）和（18）的具有初始条件，的每一个解，都具有这样的性质：，，其中是某一正数。D'Ancona所用的数据实际上是捕食者的百分比在每一年中的平均值。因此，为了把这些数据同方程组（17）和（18）的结果进行比较，对于（17）和（18）的任何解，，我们必须算出和的“平均值”。值得注意的是，即使还没有准确地求得和，我们仍然能够算出这些平均值。引理2设，是（17）和（18）的周期解，其周期，和的平均值定义为，这时，，。换句话说，和的平均值是平衡解。证明把（17）的两端除以，得到，于是由于因此，，于是，。类似地，把（18）的两端除以，由0到积分，我们得到。下面，我们考虑渔业对于上述数学模型的影响。注意到渔业使得食饵总数以速率减少，而使得捕食者的总数以速率减少。常数反映渔业的水平；即反映了海上的渔船数和下水的网数。因此，真实的状态由下列修正的微分方程组来描述：（23）这个方程组同（17），（18）完全一样（当时），只是其中换成，而换成。因此，现在和的平均值是，平均说来，中等捕鱼量实际上会增加食饵的数目，而减少捕食者的数目。相反，捕鱼量的降低，平均说来，会增加捕食者的数目，而减少食饵的数目。这个值得注意的结果称为Volterra原理，它解释了D'Ancona的数据。值得注意的是Volterra模型是非常粗糙的，有兴趣的读者可以作进一步的探讨。习题十四1.单棵树木的商品价值是由这棵树能够生产的木材体积和质量所决定的。显然依赖于树木的年龄。假设曲线已知，为树木砍伐成本。试给出砍伐树木（更确切地说砍伐相同年龄的树木）的最优年龄。如果考虑到森林轮种问题，即一旦树木从某一处砍掉，这块土地便可以种植新树，假定各轮种周期具有相等的长度，试建模讨论最优砍伐轮种的森林管理策略的问题。2.如果两个种群都能独立生存，共处时又能相互提供食物，试建立种群依存模型并讨论平衡点及稳定性，解释稳定的意义。3.如果两个种群都不能独立生存，但共处时可以相互提供食物，试建模以讨论共处的可能性。4.如果在食饵—捕食者系统中，捕食者掠夺的对象只是成年的食饵，而未成年的食饵因体积太小免遭捕获。在适当的假设下建立这三者之间关系的模型，求平衡点。第十五章常微分方程的解法建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以肯定得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如，于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。§1常微分方程的离散化下面主要讨论一阶常微分方程的初值问题，其一般形式是(1)在下面的讨论中，我们总假定函数连续，且关于满足李普希兹(Lipschitz)条件，即存在常数，使得这样，由常微分方程理论知，初值问题(1)的解必定存在唯一。所谓数值解法，就是求问题（1）的解在若干点处的近似值的方法，称为问题（1）的数值解，称为由到的步长。今后如无特别说明，我们总取步长为常量。建立数值解法，首先要将微分方程离散化，一般采用以下几种方法：（=1\*romani）用差商近似导数若用向前差商代替代入（1）中的微分方程，则得化简得如果用的近似值代入上式右端，所得结果作为的近似值，记为，则有（2）这样，问题（1）的近似解可通过求解下述问题（3）得到，按式（3）由初值可逐次算出。式（3）是个离散化的问题，称为差分方程初值问题。需要说明的是，用不同的差商近似导数，将得到不同的计算公式。（=2\*romanii）用数值积分方法将问题（1）的解表成积分形式，用数值积分方法离散化。例如，对微分方程两端积分，得(4)右边的积分用矩形公式或梯形公式计算。（=3\*romaniii）Taylor多项式近似将函数在处展开，取一次Taylor多项式近似，则得再将的近似值代入上式右端，所得结果作为的近似值，得到离散化的计算公式以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式。其中的Taylor展开法，不仅可以得到求数值解的公式，而且容易估计截断误差。§2欧拉（Euler）方法Euler方法Euler方法就是用差分方程初值问题（3）的解来近似微分方程初值问题（1）的解，即由公式（3）依次算出的近似值。这组公式求问题（1）的数值解称为向前Euler公式。如果在微分方程离散化时，用向后差商代替导数，即，则得计算公式（5）用这组公式求问题（1）的数值解称为向后Euler公式。向后Euler法与Euler法形式上相似，但实际计算时却复杂得多。向前Euler公式是显式的，可直接求解。向后Euler公式的右端含有，因此是隐式公式，一般要用迭代法求解，迭代公式通常为（6）2.2Euler方法的误差估计对于向前Euler公式（3）我们看到，当时公式右端的都是近似的，所以用它计算的会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的所谓局部截断误差。假定用（3）式时右端的没有误差，即，那么由此算出（7）局部截断误差指的是，按（7）式计算由到这一步的计算值与精确值之差。为了估计它，由Taylor展开得到的精确值是（8）（7）、（8）两式相减（注意到）得（9）即局部截断误差是阶的，而数值算法的精度定义为：若一种算法的局部截断误差为，则称该算法具有阶精度。显然越大，方法的精度越高。式（9）说明，向前Euler方法是一阶方法，因此它的精度不高。§3改进的Euler方法3.1梯形公式利用数值积分方法将微分方程离散化时，若用梯形公式计算式(4)中之右端积分，即并用代替，则得计算公式这就是求解初值问题（1）的梯形公式。直观上容易看出，用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法。梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为(10)由于函数关于满足Lipschitz条件，容易看出其中为Lipschitz常数。因此，当时，迭代收敛。但这样做计算量较大。如果实际计算时精度要求不太高，用公式（10）求解时，每步可以只迭代一次，由此导出一种新的方法—改进Euler法。3.2改进Euler法按式（5）计算问题（1）的数值解时，如果每步只迭代一次，相当于将Euler公式与梯形公式结合使用：先用Euler公式求的一个初步近似值，称为预测值，然后用梯形公式校正求得近似值，即（11）式（11）称为由Euler公式和梯形公式得到的预测—校正系统，也叫改进Euler法。为便于编制程序上机，式（11）常改写成（12）改进Euler法是二阶方法。§4龙格—库塔（Runge—Kutta）方法回到Euler方法的基本思想—用差商代替导数—上来。实际上，按照微分中值定理应有注意到方程就有（13）不妨记，称为区间上的平均斜率。可见给出一种斜率，（13）式就对应地导出一种算法。向前Euler公式简单地取为，精度自然很低。改进的Euler公式可理解为取，的平均值，其中，这种处理提高了精度。如上分析启示我们，在区间内多取几个点，将它们的斜率加权平均作为，就有可能构造出精度更高的计算公式。这就是龙格—库塔方法的基本思想。4.1首先不妨在区间内仍取2个点，仿照（13）式用以下形式试一下（14）其中为待定系数，看看如何确定它们使（14）式的精度尽量高。为此我们分析局部截断误差，因为，所以（14）可以化为（15）其中在点作了Taylor展开。（15）式又可表为注意到中，，可见为使误差，只须令，，（16）待定系数满足（16）的（15）式称为2阶龙格—库塔公式。由于（16）式有4个未知数而只有3个方程，所以解不唯一。不难发现，若令，，即为改进的Euler公式。可以证明，在内只取2点的龙格—库塔公式精度最高为2阶。4.24阶龙格—库塔公式要进一步提高精度，必须取更多的点，如取4点构造如下形式的公式：（17）其中待定系数共13个，经过与推导2阶龙格—库塔公式类似、但更复杂的计算，得到使局部误差的11个方程。取既满足这些方程、又较简单的一组，可得（18）这就是常用的4阶龙格—库塔方法（简称RK方法）。§5线性多步法以上所介绍的各种数值解法都是单步法，这是因为它们在计算时，都只用到前一步的值，单步法的一般形式是（19）其中称为增量函数，例如Euler方法的增量函数为，改进Euler法的增量函数为如何通过较多地利用前面的已知信息，如，来构造高精度的算法计算，这就是多步法的基本思想。经常使用的是线性多步法。让我们试验一下，即利用计算的情况。从用数值积分方法离散化方程的（4）式出发，记，，式中被积函数用二节点，的插值公式得到（因，所以是外插。（20）此式在区间上积分可得于是得到（21）注意到插值公式（20）的误差项含因子，在区间上积分后得出，故公式（21）的局部截断误差为，精度比向前Euler公式提高1阶。若取可以用类似的方法推导公式，如对于有（22）其局部截断误差为。如果将上面代替被积函数用的插值公式由外插改为内插，可进一步减小误差。内插法用的是，取时得到的是梯形公式，取时可得（23）与（22）式相比，虽然其局部截断误差仍为，但因它的各项系数（绝对值）大为减小，误差还是小了。当然，（23）式右端的未知，需要如同向后Euler公式一样，用迭代或校正的办法处理。§6一阶微分方程组与高阶微分方程的数值解法6.1一阶微分方程组的数值解法设有一阶微分方程组的初值问题（24）若记，，，则初值问题（24）可写成如下向量形式（25）如果向量函数在区域：连续，且关于满足Lipschitz条件，即存在，使得对，，都有那么问题（25）在上存在唯一解。问题（25）与（1）形式上完全相同，故对初值问题（1）所建立的各种数值解法可全部用于求解问题（25）。6.2高阶微分方程的数值解法高阶微分方程的初值问题可以通过变量代换化为一阶微分方程组初值问题。设有阶常微分方程初值问题（26）引入新变量，问题（26）就化为一阶微分方程初值问题（27）然后用6.1中的数值方法求解问题（27），就可以得到问题（26）的数值解。最后需要指出的是，在化学工程及自动控制等领域中，所涉及的常微分方程组初值问题常常是所谓的“刚性”问题。具体地说，对一阶线性微分方程组（28）其中为阶方阵。若矩阵的特征值满足关系则称方程组（28）为刚性方程组或Stiff方程组，称数为刚性比。对刚性方程组，用前面所介绍的方法求解，都会遇到本质上的困难，这是由数值方法本身的稳定性限制所决定的。理论上的分析表明，求解刚性问题所选用的数值方法最好是对步长不作任何限制。§7Matlab解法7.1Matlab数值解7.1.1非刚性常微分方程的解法Matlab的工具箱提供了几个解非刚性常微分方程的功能函数，如ode45，ode23，ode113，其中ode45采用四五阶RK方法，是解非刚性常微分方程的首选方法，ode23采用二三阶RK方法，ode113采用的是多步法，效率一般比ode45高。Matlab的工具箱中没有Euler方法的功能函数。（=1\*ROMANI）对简单的一阶方程的初值问题改进的Euler公式为我们自己编写改进的Euler方法函数eulerpro.m如下：function[x,y]=eulerpro(fun,x0,xfinal,y0,n);ifnargin<5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;fori=1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i),y(i));y2=y(i)+h*feval(fun,x(i+1),y1);y(i+1)=(y1+y2)/2;end例1用改进的Euler方法求解，，解编写函数文件doty.m如下：functionf=doty(x,y);f=-2*y+2*x^2+2*x;在Matlab命令窗口输入：[x,y]=eulerpro('doty',0,0.5,1,10)即可求得数值解。（=2\*ROMANII）ode23，ode45，ode113的使用Matlab的函数形式是[t,y]=solver('F',tspan，y0)这里solver为ode45，ode23，ode113，输入参数F是用M文件定义的微分方程右端的函数。tspan=[t0,tfinal]是求解区间，y0是初值。例2用RK方法求解，，解同样地编写函数文件doty.m如下：functionf=doty(x,y);f=-2*y+2*x^2+2*x;在Matlab命令窗口输入：[x,y]=ode45('doty',0,0.5,1)即可求得数值解。7.1.2刚性常微分方程的解法Matlab的工具箱提供了几个解刚性常微分方程的功能函数，如ode15s，ode23s，ode23t，ode23tb，这些函数的使用同上述非刚性微分方程的功能函数。7.1.3高阶微分方程解法例3考虑初值问题解（=1\*romani）如果设，那么初值问题可以写成的形式，其中。（=2\*romanii）把一阶方程组写成接受两个参数和，返回一个列向量的M文件F.m：functiondy=F(t,y);dy=[y(2);y(3);3*y(3)+y(2)*y(1)];注意：尽管不一定用到参数和，M—文件必须接受此两参数。这里向量必须是列向量。（=3\*romaniii）用Matlab解决此问题的函数形式为[T,Y]=solver('F',tspan,y0)这里solver为ode45、ode23、ode113，输入参数F是用M文件定义的常微分方程组，tspan=[t0tfinal]是求解区间，y0是初值列向量。在Matlab命令窗口输入[T,Y]=ode45('F',[01],[0;1;-1])就得到上述常微分方程的数值解。这里Y和时刻T是一一对应的，Y(:,1)是初值问题的解，Y(:,2)是解的导数，Y(:,3)是解的二阶导数。例4求vanderPol方程的数值解，这里是一参数。解（=1\*romani）化成常微分方程组。设，则有（=2\*romanii）书写M文件（对于）vdp1.m:functiondy=vdp1(t,y);dy=[y(2);(1-y(1)^2)*y(2)-y(1)];（=3\*romaniii）调用Matlab函数。对于初值，解为[T,Y]=ode45('vdp1',[020],[2;0]);（=4\*romaniv）观察结果。利用图形输出解的结果：plot(T,Y(:,1),'-',T,Y(:,2),'--')title('SolutionofvanderPolEquation,mu=1');xlabel('timet');ylabel('solutiony');legend('y1','y2');例5vanderPol方程，（刚性）解(=1\*romani)书写M文件vdp1000.m:functiondy=vdp1000(t,y);dy=[y(2);1000*(1-y(1)^2)*y(2)-y(1)];（=2\*romanii）观察结果[t,y]=ode15s('vdp1000',[03000],[2;0]);plot(t,y(:,1),'o')title('SolutionofvanderPolEquation,mu=1000');xlabel('timet');ylabel('solutiony(:,1)');7.2常微分方程的解析解在Matlab中，符号运算工具箱提供了功能强大的求解常微分方程的符号运算命令dsolve。常微分方程在Matlab中按如下规定重新表达：符号D表示对变量的求导。Dy表示对变量y求一阶导数，当需要求变量的n阶导数时，用Dn表示，D4y表示对变量y求4阶导数。由此，常微分方程在Matlab中，将写成'D2y+2*Dy=y'。7.2.1求解常微分方程的通解无初边值条件的常微分方程的解就是该方程的通解。其使用格式为：dsolve('diff_equation')dsolve('diff_equation'，'var')式中diff_equation为待解的常微分方程，第1种格式将以变量t为自变量进行求解，第2种格式则需定义自变量var。例6试解常微分方程解编写程序如下：symsxydiff_equ='x^2+y+(x-2*y)*Dy=0';dsolve(diff_equ,'x')7.2.2求解常微分方程的初边值问题求解带有初边值条件的常微分方程的使用格式为：dsolve('diff_equation'，'condition1，condition2，…'，'var')其中condition1，condition2，…即为微分方程的初边值条件。例7试求微分方程，的解。解编写程序如下：y=dsolve('D3y-D2y=x','y(1)=8,Dy(1)=7,D2y(2)=4','x')7.2.3求解常微分方程组求解常微分方程组的命令格式为：dsolve('diff_equ1，diff_equ2，…'，'var')dsolve('diff_equ1，diff_equ2，…'，'condition1，condition2，…'，'var')第1种格式用于求解方程组的通解，第2种格式可以加上初边值条件，用于具体求解。例8试求常微分方程组：的通解和在初边值条件为的解。解编写程序如下：clc,clearequ1='D2f+3*g=sin(x)';equ2='Dg+Df=cos(x)';[general_f,general_g]=dsolve(equ1,equ2,'x')[f,g]=dsolve(equ1,equ2,'Df(2)=0,f(3)=3,g(5)=1','x')7.2.4求解线性常微分方程组（=1\*romani）一阶齐次线性微分方程组，，这里的’表示对求导数。是它的基解矩阵。，的解为。例9试解初值问题，解编写程序如下：symsta=[2,1,3;0,2,-1;0,0,2];x0=[1;2;1];x=expm(a*t)*x0（=2\*romanii）非齐次线性方程组由参数变易法可求得初值问题，的解为.例10试解初值问题，。解编写程序如下：clc,clearsymstsa=[1,0,0;2,1,-2;3,2,1];ft=[0;0;exp(t)*cos(2*t)];x0=[0;1;1];x=expm(a*t)*x0+int(expm(a*(t-s))*subs(ft,s),s,0,t);x=simple(x)习题十五1.用欧拉方法和龙格—库塔方法求微分方程数值解，画出解的图形，对结果进行分析比较。（=1\*romani）（Bessel方程，令，精确解。（=2\*romanii），幂级数解2.一只小船渡过宽为的河流，目标是起点正对着的另一岸点。已知河水流速与船在静水中的速度之比为。（=1\*romani）建立小船航线的方程，求其解析解。（=2\*romanii）设m，m/s，m/s，用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较。第二十三章现代优化算法简介§1现代优化算法简介现代优化算法是80年代初兴起的启发式算法。这些算法包括禁忌搜索（tabusearch），模拟退火（simulatedannealing），遗传算法（geneticalgorithms），人工神经网络（neuralnetworks）。它们主要用于解决大量的实际应用问题。目前，这些算法在理论和实际应用方面得到了较大的发展。无论这些算法是怎样产生的，它们有一个共同的目标－求NP-hard组合优化问题的全局最优解。虽然有这些目标，但NP-hard理论限制它们只能以启发式的算法去求解实际问题。启发式算法包含的算法很多，例如解决复杂优化问题的蚁群算法（AntColonyAlgorithms）。有些启发式算法是根据实际问题而产生的，如解空间分解、解空间的限制等；另一类算法是集成算法，这些算法是诸多启发式算法的合成。现代优化算法解决组合优化问题，如TSP（TravelingSalesmanProblem）问题，QAP（QuadraticAssignmentProblem）问题，JSP（Job-shopSchedulingProblem）问题等效果很好。本章我们只介绍模拟退火算法，初步介绍一下蚁群算法，其它优化算法可以参看相关的参考资料。§2模拟退火算法2.1算法简介模拟退火算法得益于材料的统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下，粒子的能量较高，可以自由运动和重新排列。在低温条件下，粒子能量较低。如果从高温开始，非常缓慢地降温（这个过程被称为退火），粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。如果用粒子的能量定义材料的状态，Metropolis算法用一个简单的数学模型描述了退火过程。假设材料在状态之下的能量为，那么材料在温度时从状态进入状态就遵循如下规律：（1）如果，接受该状态被转换。（2）如果，则状态转换以如下概率被接受：其中是物理学中的波尔兹曼常数，是材料温度。在某一个特定温度下，进行了充分的转换之后，材料将达到热平衡。这时材料处于状态的概率满足波尔兹曼分布：其中表示材料当前状态的随机变量，表示状态空间集合。显然其中表示集合中状态的数量。这表明所有状态在高温下具有相同的概率。而当温度下降时，其中且。上式表明当温度降至很低时，材料会以很大概率进入最小能量状态。假定我们要解决的问题是一个寻找最小值的优化问题。将物理学中模拟退火的思想应用于优化问题就可以得到模拟退火寻优方法。考虑这样一个组合优化问题：优化函数为，其中，它表示优化问题的一个可行解，，表示函数的定义域。表示的一个邻域集合。首先给定一个初始温度和该优化问题的一个初始解，并由生成下一个解，是否接受作为一个新解依赖于下面概率：换句话说，如果生成的解的函数值比前一个解的函数值更小，则接受作为一个新解。否则以概率接受作为一个新解。泛泛地说，对于某一个温度和该优化问题的一个解，可以生成。接受作为下一个新解的概率为：（1）在温度下，经过很多次的转移之后，降低温度，得到。在下重复上述过程。因此整个优化过程就是不断寻找新解和缓慢降温的交替过程。最终的解是对该问题寻优的结果。我们注意到，在每个下，所得到的一个新状态完全依赖于前一个状态，可以和前面的状态无关，因此这是一个马尔可夫过程。使用马尔可夫过程对上述模拟退火的步骤进行分析，结果表明：从任何一个状态生成的概率，在中是均匀分布的，且新状态被接受的概率满足式（1），那么经过有限次的转换，在温度下的平衡态的分布由下式给出：（2）当温度降为0时，的分布为：并且这说明如果温度下降十分缓慢，而在每个温度都有足够多次的状态转移，使之在每一个温度下达到热平衡，则全局最优解将以概率1被找到。因此可以说模拟退火算法可以找到全局最优解。在模拟退火算法中应注意以下问题：（1）理论上，降温过程要足够缓慢，要使得在每一温度下达到热平衡。但在计算机实现中，如果降温速度过缓，所得到的解的性能会较为令人满意，但是算法会太慢，相对于简单的搜索算法不具有明显优势。如果降温速度过快，很可能最终得不到全局最优解。因此使用时要综合考虑解的性能和算法速度，在两者之间采取一种折衷。（2）要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续次的转换过程没有使状态发生变化时结束该温度下的状态转换。最终温度的确定可以提前定为一个较小的值，或连续几个温度下转换过程没有使状态发生变化算法就结束。（3）选择初始温度和确定某个可行解的邻域的方法也要恰当。2.2应用举例例已知敌方100个目标的经度、纬度如下：经度纬度经度纬度经度纬度经度纬度53.712115.304651.17580.032246.325328.275330.33136.934856.543221.418810.819816.252922.789123.104510.158412.481920.105015.45621.94510.205726.495122.122131.48478.964026.241818.176044.035613.540128.983625.987938.472220.173128.269429.001132.19105.869936.486329.72840.971828.14778.958624.663516.561823.614310.559715.117850.211110.29448.15199.532522.107518.55690.121518.872648.207716.888931.949917.63090.77320.465647.413423.778341.86713.566743.54743.906153.352426.725630.816513.459527.71335.070623.92227.630651.961222.851112.793815.73074.95688.366921.505124.090915.254827.21116.20705.144249.243016.704417.116820.035434.168822.75719.44023.920011.581214.567752.11810.40889.555911.421924.45096.563426.721328.566737.584816.847435.66199.933324.46543.16440.77756.957614.470313.636819.866015.12243.16164.242818.524514.359858.684927.148539.516816.937156.508913.709052.521115.795738.43008.464851.818123.01598.998323.644050.115623.781613.79091.951034.057423.396023.06248.431919.98575.790240.880114.297858.828914.522918.66356.743652.842327.288039.949429.511447.509924.066410.112127.266228.781227.66598.083127.67059.155614.130453.79890.219933.64900.39801.349616.835949.98166.082819.363517.662236.954523.026515.732019.569711.511817.388444.039816.263539.713928.42036.990923.180438.339219.995024.654319.605736.998024.39924.15913.185340.140020.303023.98769.403041.108427.7149我方有一个基地，经度和纬度为（70,40）。假设我方飞机的速度为1000公里/小时。我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。这是一个旅行商问题。我们依次给基地编号为1，敌方目标依次编号为2，3，…，101，最后我方基地再重复编号为102（这样便于程序中计算）。距离矩阵，其中表示表示两点的距离，，这里为实对称矩阵。则问题是求一个从点1出发，走遍所有中间点，到达点102的一个最短路径。上面问题中给定的是地理坐标（经度和纬度），我们必须求两点间的实际距离。设两点的地理坐标分别为，，过两点的大圆的劣弧长即为两点的实际距离。以地心为坐标原点，以赤道平面为平面，以0度经线圈所在的平面为平面建立三维直角坐标系。则两点的直角坐标分别为：其中为地球半径。两点的实际距离，化简得。求解的模拟退火算法描述如下：（1）解空间解空间可表为{}的所有固定起点和终点的循环排列集合，即其中每一个循环排列表示侦察100个目标的一个回路，表示在第次侦察点，初始解可选为，本文中我们使用MonteCarlo方法求得一个较好的初始解。（2）目标函数此时的目标函数为侦察所有目标的路径长度或称代价函数。我们要求而一次迭代由下列三步构成：（3）新解的产生=1\*GB3①2变换法任选序号（）交换与之间的顺序，此时的新路径为：②3变换法任选序号和，将和之间的路径插到之后，对应的新路径为（设）（4）代价函数差对于2变换法，路径差可表示为（5）接受准则如果，则接受新的路径。否则，以概率接受新的路径，即若大于0到1之间的随机数则接受。（6）降温利用选定的降温系数进行降温即：，得到新的温度，这里我们取。（7）结束条件用选定的终止温度，判断退火过程是否结束。若，算法结束，输出当前状态。我们编写如下的matlab程序如下：clc,clearloadsj.txt%加载敌方100个目标的数据，数据按照表格中的位置保存在纯文本文件sj.txt中x=sj(:,1:2:8);x=x(:);y=sj(:,2:2:8);y=y(:);sj=[xy];d1=[70,40];sj=[d1;sj;d1];sj=sj*pi/180;%距离矩阵dd=zeros(102);fori=1:101forj=i+1:102temp=cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2));d(i,j)=6370*acos(temp);endendd=d+d';S0=[];Sum=inf;rand('state',sum(clock));forj=1:1000S=[11+randperm(100),102];temp=0;fori=1:101temp=temp+d(S(i),S(i+1));endiftemp<SumS0=S;Sum=temp;endende=0.1^30;L=20000;at=0.999;T=1;%退火过程fork=1:L%产生新解c=2+floor(100*rand(1,2));c=sort(c);c1=c(1);c2=c(2);%计算代价函数值df=d(S0(c1-1),S0(c2))+d(S0(c1),S0(c2+1))-d(S0(c1-1),S0(c1))-d(S0(c2),S0(c2+1));%接受准则ifdf<0S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];Sum=Sum+df;elseifexp(-df/T)>rand(1)S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];Sum=Sum+df;endT=T*at;ifT<ebreak;endend%输出巡航路径及路径长度S0,Sum计算结果为44小时左右。其中的一个巡航路径如下图所示：§3蚁群算法3.1蚁群算法简介蚁群是自然界中常见的一种生物，人们对蚂蚁的关注大都是因为“蚁群搬家，天要下雨”之类的民谚。然而随着近代仿生学的发展，这种似乎微不足道的小东西越来越多地受到学者们地关注。1991年意大利学者M.Dorigo等人首先提出了蚁群算法，人们开始了对蚁群的研究：相对弱小，功能并不强大的个体是如何完成复杂的工作的（如寻找到食物的最佳路径并返回等）。在此基础上一种很好的优化算法逐步发展起来。蚁群算法的特点是模拟自然界中蚂蚁的群体行为。科学家发现，蚁群总是能够发现从蚁巢到食物源的最短路径。经研究发现，蚂蚁在行走过的路上留下一种挥发性的激素，蚂蚁就是通过这种激素进行信息交流。蚂蚁趋向于走激素积累较多的路径。找到最短路径的蚂蚁总是最早返回巢穴，从而在路上留下了较多的激素。由于最短路径上积累了较多的激素，选择这条路径的蚂蚁就会越来越多，到最后所有的蚂蚁都会趋向于选择这条最短路径。基于蚂蚁这种行为而提出的蚁群算法具有群体合作，正反馈选择，并行计算等三大特点，并且可以根据需要为人工蚁加入前瞻、回溯等自然蚁所没有的特点。在使用蚁群算法求解现实问题时，先生成具有一定数量蚂蚁的蚁群，让每一只蚂蚁建立一个解或解的一部分，每只人工蚁从问题的初始状态出发，根据“激素”浓度来选择下一个要转移到的状态，直到建立起一个解，每只蚂蚁根据所找到的解的好坏程度在所经过的状态上释放与解的质量成正比例的“激素”。之后，每只蚂蚁又开始新的求解过程，直到寻找到满意解。为避免停滞现象，引入了激素更新机制。3.2解决TSP问题的蚁群算法描述现以TSP问题的求解为例说明蚁群系统模型。首先引进如下记号：为城市的个数；为蚁群中蚂蚁的数量；为两城市和之间距离；为时刻位于城市的蚂蚁的个数，；为时刻边弧的轨迹强度（即连线上残留的信息量），且设（为常数），，；为时刻边弧的能见度，反映由城市转移到城市的期望程度。根据上述原理，蚂蚁在运动过程中根据各条路径上的信息量决定转移方向。与真实蚁群系统不同，人工蚁群系统具有一定的记忆功能。随着时间的推移，以前留下的信息逐渐消逝，经个时刻，蚂蚁完成一次循环，各路径上信息量要作调整。由此得到下述的人工蚁群系统模型：1）设人工蚁群在并行地搜索TSP的解，并通过一种信息素做媒介相互通信，在每个结点上且和该结点相连的边上以信息素量做搜索下一结点的试探依据，直到找到一个TSP问题的可行解。2）在时刻人工蚁由位置转移至位置的转移概率为（3）其中参数为轨迹的相对重要性（）；为能见度的相对重要性（）；为可行点集，即蚂蚁下一步允许选择的城市。分别反映了蚂蚁在运动过程中所积累的信息及启发式因子在蚂蚁选择路径中所起的不同作用。3）当个人工蚁按（3）式找到了可行解，则将各边的信息量用下式修改。即调整信息量的轨迹强度更新方程为，（4）其中为第只蚂蚁在本次循环中留在路径上的信息量；为本次循环中路径上的信息量的增量；参数为轨迹的持久性；为轨迹衰减度，表示信息消逝程度。对上述系统模型，采用人工蚁群方法求解的算法步骤可归结为：step1：（为迭代步数或搜索次数）；各和的初始化；将个蚂蚁置于个顶点上。step2：将各蚂蚁的初始出发点置于当前解集中；对每个蚂蚁（）按概率转移至下一顶点；将顶点置于当前解集。step3：计算各蚂蚁的目标函数值，记录当前的最好解。step4：按更新方程修改轨迹强度。step5：，若预定的迭代次数且无退化行为（即找到的都是相同解），则转step2。若为了简化计算，增加处理较大规模的TSP问题的能力，则可将（4）式修改为：

，其中此处为本次最优路线上的边集。

3.3人工蚁群算法性能的讨论人工蚁群算法是一种基于种群的进化算法。作为一个新兴的研究领域，虽它还远未像GA、SA等算法那样形成系统的分析方法和坚实的数学基础，但目前已有一些基本结果。在M.Dorigo三种不同的模型中，循环路径上信息量的增量不同。1）Ant-quantitysystem模型中，2）在Ant-densitysystem模型中，3）在Ant-cyclesystem模型中，其中是反映蚂蚁所留轨迹数量的常数，表示第只蚂蚁在本次循环中所走路径的长度；且时，，。算法中模型1）、2）利用的是局部信息，模型3）利用的是整体信息。人工蚁群算法中，等参数对算法性能也有很大的影响。值的大小表明留在每个结点上的信息量受重视的程度，值越大，蚂蚁选择以前选过的点的可能性越大，但过大会使搜索过早陷于局部极小点；的大小表明启发式信息受重视的程度；值会影响算法的收敛速度，过大会使算法收敛于局部极小值，过小又会影响算法的收敛速度，随问题规模的增大的值也需要随之变化；蚂蚁的数目越多，算法的全局搜索能力越强，但数目加大将使算法的收敛速度减慢。第二十四章时间序列模型时间序列是按时间顺序排列的、随时间变化且相互关联的数据序列。分析时间序列的方法构成数据分析的一个重要领域，即时间序列分析。时间序列根据所研究的依据不同，可有不同的分类。1．按所研究的对象的多少分，有一元时间序列和多元时间序列。2．按时间的连续性可将时间序列分为离散时间序列和连续时间序列两种。3．按序列的统计特性分，有平稳时间序列和非平稳时间序列。如果一个时间序列的概率分布与时间无关，则称该序列为严格的（狭义的）平稳时间序列。如果序列的一、二阶矩存在，而且对任意时刻满足：（1）均值为常数（2）协方差为时间间隔的函数。则称该序列为宽平稳时间序列，也叫广义平稳时间序列。我们以后所研究的时间序列主要是宽平稳时间序列。4．按时间序列的分布规律来分，有高斯型时间序列和非高斯型时间序列。§1确定性时间序列分析方法概述时间序列预测技术就是通过对预测目标自身时间序列的处理，来研究其变化趋势的。一个时间序列往往是以下几类变化形式的叠加或耦合。（1）长期趋势变动。它是指时间序列朝着一定的方向持续上升或下降，或停留在某一水平上的倾向，它反映了客观事物的主要变化趋势。（2）季节变动。（3）循环变动。通常是指周期为一年以上，由非季节因素引起的涨落起伏波形相似的波动。（4）不规则变动。通常它分为突然变动和随机变动。通常用表示长期趋势项，表示季节变动趋势项，表示循环变动趋势项，表示随机干扰项。常见的确定性时间序列模型有以下几种类型：（1）加法模型（2）乘法模型（3）混合模型其中是观测目标的观测记录，，。如果在预测时间范围以内，无突然变动且随机变动的方差较小，并且有理由认为过去和现在的演变趋势将继续发展到未来时，可用一些经验方法进行预测，具体方法如下：1.1移动平均法设观测序列为，取移动平均的项数。一次移动平均值计算公式为：二次移动平均值计算公式为：当预测目标的基本趋势是在某一水平上下波动时，可用一次移动平均方法建立预测模型：，，其预测标准误差为：，最近期序列值的平均值作为未来各期的预测结果。一般取值范围：。当历史序列的基本趋势变化不大且序列中随机变动成分较多时，的取值应较大一些。否则的取值应小一些。在有确定的季节变动周期的资料中，移动平均的项数应取周期长度。选择最佳值的一个有效方法是，比较若干模型的预测误差。均方预测误差最小者为好。当预测目标的基本趋势与某一线性模型相吻合时，常用二次移动平均法，但序列同时存在线性趋势与周期波动时，可用趋势移动平均法建立预测模型：，其中，。例1某企业1月～11月份的销售收入时间序列如下表所示。取，试用简单一次滑动平均法预测第12月份的销售收入，并计算预测的标准误差。月份123456销售收入533.8574.6606.9649.8705.1772.0月份789101112销售收入816.4892.7963.91015.11102.7解：首先计算出，由于，，则，预测的标准误差为计算的Matlab程序如下：y=[533.8574.6606.9649.8705.1772.0...816.4892.7963.91015.11102.7];temp=cumsum(y);mt=(temp(4:11)-[0temp(1:7)])/4y12=mt(end)ythat=mt(1:end-1);fangcha=mean((y(5:11)-ythat).^2);sigma=sqrt(fangcha)1.2指数平滑法一次移动平均实际上认为最近期数据对未来值影响相同，都加权；而期以前的数据对未来值没有影响，加权为0。但是，二次及更高次移动平均数的权数却不是，且次数越高，权数的结构越复杂，但永远保持对称的权数，即两端项权数小，中间项权数大，不符合一般系统的动态性。一般说来历史数据对未来值的影响是随时间间隔的增长而递减的。所以，更切合实际的方法应是对各期观测值依时间顺序进行加权平均作为预测值。指数平滑法可满足这一要求，而且具有简单的递推形式。设观测序列为，为加权系数，，一次指数平滑公式为：（1）假定历史序列无限长，则有（2）（2）式表明是全部历史数据的加权平均，加权系数分别为；显然有由于加权系数序列呈指数函数衰减，加权平均又能消除或减弱随机干扰的影响，所以（2）称为一次指数平滑，类似地，二次指数平滑公式为：（3）同理，三次指数平滑公式为：（4）一般次指数平滑公式为：