常微分方程求解的高阶方法毕业论文.doc 29页

'　常微分方程求解的高阶方法毕业论文目录第一章前言11.1案例引入微分方程概念11.2微分方程的基本概念11.2.1微分方程及微分方程的阶11.2.2微分方程的解、通解与特解11.2.3微分方程的初值条件及其提法21.2.4微分方程的解的几何意义.21.3从解析方法到数值方法概述31.4常温分方程的离散化4第二章数值解法公共程序模块分析5第三章欧拉（Euler）方法73.1Euler方法思想73.2Euler方法的误差估计83.3改进的Euler方法83.3.1梯形公式83.3.2改进Euler法9第四章休恩方法104.1休恩方法思想104.2休恩方法的步长和误差10第五章泰勒级数法115.1泰勒定理115.2N次泰勒方法12第六章龙格-库塔（Runge—Kutta法）136.1龙格-库塔（Runge—Kutta）方法基本思想136.2阶龙格-库塔（Runge—Kutta）方法公式14第七章预报-校正方法157.1Milne-Simpon方法167.2误差估计于校正167.3正确的步长17 第八章一阶微分方程组与高阶微分方程的数值解法178.1一阶微分方程组的数值解法178.2高阶微分方程的数值解法18第九章常微分方程模型数值解法在数学建模中的应用199.1耐用消费新产品的销售规律模型199.1.1问题的提出199.1.2模型的构建199.1.3模型的求解209.2司机饮酒驾车防避模型的数值解法219.2.1模型假设229.2.2模型建立229.2.3模型求解249.2.4模型评价259.2.5诚恳建议259.2.6模型推广26主要参考文献26致谢27 第一章前言1.1案例引入微分方程概念在科技、工程、经济管理、生态、生态、刑侦等各个领域微分方程有着广泛的应用。我们看一实例。案例：一次谋杀案，在某天下午四点发现尸体，尸体的体温为30℃，假设当时屋内空间的温度保护20℃不变，现判断谋杀是何时发生的？解决此问题首先必须要从尸体温度的变化寻求关系式，这就需要知道物理学中的加热与冷却规律。物理学家牛顿（Newton）曾提出，一块热的物体，其温度下降的速度是与它自身温度的差值成正比。同样，一块冷的物体，其温度上升的速度是与他自身温度同外界温度的差值成正比。据此我们可找到温度与时间之间的函数关系式，这事实上就是一个微分方程的建立问题。再如传染病传染问题（人口增长模型问题）也要用到微分方程的知识。通过求解微分方程，可以得到所需求的函数。1.2微分方程的基本概念1.2.1微分方程及微分方程的阶含未知函数的导数(或微分)的方程称为微分方程；未知函数是一元函数的微分方程，称为常微分方程；未知函数是多元函数的微分方程，称为偏微分方；(1.1)和(1.5)式均是微分方程.微分方程中未知函数的导数的最高阶数，称为微分方程的阶.微分方程(1.1)是一阶的，微分方程(1.2)是二阶的.1.2.2微分方程的解、通解与特解能使微分方程成为恒等式的函数，称为微分方程的解.例如和都是的解. 又如和都是的解.如果微分方程的解中含任意常数,且独立的(即不可合并而使个数减少的)任意常数的个数与微分方程的阶数相同，这样的解为微分方程的通解.不包含任意常数的解为微分方程特解.1.2.3微分方程的初值条件及其提法用以确定微分方程解中任意常数的特定条件，称为微分方程的初值条件.初值条件的提法：当x=x0时，y=y0，1.2.4微分方程的解的几何意义.微分方程的解的图形称为微分方程的积分曲线.通解的图形是一族积分曲线，称为微分方程的积分曲线族.微分方程的某个特解的图形就是积分曲线族中满足给定初值条件的某一特定的积分曲线. 所以函数是所给微分方程(1.3)的解.又因为这个解中含有两个独立的任意常数，任意常数的个数与微分方程(1.3)的阶数相同，所以它是该方程的通解.1.3从解析方法到数值方法概述求解常微分方程的解析方法很多，像变量分离法，积分因子法，遗憾的是实际上得到的大部分常微分方程都不能使用这些理论上的方法。数值求解微分方程的方法基于有限维近似，这个过程称为离散化，我们将用代数方程代替微分方程，用代数方程的解近似微分方程的解，对初值问题来说，近似解的值是在求解区间上一步步地产生的，因此求解常微分方程的数值方法也称为离散变量法，在由一个离散点的值计算下一个点的值时，一般会产生一定的误差，这样新的近似解将落在常微分方程的另一个解上，而这个解与开始所求的解是不同的，解的稳定性决定了这类误差将随时间的增大而放大或缩小。 1.4常温分方程的离散化下面主要讨论一阶常微分方程的初值问题，其一般形式是在下面的讨论我们总假定函数f(x,y)连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L，使得这样，由常微分方程理论知，初值问题(1)的解必定存在唯一。所谓数值解法，就是求问题（1.5）的解y(x)在若干点处的近似值的方法，称为问题（1.5）的数值解，称为由到的步长。今后如无特别说明，我们总取步长为常量h。建立数值解法，首先要将微分方程离散化，一般采用以下几种方法：（i）用差商近似导数若用向前差商代替代入（1.5）中的微分方程，则得化简得如果用的近似值代入上式右端，所得结果作为的近似值，记为，则有这样，问题（1.5）的近似值可通过求解下述问题得到，按式（1.7）由初值可逐次算出。式（1,7）是个离散化的问题，称为差分方程初值问题。需要说明的是，用不同的差商近似导数，将得到不同的计算公式。 （ii）用数值积分方法将问题（1.5）的解表成积分形式，用数值积分方法离散化。例如，对微分方程两端积分，得右边的积分用矩形公式或梯形公式计算。（iii）Taylor多项式近似将函数在处展开，取一次Taylor多项式近似，则得再将的近似值代入上式右端，所得结果作为的近似值，得到离散化的计算公式以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式。其中的Taylor展开法，不仅可以得到求数值解的公式，而且容易估计截断误差。第二章数值解法公共程序模块分析编程选择:由于并不需要采用STL等泛型程序设计的方法，采用C++并不会比采用C减少太多代码，况且这里的实际代码比较简单，所以为了减少系统的开销，采用TubroC来实验。编程风格:按照常微分方程数值解三个基本步骤：将问题离散化；建立递推格式；按步进法计算，所以求微分方程的数值解的算法框架都是相同的，不同的是所使用的递推形式不同，则可以用公共子程序来代替，对不同的方法的计算结果用统一的格式来显示，同时也可以比较不同方法的精确度[4]。公共程序模块如下:这里为了良好地比较，选用可求解析解的一阶常微分方程作为讨论:(2.1)其解析式为/*Filename:numerical_base.c*/#include#include #include#defineMAX100intreal=1;doubleX[MAX],Y[MAX],Z[MAX],F[MAX],G[MAX],CY[MAX],E[MAX];/******以下代码根据待求解的对象的特殊性进行赋值和在main()中选取*****/doublefunc(doublex,doubley)//计算各离散点处导数值{returny-2.0*x/y;}doubleexact_value(doublex)//计算各离散点解析解以测数值解精度{returnsqrt(2.0*x+1.0);}/*********************changingpartend**************************/voidcal_error()//计算误差值以对各种方法进行比较{inti;for(i=0;i<=real;i++){CY[i]=exact_value(X[i]);E[i]=fabs(CY[i]-Y[i]);}return;}voidshowtable_s()//微分方程组输出时用{//内容与showtable()类似}//输出各离散点处的X值，Y值，导数值，精确值，误差值//分别对应于X[k]，Y[k]，F[k]，CY[k]，E[k]voidshowtable()//优化输出显示{inti,j;printf("n");for(i=0;i<=78;i++)printf("=");printf("n");printf("%5s%8s%15s%15s%18s%14s","k","X[k]","F[k]","Y[k]","CY[k]","E[k]");for(i=0;i<=78;i++)printf("-");for(i=0;i<=real;i++){printf("n%5d",i);printf("%11.4E%17.10E%17.10E%17.10E%11.4E",X[i],F[i], Y[i],CY[i],E[i]);}printf("n");for(i=0;i<=78;i++)printf("=");getchar();}对程序中的诸名称解释如下：符号常量MAX表示允许取的离散点的个数的最大值，用于初始化数组；变量real表示实际取的离散点的个数；数组X,Y,F用于存放各,,,数组CY,E用于存放离散点处的精确解的值和的整体误差变量a,b,x0,y0用于表示求解区间的左右端点和给定的初始点；函数exact_value()用于求解问题的精确解；函数cal_error()用于计算各离散点的截断误差；函数showtable()用于在屏幕上显示计算结果；第三章欧拉（Euler）方法3.1Euler方法思想Euler方法就是用差分方程初值问题（1.7）的解来近似微分方程初值问题（1.5）的解，即由公式（1.7）依次算出的近似值。这组公式求问题（1.5）的数值解称为向前Euler公式。如果在微分方程离散化时，用向后差商代替导数，即，则得计算公式用这组公式求问题（1.5）的数值解称为向后Euler公式。向后Euler法与Euler法形式上相似，但实际计算时却复杂得多。向前Euler公式是显式的，可直接求解。向后Euler公式的右端含有，因此是隐式公式，一般要用迭代法求解，迭代公式通常为 3.2Euler方法的误差估计对于向前Euler公式（1.7）我们看到，当时公式右端的都是近似的，所以用它计算的会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的，所谓局部截断误差。假定用（1.7）式时右端的没有误差，即，那么由此算出局部截断误差指的是，按（3,3）式计算由到这一步的计算值与精确值之差。为了估计它，由Taylor展开得到的精确值是（3.3）、（3.4）两式相减（注意到）得即局部截断误差是阶的，而数值算法的精度定义为：若一种算法的局部截断误差为，则称该算法具有p阶精度。显然p越大，方法的精度越高。式（3.5）说明，向前Euler方法是一阶方法，因此它的精度不高。3.3改进的Euler方法3.3.1梯形公式利用数值积分方法将微分方程离散化时，若用梯形公式计算式(1.8)中之右端积分，即 这就是求解初值问题（1.5）的梯形公式。直观上容易看出，用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法。梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为由于函数关于y满足Lipschitz条件，容易看出其中L为Lipschitz常数。因此，当时，迭代收敛。但这样做计算量最大。如果实际计算时精度要求不太高，用公式（3.6）求解时，每步可以只迭代一次，由此导出一种新的方法—改进Euler法。3.3.2改进Euler法按式（3.1）计算问题（1.5）的数值解时，如果每步只迭代一次，相当于将Euler公式与梯形公式结合使用：先用Euler公式求的一个初步近似值，称为预测值，然后用梯形公式校正求得近似值，即式（3.7）称为由Euler公式和梯形公式得到的预测—校正系统，也叫改进Euler法。为便于编制程序上机，式（3.7）常改写成 改进Euler法是二阶方法。第四章休恩方法4.1休恩方法思想休恩（Heun）其实也是欧拉方法的改进型，因此也叫做欧拉改进法。该方法引入一种新的思路，来构造求解上的初值问题，（4.1）要得到解，可以用微积分基本定理，在上对积分得（4.2）其中的不定积分为待定求函数。对求解方程（4.2）（4.3）然后可用数值积分方法逼近（4.3）中的定积分，如果采用步长的梯形公式，则结果为（4.4）注意公式（4.4）的右端包含了待定值。要继续求解，需要的一个估计值，欧拉方法的解能够满足这一目的，将它代入（4.4）后，得到求解的公式，称为休恩（Heun）方法：（4.5）重复这个过程，得到逼近解曲线的一系列点，在每一步中都用欧拉方法作为预报，然后用梯形公式进行校正，得到最终的值。休恩方法的一般步骤为：(4.6)4.2休恩方法的步长和误差积分公式（4.3）中梯形公式的误差项为 （4.7）如果每步中的误差仅由式（4.7）给出，则步后休恩方法的累积误差将是（4.8）所以休恩方法的精确度可以这样计算：设是初值问题（4.1）的解，如果且是有欧拉方法计算的近似值序列，则（4.9）（4.10）其中，特别地，区间终点处的最终全局误差满足，（4.11）这就是说，步长如果减小为原来的（为整数），则可期望最终全局误差将降至大约。第五章泰勒级数法泰勒级数法有着广泛的应用，并且是比较求解初值问题的各种不同数值方法的标准，它可设计为具有任意指定的精度。5.1泰勒定理设，且在不动点处有次泰勒级数展开：（5.1）其中， （5.2）表示函数关于的次全导数。求导公式可以递归得计算：(5.3)并且一般有，（5.4）其中为导数算子区间上的初值的近似数值解可由各子区间上的公式（5.1）来推导。5.2N次泰勒方法次泰勒方法的一般步骤为：（5.5）其中在个步有，。次泰勒方法的最终全局误差是阶的，因此可选择所需大小的，使得误差足够小。如果固定，则理论上可以推导出步长，使之满足任意想要的最终全局误差。然而在实际运算中，通常用和计算两个近似结果集，然后比较其结果。次泰勒方法的精度：设为初值问题的解,如果，并且是休恩方法产生的一个近似值序列，则 （5.6）特别得，区间重点处的最终全局误差满足，（5.7）这就是说，步长如果减小为原来的，则可期望最终全局误差将降至大约，为整数。第六章龙格-库塔（Runge—Kutta法）6.1龙格-库塔（Runge—Kutta）方法基本思想回到Euler方法的基本思想—用差商代替导数—上来。实际上，按照微分中值定理应有注意到方程就有不妨记，称为区间上的平均斜率。可是给出一种斜率,(6.1)式就对应地导出一种算法。向前Euler公式简单地取为，精度自然很低。改进的Euler公式可理解为取，的平均值，其中，这种处理提高了精度。如上分析启示我们，在区间内多取几个点，将它们的斜率加权平均作为，就有可能构造出精度更高的计算公式。这就是龙格—库塔方法的基本思想。首先不妨在区间内仍取2个点，仿照（6.1）式用以下形式试一下 其中，，，为待定系数，看看如何确定它们使（6.2）式的精度尽量高。为此我们分析局部截断误差，因为，所以（6.2）可以化为其中在点作了Taylor展开。（6.3）式又可表为注意到中，可见为使误差，只须令待定系数满足（6.4）的（6.3）式称为2阶龙格—库塔公式。由于（6.4）式有4个未知数而只有3个方程，所以解不唯一。不难发现，若令即为改进的Euler公式。可以证明，在内只取2点的龙格—库塔公式精度最高为2阶。6.2阶龙格-库塔（Runge—Kutta）方法公式要进一步提高精度，必须取更多的点，如取4点构造如下形式的公式： 其中待定系数，，共13个，经过与推导2阶龙格—库塔公式类似、但更复杂的计算，得到使局部误差的11个方程。取既满足这些方程、又较简单的一组可得这就是常用的4阶龙格—库塔方法（简称RK方法）。第七章预报-校正方法欧拉方法、休恩方法、泰勒方法以及龙格—库塔方法都称为单步方法，因为他们只利用前一个点的信息计算下一个点，即计算时只使用了初始点。一般地，只有用来计算。当计算下一个点，就可以利用几个已计算出的点来计算下一个点。多步法的一个优点是，可以确定它的局部截断误差，并可以包含一个校正项，用于在每一步计算中改变解的精确度。该方法还可以确定步长是否小到能得到的精确值，同时又大到能够免除不必要的和费时的计算。使用预估子和校正子的组合在每一步只与要进行两次函数求值。对于多步法，常用的预估—校正公式有Milne-Hamning方法、Milne-Simpon方法、4阶隐预估—校正公式等。下面，我们向大家介绍Milne-Simpon方法。 7.1Milne-Simpon方法Milne-Simpon方法的预估子基于在区间上对的积分：（7.1）预报子使用的基于点，，和的拉格朗日多项式逼近，在区间上对他积分，得到米尔恩预估子：（7.2）校正子的推导类似。此时值已知，基于点，和新点构造的新点的拉格朗日多项式，然后在区间上对该多项式积分，结果得Simpon公式：（7.3）7.2误差估计于校正计算预估子和校正子的数值积分公式的误差项都是的，公式（7.2）和（7.3）的局部截断误差为（预估子的截断误差）（7.4）（校正子的截断误差）（7,5）设足够小，使得在区间上近乎为常数，则可消去式（7.4）和式（7.5）中的5阶导数项，结果为（7.6）公式（7.6）给出的预估子误差聚集基于两个计算值和，而没有使用高阶导数，可用它来改进预估值。假设每步中预估和校正值的差缓慢变化，则在式（7.6）中可用和分解替代和，得到如下的修正为：（7.7）在校正过程中用该修正值替代，公式（2.44）变为：（7.8） 因此，改进（修正）的Milne-Simpon方法为（7.9）（7.10）（7.11）由（7.9）和（7.10）知预估—校正方法是休恩方法的一个改进。7.3正确的步长用预估—校正方法在大区间上求解初值问题时，有时也出现问题。如果，而步长过大，则预计估—校正方法可能不稳定。根据经验，当小误差递减传播时，结果稳定；当小误差递增传播时候结果不稳定。当在区间上使用的步长太大时，会出现不稳定，有时表现为计算解的振荡性。采用较小的步长可使振幅减小。如果使用步长控制，则Milne-Simpon方法应该使用如下的误差估计：（7.12）这是一类不动点迭代过程。可以证明，其的步长应该满足以下条件（7.13）第八章一阶微分方程组与高阶微分方程的数值解法8.1一阶微分方程组的数值解法设有一阶微分方程组的初值问题若记，，，则初值问题（8.1）可写成如下向量形式 如果向量函数在区域：连续，且关于满足Lipschitz条件，即存在L>0，使得对，，都有那么问题（8.2）在上存在唯一解。问题（8.2）与（1.5）形式上完全相同，故对初值问题（1.5）所建立的各种数值解法可全部用于求解问题（8.2）。8.2高阶微分方程的数值解法高阶微分方程的初值问题可以通过变量代换化为一阶微分方程组初值问题。设有m阶常微分方程初值问题引入新变量问题(8.3)就化为一阶微分方程初值问题然后用6.1中的数值方法求解问题（8.4），就可以得到问题（8.3）的数值解。最后需要指出的是，在化学工程及自动控制等领域中，所涉及的常微分方程组初值问题常常是所谓的“刚性”问题。具体地说，对一阶线性微分方程组其中，,A为m阶方阵。若矩阵A的特征值满足关系 则称方程组（8.5）为刚性方程组或Stiff方程组，称数为刚性比。对刚性方程组，用前面所介绍的方法求解，都会遇到本质上的困难，这是由数值方法本身的稳定性限制所决定的。理论上的分析表明，求解刚性问题所选用的数值方法最好是对步长h不作任何限制。第九章常微分方程模型数值解法在数学建模中的应用9.1耐用消费新产品的销售规律模型9.1.1问题的提出新产品进入市场后，一般会经历一个销售量逐渐增加然后逐渐下降的过程。据此在时间一销售坐标系给出的曲线称为产品的生命曲线，其形状呈钟型。然而对于耐用消费品，情况有所不同，其生命曲线在开始有一个小的高峰，然后是一段平坦的曲线，甚至会下降，而后再次上升，达到高峰，从而呈双峰形曲线。如何解释这一似乎与传统的产品生命曲线理论相矛盾的现象昵?澳大利亚的斯蒂芬斯和莫赛观察到购买耐用消费品的人大致可以分为两类：一类是十分善于接受新事物的，称为“创新型”顾客，他们往往从产品的广告，制造商提供的产品说明书和商店的样品了解了产品的功能和性能后立即决定是否购买；另一类顾客则相对比较保守，称为“模仿型"顾客，他们要根据若干已购买该商品的用户的实际使用经验所提供的口头信息来决定是否购买。本节经过细致的分析，建立数学模型，对这一现象做出了科学的解释。9.1.2模型的构建将消费者获得的信息分为两类，一类称为“搜集型”的，来自广告、产品说明、样品，“创新型”的顾客在获得此类信息就可以做出是否购买的决定：另一类信息称为“体验型”的，即用户使用后获得的实际体验，经常以口头形式传播，“模仿型"顾客在获得此类信息后方能决定购买与否。设K为潜在的用户总数，K和置分别为其中的“创新型”和“模仿型”人数，又设为时刻已购买商品的顾客数，而和分别表示其中的“创新型”和“模仿型”顾客数，设为时刻中已经获得“搜集型”信息的人数，那么由于这部分信息可以直接从外部获得，也可以已经获得这种信息的人群中获得，于是有类似于巴斯模型的建立有 由于获得了“搜集型”信息的“创新型”顾客立即决定是否购买，于是应有对“模仿型”顾客，可以从已购买该商品的“创新型”或“模仿型”顾客中得到信息，因此有这里，忽略了顾客购买该商品后需要有一段短暂的试用才会传播体验信息的滞后作用。综上，斯蒂芬斯一莫赛模型是一常微分方程组的初值问题模型：两为时刻购买该商品的总人数。9.1.3模型的求解很容易求出斯蒂芬斯一莫赛模型中的解析解。其中，表示外部信息使“创新型”顾客购买新产品的比率；表示口传信息使“创新型”顾客购买新产品的比率；表示口传信息使“模仿型’’顾客购买新产品的比率。对于斯蒂芬斯一莫赛模型中的解析解则不能求出，于是可以用Adams四阶预测—校正公式求得，即使用求得，且在它的精度要求达到很高情形下求出。利用上述公式给出的数值算法，通过数学软件实现。具体程序如下： 设方程（9.3）中的于是有下面程序：,y(0)=0`)S=(K2*exp(t*a*b+t*a*K2+log(b/K2)/(b+K2)*b+log(b/K2)/(b+K2)*K2)-b)/(1+exp(t*a*b+t*a*K2+log(b/K2)/(b+K2)*b+log(b/K2)/(b+K2)*K2))9.2司机饮酒驾车防避模型的数值解法在2004年全国大学生数学建模竞赛题中有一个关于司机饮酒驾车模型。问题的提出《车辆驾驶人员血液、呼气酒精含量闽值与检验》国家新标准规定，车辆驾驶人员血液中的酒精含量大于或等于20毫克／百毫升，小于80毫克／百毫升为饮酒驾车，血液中的酒精含量大于或等于80毫克／百毫升为醉酒驾车。大李在中午12点喝了一瓶啤酒，下午6点检查时符合新的驾车标准，紧接着他在吃晚饭时又喝了一瓶啤酒，为了保险起见他呆到凌晨2点才驾车回家，又一次遭遇检查时却被定为饮酒驾车，这让他既懊恼又困惑，为什么喝同样多的酒，两次检查结果会不一样呢?请你参考下面给出的数据(或自己收集资料)建立饮酒后血液中酒精含量的数学模型，并讨论以下问题：1．对大李碰到的情况做出解释：2．在喝了3瓶啤酒或者半斤低度白酒后多长时间内驾车就会违反上述标准，在以下情况下回答：1)酒是在很短时间内喝的；2)酒是在较长一段时间(比如2小时)内喝的。3．怎样估计血液中的酒精含量在什么时间最高。4．根据你的模型论证：如果天天喝酒，是否还能开车?5．根据你做的模型并结合新的国家标准写一篇短文，给想喝一点酒的司机如何驾车提出忠告。参考数据1.人的体液占人的体重的65％至70％，其中血液只占体重的7％左右：而药物(包括酒精)在血液中的含量与在体液中的含量大体是一样的。2．体重约70kg的某人在短时间内喝下2瓶啤滔后，隔一定时间测量他的血液中滔精含量(毫克／百毫升)，得到数据如下:时间(小时)0.250.50.7511.522.533.544.55酒精含量306875828277686858515041时间(小时)678910111212141516酒精含量3835282518151212774 9.2.1模型假设l、驾驶司机没有其他疾病，消化系统良好，属于健康人群，其体重为70kg左右。2、酒精在人体内的循环系统分为胃腔系统(系统I)和体液系统(系统II)，两个系统的容积(即血液体积或酒精分布容积)在过程中保持不变。3、酒精从系统I向系统II的转移的速率系数，及向体外的排出的速率系数，与该系统的酒精浓度成正比，这两个速率系数、是由人体的身体机能所决定的常数。4、循环过程只考虑由体外进入胃腔，再由胃腔进入体液，最后由体液排除体外。不考虑人体其他机体对酒精的吸收，体液的变化可以忽略而保持一定。5、符号说明：：酒精进入胃腔的速率，设为常数：测试时间(小时)：饮酒时间(小时)：时刻人体胃腔中的酒精含量(毫克／百毫升)：胃腔中初始酒精量(毫克)：刚喝完酒时胃腔中的酒精量(毫克)：酒精由胃腔转移至体液的速率系数：酒精由体液排出体外的速率系数：酒精由胃腔转移至体液的转移速率(毫克／小时)：时刻人体体液中的酒精含量(毫克／百毫升)：时刻人体体液中酒精浓度(毫克／百毫升)：人体体液的体积(百毫升)：体液系统中初始酒精浓度(毫克／百毫升)记为：酒精排出体外的速率(毫克／小时)9.2.2模型建立 由酒精在人体内吸收、转移规律的特点，应用药物动力学原理建立人体内胃腔与体液循环系统模型，可用微分方程描述其动态过程。一般情况---长时间饮酒，原身体内有残余酒精1)胃腔系统过程：胃腔系统I进入体外，酒精含量，初始.当时，有求解到当时，有解得综合（9.5）、（9.6）得到于是2)体液系统过程：体液系统Ⅱ体外酒精含量，体液,酒精浓度，初始，则有 当时，有令，，则当时，令,则于是对于短时间饮酒，体内有残余酒精可以解释如下：,从而即9.2.3模型求解上述模型的表达式(9.5)，(9.6)，(9.8)均可归结为常微分方程初值问题，对于其解可用上面介绍的数值解法的方法给出。在这里给出了一个计算数值解的程序。在模型中考虑长时间饮酒的情况，用MATLAB 计算出，当大李饮酒的时间达到1．865个小时，检测时其酒精含量是20．0276毫克／百毫升，正好超标。大李短时间继续饮酒8小时后体内酒精含量(对于上式(6)的求解程序)k1=1.8100;k2=0.2100;D0=51200/2+0.0132;c0=18.3404;t=8;a=(k1*D0)/(v*(k1-k2));c=c0*exp(-k2*t)+a*(exp(-k2*t)-exp(-k1*t))运行结果：c=15.4695大李长时间饮酒后体内酒精含量（对于上式（9.9）的求解程序）k1=1.8100k2=0.2100v=447.867t0=1.865d0=0.0132f0=51200/4c0=18.3404a=f0/(v*k2)b=(f0-k1)/(v*(k1-k2))d1=f0/k1+(d0-f0/k1)*exp(-k1*t)c1=c0+a+b*exp(-k1*t)+(c0-a-b)*exp(-k2*t0)a1=d1/(v*(k1-k2))c=(c1+a1)*exp(-k2*(t-t0))-a1*exp(-k1*(t-t0))运行结果：c=20.02769.2.4模型评价1、本模型成功剖析了一部分想喝酒驾车的司机人员的心理。他们总想侥幸，然而事实不允许他们这么做，我们所做的工作让他们的这种心理无迹可遁，对促进交通安全也不无贡献。2、缺点：没有考虑其他可能的因素给饮酒驾车问题带来的影响，比如人的体重、司机的健康状况、交警检验程序不够科学等。求得的方案也许并不是最优的，但是相比之下比较满意的。9.2.5诚恳建议 广大的司机朋友们，为了您和他人的安全。请不要酒后驾车。但适量饮酒有助于健康。如果您是一位酒精爱好者，在一定的条件下，只要符合新的检验标准，饮酒也是无可厚非的，在这里根据我们所建立的饮酒驾车模型，得出血液的酒精随时间变化的关系。经分析，计算，检验基本符合实际情况。特向您诚肯地提供一些建议：当您辛苦了一天，晚上归来时，在保证至少6小时的休息时间的前提下，适当喝些酒，是不影响第二天工作的，但不要连续喝酒，更不要酒后驾车。有关数据附下，供您参考。以喝啤酒为依据，经过n小时后可以驾车，其与瓶数的关系如下表：饮酒量（瓶）123456789101112时间n(小时)61012131415161617171818备注：3瓶啤酒相当于半斤低度白酒9.2.6模型推广严禁酒后驾车?现有动力系统模型基本解决驾驶员饮酒量与停驾时间量化分析的交通难题，对驾驶员掌握驾驶时机有重要意义；模型的实际应用是当今社会非常急需，酒后驾车者被视为公路第一杀手；应用课题：如驾驶员饮酒量与停驾时间量化分析，驾驶员理论培训．肇事时血液中酒精浓度的反推算，车保赔偿等的研究。我们将研究初步结果送到相关单位专家手中，听取他们的意见。他们是本项目涉及到的实际应用领域的执行者和评判者。确切地说，他们的意见对我们进一步如何完善模型是非常有积极意义的。根据他们对该研究初步结果提出的宝贵意见：1、对于酒后驾驶的安全性，保险对酒后肇事的赔付等有着指导作用。2、对于法医学中所用的血中乙醇浓度反推生前饮酒量有意义。3、实验严谨，结论有明显的对比性．对于酒精在人体内的代谢浓度，有较完整数据。4、在“严禁酒后驾车”、“酒后驾车肇事不予赔偿”的规定和现实之间寻求一种合情合理又合法的新途径，提出了“安全饮酒”的新概念。5、“酒后安全驾车时刻表”，对于有效地预防和避免交通事故的发生有者一定的积极意义。6、研究提供了更科学、数字化地判断驾驶员是否应该驾车的依据，有利于解决驾驶员饮酒量与停驾时间量化分析的交通执法难题。主要参考文献[1]徐士良.数值分析与算法[M].北京:机械工业出版社,2007.[2]JohnH.Mathews,KurtisD.Fink.NumericalMethodsUsingMATLABFourthEedition[M]. PublishingHouseofElectronicsIndusty,2005.[3]刘萍,数值计算方法[M].北京:人民邮电出版社,2002.[4]胡健伟,汤怀民.微分方程数值方法[M].北京:科学出版社,2000.[5]陈志敏.龙格_库塔法及其Mathematica实现[J].武汉工程职业技术学报,2006.[6]张丽娟.常微分方程的Euler解法及其计算机实现[J].长春师范学院学报(自然科学版),2005.致谢论文得以完成，要感谢的人实在太多了。首先要感谢汪继文老师，本论文从选题到完成，每一步都是在汪老师的指导下完成的，倾注了汪老师大量的心血。汪老师指引我的论文的写作的方向和架构，并对本论文初稿进行逐字批阅，指正出其中误谬之处，使我有了思考的方向，他的循循善诱的教导和不拘一格的思路给予我无尽的启迪。汪老师要指导很多同学的论文，加上本来就有的教学任务，工作量之大可想而知，但在一次次的回稿中，精确到每一个字的的批改给了我深刻的印象，使我在论文之外明白了做学问所应有的态度。在此，谨向汪老师表示崇高的敬意和衷心的感谢！谢谢汪老师在我撰写论文的过程中给与我的极大地帮助。其次，论文的顺利完成，离不开其它各位老师、同学和朋友的关心和帮助。在整个的论文写作中，各位老师、同学和朋友积极的帮助我查资料和提供有利于论文写作的建议和意见，在他们的帮助下，论文得以不断的完善，最终帮助我完整的写完了整个论文。另外，要感谢在大学期间所有传授我知识的老师，是你们的悉心教导使我有了良好的专业课知识，这也是论文得以完成的基础。感谢所有给我帮助的老师和同学，谢谢你们！'