（基础数学专业论文）微分方程的数值模拟.pdf

河南大学硕士论文 摘要 对物理模型的数值模拟是判断模型理论分析是否正确的一个重要手段，其模拟 结果对实际应用具有一定的指导作用，因此，对研究的微分方程寻找稳定高效且精 度高的数值计算方法具有重要的理论意义和价值。 本文在第一章中简要介绍了m e m s 电子闸模型和生态植被模型的背景，然后对 有限差分方法进行了简单介绍。在第二章，针对m e m s 电子闸模型，用向前的e u l e r 差 分格式直接对方程进行离散，发现数值计算结果产生震荡。改用一阶精度的半隐半 显的筹分格式进行离散，发现数值计算结果相当好。在文章中，我们还直接对原方程 用两阶精度的中心差分格式进行离散，之后对方程中给出的初边值条件进行变形， 得到的数值计算结果与实际背景相符合，而且和半隐半显格式得出的数值结果相符 合。 在第三章，针对生态的植被模型，我们主要研究了在死亡率和降水率改变的情 况下，植被的变化情况。对这个无量纲化的常微分方程组，采用显式差分格式进行 离散，数值计算结果与实际结果相符。 本文，针对每种计算格式得出的数值计算结果都进行了详细的分析，每个数值 计算结果都是对理论结果是否正确的验证。 关键词：差分格式，数值模拟，m e m s ，植被模型 河南大学硕士论文 a b s t r a c t n u m e r i c a ls i m u l a t i o no ft h ep h y s i c a lm o d e li sa ni m p o r t a n ta n du s e f u lw a yt od e t e r - m i n et h e o r e t i c a la n a l y s i so ft h em o d e lt oc o r r e c to rn o t n u m e r i c a lr e s u l t sp l a yag r e a t g u i d i n gr o l e i np r a c t i c a la p p l i c a t i o n s ol o o k i n gf o rh i g h - p r e c i s i o n ，h i g he f f i c i e n c ya n d s t a b l ec a l c u l a t i o nm e t h o dh a si m p o r t a n tt h e o r e t i c a ls i g n i f i c a n c ea n dv a l u e w eb r i e f l y d e s c r i b e st h em e m sm o d e la n dt h ev e g e t a t i o nm o d e li nt h i sp a p e r ，t h e nt h eb a c k g r o u n d a b o u tt h ef i n i t ed i f f e r e n c em e t h o d so fk n o w l e d g ei si n t r o d u c e d s e c o n d l y ，w ed i r e c t l yu s ef o r w a r de u l e rd i f f e r e n c es c h e m ea b o u te q u a t i o n so fm e m s m o d e l ，b u tn u m e r i c a lr e s u l t ss h o c k s ow ec h a n g e dt h eo r i g i n a le q u a t i o ni n t oe q u i v a l e n t e q u a t i o n sg r o u pf o r m s ，u s i n gs e m i i m p l i c i td i f f e r e n c es c h e m e ，a n dt h en u m e r i c a lr e s u l t si s v e r yg o o d i nt h ep a p e r ，w ea l s od i r e c t l yu s ec e n t r a ld i f f e r e n c es c h e m ew i t ht h eo r i g i n a l e q u a t i o ni nt w o - o r d e ra c c u r a c y , a n dt h er e s u l ti sa l s oq u i t eg o o d ，e v e nt h en u m e r i c a l r e s u l t so fb o t hf o r m a t sa r et h es a l t l e i nv e g e t a t i o nm o d e lp a p e r ，w em a j r a ys t u d i e s v e g e t a t i o nc h a n g e ss t a t ew i t ht h ec h a n g i n gi nt h er a t eo fm o r t a l i t ya n dt h er a t eo fr a i n f a l l o nt h ed i m e n s i o n l e s so r d i n a r yd i f f e r e n t i a le q u a t i o n s ，p a p e r su s i n gt h ee x p f i c i td i f f e r e n c e s c h e m ef o rd i s c r e t e a tl a s t ，w ef i n dt h en u m e r i c a lr e s u l t sa n da c t u a ls t a t ea r et h es a m e t h i r d l y , w eh a v ec a r r i e do u tad e t a i l e da n a l y s i sf o rt h en u m e r i c a lr e s u l t so fe a c h d i f f e r e n c es c h e m e ，w h i l ee a c hn u m e r i c a lr e s u l t sa r et h ev e r i f i c a t i o no ft h e o r e t i c a lr e s u l t s t oc o r r e c to rn o t k e y w o r d s ：d i f f e r e n c es c h e m e ，n u m e r i c a ls i m u l a t i o n ，m e m s ，v e g e t a t i o nm o d e l i i 关于学位论文独立完成和内容创新的声明 本人向河南大学提出硕士学位中请。本人郑重声明：所呈交6 勺学住论文是 本人在导师酌指导下独立完成的。对所研究的课题有新的见解。据我所知，除 丈中特别加以说明、标注和致谢的地方外，论文中不包括其他人已经发袁或撰 写过的研究成果。也不包括其他人为获得任何教育、科研机构的学位或证书而 使用过的材料。与我一同工作的同，李黠奉柳流所做的任何贡南爰均已在论文中作 了明确的说明并表示了谢意了二：；：筑。j 。i ，“夸i ：! 学位聋蔷采海位淼? 墨? 。i 三：主亡垂佥 学位睁讯，i 烤住论疲雌渚堵；名：i 芝! l至! 垒! 誓曩砜觏j | ，奠奠7 疆鬈z 荔鞭；茹毒。瓣。7in ：。 x 菇澎毪 攀。 锄l 摒：。多薅，受日 磐：爹搭：未“二i 豢芝鼍 琴，基于学位论文著作极使讯授权书一零 垆f ”露i ：麓，聚波蕤磊 掌 本人经河韵大学审核批糖缝子硕士擘拉。作为学位论文酚作者，本人完全 7 解并同意河南诙擘有关保留绶罔学隹穗囊韵要求，即河满大学有权向国家 图书馆、科研信息机构、数据收蓑牟几掏和本校图书馆等提供擎位论文( 甄质文 本和电子文本) 以供公硷检索、查阅：；：漆认援权河赢失学出于宣扬、展览学校 学术发展和进行学术交流等7 酗的，。可瞒苯取焉卿：缩印、扫描和拷贝等复制手 段保存、汇编学位论文( 甄质文本和电子文本) 。 ( 涉及保密内容的学位论文在解密后适用本授权书) ；i王瓜 学位获得者( 学住论文作者) 釜名：2 1 l尘! 垒1 2 q 学位论文指导教师釜名： 2 d co 年f 月爹曰 曼邀。 jo - q - 月劈曰 第一章预备知识 1 1背景知识 文章的：i ：作之一是以m e m s 为背景的微电子机械系统动态电子闸模型的数值模 拟，m e m s 是指集微型传感器、执行器以及信号处理和控制电路、接口电路、通信和 电源于一体的微型机电系统。近些年来，m e m s 在实际中被广泛的应用1 ，2 】。概括 起来，m e m s 具有以下几个基本特点，微型化、智能化、多功能、高集成度和适于大 批量生产。m e m s 技术的目标是通过系统的微型化、集成化来探索具有新原理、新 功能的元件和系统。这是一种典型的多学科交叉的前沿性研究领域，几乎涉及到自 然及工程科学的所有领域，如电子技术、机械技术、物理学、化学、生物医学、材料 科学、能源科学等。由于尺寸缩小带来的影响( s c a l i n ge f f e c t s ) ，许多物理现象宏观 世界与微观世界有很大区别，因此许多原来的理论基础都会发生变化，宏观世界的 理论不能完全适用于微观世界。微观世界的理论是很多学科交叉的前沿学科，微观 世界理论的成熟，会给我们的生活带来很多便利，m e m s 技术与航空航天、信息通 信、生物化学、医疗、自动控制、电子消费以及兵器等应用领域相结合，制作出符合 各领域要求的微传感器、微执行器、微结构m e m s 器件等，微电子动力系统系统能 给我们的生活带来很多便利，然而微观世界的理论研究相当困难。因此，我们用数 值模拟的方法来研究该模型。本文所做的工作之一就是从数值计算的角度出发，来 模拟已有的理论模型 3 】，数值计算出影响动态电闸工作的某些元素的数据，对于各 种微型元器件的实际制造提供量化的理论指导。 近年来，随着全球变暖的日益加速，环境问题已经成为人们关注的焦点，而生态 环境沙漠化成为当前亟待解决的热点问题之一。沙漠化是中国北方广大地区主要的 土地退化形式，沙漠化是干旱、半干旱及部分半湿润地区由于人地关系不相协调所 造成的以风沙活动为主要标志的土地退化现象，沙漠的形成演变主要受控于气候变 化，而人类既是沙漠化的导致者，也是沙漠化的受害者，更是沙漠化的防治者，所以 自觉控制人的行动，采取合理的防治措旌，在一定程度上能解决植被沙漠化的问题。 关于植被态与水源和地形的关系，在参考文献 4 中已经体现，关于干旱，半干旱生 态的问题，以及生长于山坡的植被等问题已经被土壤植被专家广泛的研究f 5 ，6 ，7 】。 本文从数值模拟的角度去验证理论工作是否有实际意义，给如何治理沙漠化提供一 河南大学硕士论文 些理论的指导。关于植被态和沙漠态的模型的研究已经有很多的理论和数值分析成 果，在这方面做的工作国内外都已经相当成熟了。2 0 0 4 年以 8 】为背景而做的 9 只是 在模型的基础上加以分析，在这个模型中，作者的主要工作是在考虑降水率小于死 亡率，山坡地形和水分蒸发等因素的情况下，生态系统的植被态出现沙漠化情况而 做的数值分析。我们做的工作是以 3 为背景，数值模拟无量纲化的常微分方程组 4 】。 在文章中，不但考虑f 9 】中降水率大于死亡率的情况，而且还考虑到降水率等于死亡 率，降水率小于死亡率的两种变化情况，以及植被随降水量周期变化的情况。 1 2简介有限差分格式 应用有限差分方法数值求解微分方程，必须先把连续问题进行离散化 1 0 。为 此首先要对求解区域进行网格剖分，由于求解的方程各不相同，因此求解区域也不 尽相同，所以网格剖分的方法也不同。下面用具体的例子说明不同求解区域的网格 剖分情况，并引入一些常用的术语。 例1 双曲型方程和抛物型方程的初值问题，求解区域是 d = ( z ，t ) i o 。 0 ，有简记为丁，称其为时间步长。这样两族网格线可以写作 。= 巧= j a x = j h ，j = 0 ，士1 ，士2 ， t = t n = n a t = 仃7 ，n = 0 ，士1 ，士2 ， 网格节点( ，t n ) 有时简记为( 歹，礼) 。如图1 1 例2 双曲型方程和抛物型方程的初边值问题，设其求解区域是 d = ( z ，t ) l o z z ，t o ) 2 河南大学硕士论文 图1 1 网格剖分 网格点 这个区域的网格由平行于t 轴的直线族z = x j ，j = 0 ，1 ，j ，与平行于z 轴的直线 族t = t n , n = 0 ，1 ，所构成，其中x j = j a x = j h ，a x = h = j t ，t n = n a t = n 7 - ，网格 剖分的情况，如图1 2 其中内部的点叫着内点，边界的点叫边界点。上面我们列举 了不同求解区间的网格剖分情况，以后遇到具体问题我们再做进一步的讨论。 图1 2 网穑制分 数值求解方程可以构造很多差分格式，他们都是基于t a y l o r 展开得到的，它的 基本思想就是用差商来近似等于偏导数。 虹地幽t= 【警 了+ d ( 丁) ；。 亟血崆2 t 亟型= 铷3 + d ( 丁2 ) ；。况1 1 。” 迎盟划h= 【知+ 。( ) ；o a z o j 。、一” 3 河南大学硕士论文 亟出止2 h 捌= 塞 7 + d ( 坝。a zo j 。一” 啦出盥嚏盐捌=02unh j j + 。( h 2 ) 口z 。 上面这些式子去掉右边第二项即用差商来代替微商【1 0 】。 1 3有限差分格式的截断误差，相容性，收敛性及稳定性的概念 对于扩散方程的初值问题： 百o u ：口0 丽2 u ，z r , t 0 ， ( 1 1 ) 瓦2 口面，z u ， 【l 1 ) t 正( z ，0 ) = 夕( z ) ，z r ，a 0 通过网格剖分之后，用差商来代替微商之后，可以写成如下形式： _ u f _ - 哆一口盟挚扎 ( 1 2 ) j = 0 ，士1 ，扎= 0 ，1 ， 下面给出几个记号： a + t u ( x ，t ) = u ( x ，t + 7 ) 一u ( z ，t ) ， 磋u ( z ，t ) = 钆( z + h ，t ) 一2 u ( x ，t ) + 乱( z h ，t ) 事实上，通常定义截断误差t ( x ，t ) 为： t ( x , t ) = 睾倒( 叫) 一。嘉磋札( x , t ) 上面式子右端做t a y l o r 展开后，得到的主要部分就是所求方程的截断误差。由上面的 定义可以看到，扩散方程的显式格式( 1 2 ) ，用微分方程的解乱( 巧，t 竹) 来代替( 1 2 ) 式中的全部近似解“? ，整理之后得到结果就是截断误差。并且只要网格剖分的足够 细，及空间步长h 和时间步长丁足够小，差分格式的解就会无限接近微分方程的 解。由截断误差的定义可知，要求出一个差分格式的截断误差，只要把相应的微分 方程问题的充分光滑的解代入这个差分格式，再进行t a y l o r 展开就可以了。 相容性：从偏微分方程建立差分方程时，总是要求当丁一0 ，h _ 0 时差分方程 与微分方程充分接近，这就导致了差分方程的一个基本特征，差分格式的相容性。 也即当丁，九一0 时，若t ( x j ，t n ) 一0 ，这时称微分方程和差分方程是相容的。 河南大学硕士论文 收敛性：当空间步长丁和时间步长h 无限缩小时，差分格式的解是否逼近到微 分方程问题的解，这就是差分格式收敛性问题。 稳定性：利用有限差分格式进行计算时是按时间层逐层推进的，如果考虑两层 格式，那么要算几十1 层上的值需用到n 层的值，这个过程必然有舍入误差，从而 就要分析这种误差传播的情况，希望误差的影响不至于越来越大，以致掩盖差分格 式解的面貌，这便是所谓稳定性问题。 其实，在很多情况下，收敛性和稳定性是可以放在一起讨论的。本文中，我们 仅考虑对所研究模式的数值模拟，未涉及其相应的差分格式相容性，收敛性，稳定 性的进一步讨论。 5 第二章关于动态m e m s 电子闸模型的数值模拟 2 1简介m e m s 电子闸模型 本文主要以m e m s 为背景来研究和分析微电子动力系统的动态电闸模型 1 1 ，1 2 ， 文中主要采用差分方法进行数值模拟，用经典的显示差分格式，隐式差分格式，半 隐半显格式对模型进行分析和研究。并行板电子容器f 1 1 是微电子系统最常见的模 式，有动态和静态两种，在静态击穿条件下【1 3 】，通常采用稳定电压。在动态击穿条 件下，给出的稳定电压可以逐渐减少，这样，可移动模板震动的幅度就会变化。关 于动态的m e m s 电子闸随时间变化的模拟在 1 4 】中已经用有限元方法做过一部分工 作。我们研究的这个模型是被阶梯函数电压控制，具有弹性的动态反应微电子系统 电子闸【3 ，1 3 ，1 5 】，他的模型和并行板电子容器得的模型的数值计算结果是一致的。 我针对这个方程采用不同的差分格式逐步提高精度的方式来做数值分析。 下面的这个图就是研究的模型基础，主要研究移动木板的轨迹。方程( 2 1 ) 就是 有弹性 移动木板 固定木板 图2 1 模型 以参考文献【3 ，1 6 ，1 7 为背景而建立的模型： z ”一外丢南p o ( 2 1 ) z z + 互丽研。 o ( 2 1 ) x ( o ) = 0 ，x 7 ( 0 ) = 0 其中z ( t ) 是关于t 的函数。我们主要研究( 2 1 ) 关于初值i 、, - j 题的解，从给出的初值 可以看出方程的解都是非负的。方程( 2 1 ) 中z ( t ) 是移动木板运动的轨迹，v 表示 电压，表示绝热系数。 河南大学硕士论文 2 2m e m s 电子闸模型的数值离散及数值模拟 2 2 1 经典的显式差分格式及数值模拟 下面用差分方法数值求解方程( 2 1 ) ，首先，我们进行网格剖分，先对坐标轴进 行等距划分f 在这里平行z 轴的直线是等距的，平行于y 轴也是等距的) ，然后再做 两族平行于z 轴和t 轴的直线，对求解区域进行矩形分割：在这里t 是表示时间轴， z 轴上的数据代表木板运动的轨迹。令a t = 7 = t n + l t 竹为步长。 o 图2 2 网格剖分 为了计算简单，我们采用换元技巧，给方程降阶，把方程( 2 1 ) 转化为下面的方 程组形式： ：z 7 ，x ( o ) = 0 ，(22y ) 2 z ，j5 ，l， 可7 = 一z + 去南，可( o ) = o ( 2 3 ) 可2 一z + 互i 虿_ = f 1 习可( u ) 5 u l z j j 首先对方程组中的两个式子都用最常见的向前差分格式离散， x n + l z n 5 ， 1 y n + l 一- - y n - - x n t - 互1 南，珈一o o 。_ 。_ 一 一一 7 -zl 之十工一z 礼j 。 整理之后的式子， z n + 12r y n + x n ， 1t v 2 y n + l2 鲰一t x n + 互石万哥。 8 ( 2 4 ) ( 2 5 ) ( 2 6 ) ( 2 7 ) 河南大学硕士论文 很显然差分方程与微分方程是相容的，并且截断误差为d ( 丁) ，经过数值模拟实验之后 发现这种格式是不稳定的，产生数值震荡的比较厉害，数值结果( 参看文献 1 8 ，1 9 】做 的数据图) 如图1 ，图2 ，图3 ，图中纵轴表示z 即木板运动轨迹变化情况，横轴表示时 间t ，本章图形都如此，详情参看附录。 图3 从图1 ，图2 ，图3 中，我们可以看出这种差分格式是不稳定的，它产生的数值结果产 生震荡，与实际物理背景不符合。但是从这些数据中我们还可以观察出，当绝热系 数不变的情况下，电压越大，木板运动的振幅就越大；当电压不变的情况下，绝热 系数越大，振幅反而越小，这些情况是和理论分析结果相符合的，所以，本文就对 上面的差分格式进行改进。 2 2 2 半隐半显的差分格式及数值模拟 我们把( 2 2 ) ( 2 3 ) 改成如下差分格式， y n ：x n - - ：_ x n - 1 ， = 0 x o ， ( 2 8 )= ， 2 ， 【z 子j _ y n + l _ - - y n = - x n + 互1 志，珈朝 ( 2 9 ) 整理之后， x n = r y n + z 住一1 ，( 2 1 0 ) y n + l 一- - - y n - - t x n - i - 互1 禹 ( 2 1 1 ) 这个差分格式和微分方程( 2 2 ) ，( 2 3 ) 是相容的，截断误差为d ( r ) ，数值实验结果如 图4 ，图5 ，图6 ， 河南大学硕士论文 图6 图4 ，图5 ，图6 ，反映了模型在条件口2 。， ( 2 1 2 ) j j 一百石f ：_ r = _ 可， u ， 【2 上z ) 初值为z ( 0 ) = o ，z ( o ) = 0 ，很显然这个差分格式的截断误差是d ( 丁2 ) ，差分方程的 格式与微分方程也是相容的。 对( 2 1 2 ) 式整理之后得， z n + = 一+三高+2znxnt，(213-xnt2 xn(213)z n + 1 。 + 互泵_ f r 二j ：严+ 红n 一 一1 ，) 一( 2 - - t 2 ) x n - - x n - - 1 + 互1 群岛，( 2 1 4 ) 河南大学硕士论文 初值x ( o ) = o ，z 7 ( o ) = 0 ，直接用上面的初值不能做数值计算，我们就对初值进行改 进，对 z 。( 0 ) = 0 ， 用差商来代替微商可得， z 7 ( o ) 生必：0 ， 所以就得到， x ( - 1 ) = x ( o ) = 0 这样我们就得到了关于( 2 1 ) 式的中心差分方程， t = ( 2 - r 2 ) x n - - x n - 1 + 互1 群岛，( 2 1 5 ) x ( - 1 ) = x ( o ) = 0 下面就可以进行数值计算了。 我们先考虑在条件口2 ( + 1 ) 3 下，我们让电压变化，绝热系数不变，观察木板的 运动轨迹，图1 4 ，图1 5 ，图1 6 ，显示物体运动的轨迹是不规则的，在绝热系数不变， 电压逐渐变大的情况下，从图上看，木板的振幅是随着电压的增大而逐渐增大，实 际上，当绝热系数一定，电压增大到一定的程度，元器件会自动被击穿，或者说产 品不能正常工作了，也就是电压增到一定数值后，图中出现大振幅在客观实际中是 不存在的，上面的数据图在某种程度上反映了物体的非周期运动。也就是说，做这 种类型的零件，一定要注意材料的性质，绝热系数太大，太小都不好，上面的数据 也给人们控制电压提供一定数据基础，在做数值计算时，我们实际考虑的数据都是 大于零的，小于零是符合实际情况的，但我们不去考虑。其实，从图1 4 和图1 5 中，我 们可以观察到数值计算结果有小于零的情况，那是因为，本文在做数值计算时没有 对z 限制非负这个条件。如果在做数值计算时，对z 加上非负这个条件限制，那 么图1 4 ，图1 5 中就不会有小于零的那部分值了，这时，虽然会产生击穿现象，但从 图中。町以看出木板是在做周期运动。图1 6 显示的是，木板没有任何规律的运动状态， 它说明m e m s 模型不适应电压太大的情况。这也说明两阶精度的中心差分格式能更 精确的模拟出模型的运动状态。 接着，我们就考虑在满足u 2 = 三( f + 1 ) 3 的条件下，木板运动轨迹的各种情况。 1 3 河南大学硕士论文 图1 7 图1 8 图1 9 通过观察图1 7 ，可以发现在绝热系数为1 时，数值计算结果和半隐半显的数值计算 结果完全一样，与客观实际情况基本吻合。图1 8 ，在做数值计算时绝热系数取0 ，这 是一个在不可能状态下发生的情景，很显然木板运动的轨迹杂乱无章。图1 9 中，电 压取4 ，绝热系数取3 ，木板的运动轨迹也是杂乱无章的，所以，m e m s 模型不适 应于电压太大的情况。我们从这两个图中也发现了比较有意思的结果，它们除振幅 不一样之外，其他都是一样的。 1 4 第三章生态植被模型的数值模拟 3 1简介生态植被模型 文章主要考察生态的多样性，考察生态的多样性是一个很有价值的问题。本文 考察了一个在降水率足够大，死亡率相对比较小的情况下的一个植被态模型，这是 一个o d e 方程。我主要用有限差分的方法来分析这个模型，在分析的过程中我选择 比较经典的e u l e r 向前差分格式，在精度比较低的前提下也能得出比较好的结果。在 这一章中，我分析了降水率与死亡率的变化对植被态的影响程度，并且发现当降水 率周期变化时，植被态也是周期变化的，当生态环境未沙漠化时，时间足够大后，植 被稳定在平衡态附近，而且，当死亡率和降水率相等时，在这个模型里会出现混合 的生态问题。 下面是以 8 为背景的模型，无量纲化之后的模型【9 】，可以从沙漠态和植被态两个 方面去研究，从沙漠态方面的研究主要模型是， 雨o b ：w b p b ， ( 3 1 ) 尝：d a w + v v w + r w 一, k w b ( 3 2 ) 【，l 初值条件：b ( x ，0 ) = b o ( x ) ，w ( z ，0 ) = w o ( z ) 上面的这个模型，对于讨论植被的沙 漠态化问题非常好，对于讨论生态系统的植被态情况非常困难，所以我们就不考虑 上面的p d e 模型，我们考虑【1 0 】的o d e 模型， b t = b ( w p ) ， ( 3 3 ) w t = r w ( 1 + 入b ) ，t 0 ( 3 4 ) 初值条件为： b ( 0 ) = b 0 ，w ( o ) = w o 数学符号b 表示植被密度，符号w 表示水资源，初始条件已经给出，p 代表植被在 缺水条件下的死亡率，圭表示植被利用水的效率，r 表示降水量水平。t 表示时间 坐标，z 表示空间坐标。 1 5 河南大学硕士论文 关于方程! 1 1 ( 3 3 ) ，( 3 4 ) ，我们用差分方法做数值模拟。首先，进行网格剖分( 如前 面所示均匀剖分网格) ，分别对方程组中两个方程进行数值离散，在这里t 是表示时 间轴，b 是植被密度，w 是水资源，令a t = 1 = t n + l t 竹为步长。 在做数值模拟之前我们先了解一下方程的平衡态问题， b ( w p ) = 0 ，( 3 5 ) r w ( 1 + a b ) = 0 ( 3 6 ) 整理之后得到平衡态为， 溉= p ，( 3 7 ) b = 妻( 罢一1 ) ( 3 8 ) 由( 3 7 ) ，( 3 8 ) 式，可以看出在r p 的情况下，当时间趋于无穷时，植被密度和水 资源都基本趋于稳定的状态，这个稳定的状态和平衡态的值之间是有着密切的关系 的。 定理3 1 1 对于上面给定的任意初始值( b o ，w o ) ，在满足b o 0 ，w o o ，j o 。 ，在降雨水平足够高的前提下式子p 砂和p 纠的解，接近平衡态( b 1 ，阢) 2 q 。总之， 我们观察到唯一的植被态是渐近平衡态。 3 2e u l e r 甸前差分格式及数值模拟 把方程组( 3 3 ) ，( 3 4 ) 转变为差分方程组， b j + _ i _ - 一b j ：s j ( w j p ) ， w j _ + 1 _ - w j ：r 一嵋( 1 + 入岛) ，t o ， 初值为， b ( o ) = b o ，w ( o ) = w o 整理之后得， 马+ 1 = 7 马( 一p ) + w j ， w j + 1 = t r 一7 - ( 1 + a 岛) + w j ，t 0 ， 1 6 ( 3 9 ) ( 3 1 0 ) ( 3 1 1 ) ( 3 1 2 ) 河南大学硕士论文 初僵为， b ( 0 ) = b o ，w ( o ) = 关于( 3 9 ) ( 3 1 0 ) 式的截断误差，先把岛+ 1 ，w j + 1 在马，点处分别做t a y l o r 展开， s j + l 岛+ 百o w j 7 - + 。( 下2 ) ，( 3 1 3 ) w j + 1 = + 百o w j 7 + 。( 丁2 ) ( 3 1 4 ) 把( 3 1 3 ) ( 3 1 4 ) 式代入( 3 1 1 ) ( 3 1 2 ) 式，得到的截断误差为d ( 丁) ，并且差分方程和微 分方程也是相容的。 在文章中我们可以想象到，如果死亡率很高，植被会慢慢的消亡，取而代之的 就是土壤的沙漠化。现在，我们就考虑在死亡率保持不变，降水率足够大的情况下， 植物的变化规律。现在，我们就用最经典的e u l e r 差分格式，来观察一下这个模型是否 符合自然界客观规律，本章图中纵轴分别表示植被密度b 和水量密度w ，横轴表 示时间t 。 图2 0 图2 1 图2 2 图2 3 图2 4 河南大学硕士论文 图2 5 数值计算结果表明，在计算的过程中数值解是稳定的。上面的数值计算结果， 是在降水量大于死亡率的情况下得出的，即r 弘，从图2 0 ，2 1 ，2 2 ，2 3 ，2 4 ，2 5 可 以发现，文章中死亡率一定的情况下( 死亡率取0 1 ) ，水量越大，植物的植被态密度 就会越大，并且随着水量的增多植被的变化也是有规律的变化着，在水量r 0 2 5 时， 植被密度会随着水量的增加而增加，并且增加到一定程度后也会趋于稳定，而且这 个稳定的值和平衡态的值也有关系。数值计算结果也显示出，当植被密度很快增加 的时候水量密度会很快的减少，当植被趋于稳定态的时候，水量也处于一个稳定的 值。 上面所做的分析都是在r p 条件下得到的，关于冗 p 的情况，其实，也就 是在降水率较小时，植被密度的变化情况，这是关于文中提到的p d e 模型研究的问 题，这里不再讨论。然而关于冗= 弘情况，在o d e 的模型里，植被密度，和水量密 度变化的情况，我们先看它的数值计算结果图2 6 ，图2 7 ，图2 8 。 图2 6 图2 7 图2 8 从数值模拟的结果中发现，在r = p 植被密度随着时间的发展也是有规律的， 都是从植被态慢慢演化为沙漠化形态，并且死亡率越高，植被演变为沙漠化的形态 就越快，随着植被沙漠化的演变，土壤的水资源和各种其他情况基本趋于稳定，但 是上面的模型还不能准确的表达出沙漠化之后的规律。 】8 河南大学硕士论文 3 3 数值模拟植被周期变化规律 定理3 3 1 如果降水率r 是一个关于t 的t 周期函数，常微分方程p 剀p 圳的 解也是一个t 周期的，在降雨率满足 ( 1 + 等) 一1z ? ( 1 - e - ( 1 + 譬，c r 一出p t 的条件下，当植被在周期植被状态下植被b ( ) 的周期解永远不会消失。定理的理 论证明参阅文献胆u 。 在定理中，r o ：【r ( t ) l o t t ) ，仃= r a i n 1 ，肛) ，在整个文章中，我们都是在 o 0 ，下面我们就来看它的实际数值计算结果，图2 9 是 在p ：o 1 ，r ：3 0 + s i n ( o 1 木歹) 的条件下的得到的，图3 1 ，3 2 是图2 9 中的周期变化的 一部分内容。图3 0 是在p ：0 1 ，r = 3 0 + 3 0 木s i n ( o 1 牢歹) 条件下得到的，图3 3 ，3 4 同 样为图3 0 的一部分。 图2 9 图3 0 图3 1 图3 2 圈3 3 1 9 河南大学硕士论文 图3 4 分析图2 9 ，3 0 ，我们可以发现在满足上述定理的条件下牛态都是植被形态，而且在降 水率周期变化下，植被也是周期变化的，并且植被周期变化的长短，是和周期函数 前面的系数有关系的，植被密度的大小，也和降水率，死亡率都有密切的关系，他 们的规律在数值计算结果中可以看出来。 第四章总结与展望 本文主要介绍了有关有限差分格式的一些背景知识，简单介绍了如何运用有限 差分格式对微分方程进行数值离散，进一步阐述了，怎么样运用有关数值模拟的一 些知识，对本文中已知模型进行研究。 本文的第一个工作，研究的是以m e m s 为背景而建立的模型，针对这个微分方 程，主要采用了经典的向前e u l e r 差分格式，半隐半显差分格式和中心差分格式来进 行数值模拟，从数值计算的角度来说，后两种格式得出的结果几乎相同，这说明我 们所做的工作是成功的。事实上，本文除了采用半隐半显的差分格式外，还可以采 用半显半隐的差分格式。关于本文的这个问题，还可以用r u n g e k u t t a 的方法进行讨 论，关于文中出现的绝热系数和木板运动轨迹的关系还可以再深入的进行讨论。 本文的第二个工作，是以生态系统的植被态为背景而建立的模型，关于这个微 分方程组，本文主要介绍了，运用经典的向前差分格式来进行数值模拟的结果。本 文的讨论都是在降水率大于死亡率，降水率等于死亡率的前提下展开的。在降水率 小于死亡率的情况下，植被就出现沙漠化现象了，这就是另外一个模型了，关于这个 问题本文不再讨论。其实，针对文中的方程组，我们还可以采用经典的向后e m e r 差 分格式，或者一前一后，一后一前的差分格式，还有r u i l 盼k u t t a 方法进行分析讨论， 但是采用这几种方法进行讨论比较困难，有兴趣的同学可以尝试一下，本文都不再 做介绍了。最后，我们所做的结果在附录里面都能体现。 展望：在有限差分方法中，c n 格式是无条件稳定的。文中虽然用不同的差分 格式对方程进行数值模拟，但没有尝试这种格式，我们相信c n 格式做出的结果会 有更大说服力。 附录 2 1 8o o 6o o 4o o 2o o o 2o o 一4o o - b0o 6o 4o 2 o 0 2o 4o - 60 河南大学硕士论文 020 04 0060 08o01000 图1 木板移动的轨迹( 电压1 ，绝热系数1 ) 02 004006008001000 图2 木板移动的轨迹( 电压1 ，绝热系数2 ) 河南大学硕士论文 6 ooo 4 ooo 2 ooo o 2oo0 4ooo 6ooo - 80oo o 3o o 25 o 2o o 15 o ，1o 0 o5 o 0o 02 004006 0080o1000 图3 木板移动的轨迹( 电压2 ，绝热系数1 ) 02 00 400 600800 10 00 图4 木板移动的轨迹( 电压1 ，绝热系数1 ) 河南大学硕士论文 o 12 o 1o o o8 o o6 o o4 o 02 o oo o o7 o o 6 o 0 5 o o 4 o 03 0 02 o o1 o oo o 0 1 02004 0 060o80 01000 图5 木板移动的轨迹( 电压1 ，绝热系数2 ) 0 200 4006008001000 图6 木板移动的轨迹( 电压1 ，绝热系数3 ) 2 4 河南大学硕士论文 1o0o 5o0 0 - 50o - 10o0 60 4 0 2 o 0 - 2o 40 - 6 0 020 040 060 08001000 图7 木板移动的轨迹( 电压3 ，绝热系数1 ) 0 20040o 6008001000 图8 木板移动的轨迹( 电压4 ，绝热系数1 ) 河南大学硕士论文 15o 1oo 5o o - 5 o 1oo - 15o 1 o o 8 o 6 0 4 o 2 o 0 020040 06 00800100o 图9 木板移动的轨迹( 电压5 ，绝热系数1 ) 0200 400 60080 0100o 图l o 木板移动的轨迹( 电压是扼( 数据实验采用的1 4 1 4 ) ，绝热系数1 ) 河南大学硕士论文 0 3 o o 2 5 o 2 0 0 15 o 1o o o 5 o 0 o o 12 0 o o 08 o o6 0 o 4 o 02 0 oo 0200 400 60080 01000 图1 1 木板移动的轨迹( 电压1 ，绝热系数1 ) 0200 4 00 6008 00 1000 图1 2 木板移动的轨迹( 电压1 ，绝热系数2 ) 河南大学硕士论文 o o7 o o 6 o o 5 o 0 4 o o3 o 02 0 o1 0 0 0 o 0 1 300 2o0 1oo o 1oo 2 0o 3oo 020040060 080 01000 图1 3 木板移动的轨迹( 电压1 ，绝热系数3 ) 0200400 60 08 00 10 00 图1 4 木板移动的轨迹( 电压2 ，绝热系数1 ) 河南大学硕士论文 4o oo 20 o0 o 20 oo 40 0o 1oo 50 o 50 - 1oo 0200 40 0 6008001000 图1 5 木板移动的轨迹( 电压3 ，绝热系数1 ) 020 04 0060080 010 00 图1 6 木板移动的轨迹( 电压4 ，绝热系数1 ) 河南大学硕士论文 1 0 o 8 o 6 0 4 0 2 o o 02004006008001000 图1 7 木板移动的轨迹( 电压以( 数据实验采用的是1 4 1 4 ) ，绝热系数1 ) 40 3o 20 1o 0 10 20 30 4o 02 004 006o0800100 o 图1 8 木板移动的轨迹( 电压0 5 ，绝热系数0 ) 河南大学硕士论文 15o 1oo 5o o - 5 0 1oo 150 1 5 1 4 1 3 1 2 1 1 10 0 ，9 0 8 o 7 0 6 o 5 o 4 0 3 0 2 0 1 o 0 o20040 060080 01000 图1 9 木板移动的轨迹( 电压4 ，绝热系数3 ) o 图2 0 植被密度b 和水量密度w 的时间发展图像( 死亡率0 1 ，降水率0 1 8 ) 3 1 1 5 1 4 1 3 1 2 1 1 1 o o 9 o 8 0 7 o 6 0 5 0 4 o 3 o 2 o 1 0 o 河南大学硕士论文 图2 l 植被密度b 和水量密度w 的时间发展图像( 死亡率o 1 ，降水率0 2 ) 1 5 1 4 1 3 1 2 1 1 1 0 0 9 o 8 0 _ 7 0 6 o 5 0 4 o 3 0 ，2 0 ，1 0 0 o 图2 2 植被密度b 和水量密度w 的时间发展图像( 死亡率0 1 3 2 1 0 0 0 01 2 0 0 0 降水率0 2 5 ) 河南大学硕士论文 2 o 1 8 1 ，6 1 4 1 2 1 o 0 8 o 6 0 4 0 2 o o 2 图2 3 植被密度b 和水量密度w 的时间发展图像( 死亡率o 1 ，降水率0 3 ) 3 0 2 5 2 o 1