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

常微分方程求解的高阶方法毕业论文目 录第一章 前 言11.1案例引入微分方程概念11.2微分方程的基本概念11.2.1微分方程及微分方程的阶11.2.2微分方程的解、通解与特解11.2.3微分方程的初值条件及其提法21.2.4微分方程的解的几何意义.21.3从解析方法到数值方法概述31.4常温分方程的离散化4第二章 数值解法公共程序模块分析5第三章 欧拉（Euler）方法73.1 Euler方法思想73.2 Euler方法的误差估计83.3改进的Euler方法83.3.1梯形公式83.3.2改进Euler法9第四章 休恩方法104.1 休恩方法思想104.2休恩方法的步长和误差10第五章 泰勒级数法115.1泰勒定理115.2 N次泰勒方法12第六章 龙格-库塔（RungeKutta法）136.1龙格-库塔（RungeKutta）方法基本思想136.2 阶龙格-库塔（RungeKutta）方法公式14第七章 预报-校正方法157.1 Milne-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减少太多代码，况且这里的实际代码比较简单，所以为了减少系统的开销，采用Tubro C来实验。编程风格: 按照常微分方程数值解三个基本步骤：将问题离散化；建立递推格式；按步进法计算，所以求微分方程的数值解的算法框架都是相同的，不同的是所使用的递推形式不同，则可以用公共子程序来代替，对不同的方法的计算结果用统一的格式来显示，同时也可以比较不同方法的精确度4。公共程序模块如下: 这里为了良好地比较，选用可求解析解的一阶常微分方程作为讨论: (2.1)其解析式为 /* Filename: numerical_base.c */#include#include#include#define MAX 100int real = 1;double X MAX, Y MAX, ZMAX, FMAX, GMAX, CYMAX, EMAX;/*以下代码根据待求解的对象的特殊性进行赋值和在main()中选取*/double func (double x, double y) /计算各离散点处导数值 return y - 2.0*x/y ; double exact_value (double x) /计算各离散点解析解以测数值解精度 return sqrt( 2.0*x + 1.0) ; /* changing part end*/void cal_error() /计算误差值以对各种方法进行比较 int i; for(i=0; i=real; i+) CYi = exact_value( Xi ); Ei = fabs( CYi - Yi); return; void showtable_s() /微分方程组输出时用 /内容与 showtable() 类似 / 输出各离散点处的 X值， Y值， 导数值， 精确值， 误差值 / 分别对应于 Xk， Yk， Fk， CYk， Ekvoid showtable() /优化输出显示 int i,j; printf (n); for (i=0; i=78; i+) printf(=); printf(n); printf(%5s%8s%15s%15s%18s%14s,k,Xk,Fk,Yk,CYk,Ek); 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,Xi,Fi, Yi,CYi,Ei); printf(n); for(i=0; i0，使得对，都有 那么问题（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酒精含量38352825181512127749.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计算出，当大李饮酒的时间达到1865个小时，检测时其酒精含量是200276毫克百毫升，正好超标。大李短时间继续饮酒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、在“严禁酒后驾车”、“酒后驾车肇事不予赔偿”的规定和现实之间寻求一种合情合理又合法的新途径，提出