（计算数学专业论文）分数阶微分方程的adomian方法.pdf

摘要 随着分数阶微分方程越来越多的被用来描述光学和热学系统，流变学及材料和 力学系统，信号处理和系统辨识，控制和机器人及其他应用领域中的问题，因此， 对分数阶微分方程计算方法的研究显得尤为迫切本文旨在分数阶微分方程的数 值计算方向作出努力 该论文共分四章，第一章简要地回顾分数阶微分方程的历史和几种定义 第二章对分数阶微分方程的预估校正法进行简要讨论，在此基础上，运用数值 跟踪法得到分数阶b r u s s e l a t o r 系统有极限环的有效维数 第三章介绍分数阶偏微分方程的a d o m i a n 方法，给出了求解非线性的时间分 数阶偏微分方程近似解的方法，并运用求解分数阶n a v i e r - s t o k e s 方程 第四章在分数阶微分方程预估校正法和a d o m i a n 方法的基础上，从新的观点 来设计分数阶微分方程的计算格式，即基于a d o m i a n 技巧的数值计算方法此方 法有效地解决非线性分数阶常微分方程和分数阶偏微分方程的数值计算结合基于 a d o m i a n 技巧的数值计算方法和数值跟踪法得到分数阶l o r e n z 系统有混沌吸引子 的有效维数；运用基于a d o m i a n 技巧的数值计算方法对时间分数阶b u r g e r s 方程进 行数值求解此章内容是本文的核心内容 关键词分数阶微分方程，预估校正，数值跟踪，a d o m i a n 方法 ab s t r a c t f a c t i o n a ld i f f e r e n t i a le q u a t i o n sa r ei n c r e a s i n g l yu s e dt om o d e lp r o b l e m si na c o u s t i c s a n dt h e r m a ls y s t e m s ，r h e o l o g ym a t e r i a l sa n dm e c h a n i c a ls y s t e m s ，s i g n a lp r o c e s s i n ga n d s y s t e m si d e n t i f i c a t i o n ，c o n t r o la n dr o b o t i c s ，a n do t h e ra r e a so fa p p l i c a t i o n a st h ei n t e r - d i s c i p l i n a r ya p p l i c a t i o n sa r ec h a r a c t e r i z e de l e g a n t l yb yf r a c t i o n a ld i f f e r e n t i a le q u a t i o n s ，i t g e t sm o r ea n dm o r ei m p o r t a n tt os e e kn u m e r i c a lc o m p u t a t i o n sf o rf r a c t i o n a ld i f f e r e n t i a l s y s t e m s t h i sp a p e rc o n c e n t r a t e so nn u m e r i c a lc o m p u t a t i o nf o rf r a c t i o n a ld i f f e r e n t i a l e q u a t i o n s t h i sp a p e rc o n s i s t so ff o u rc h a p t e r s m o r ed e t a i l e d ，t h ef i r s tc h a p t e rb r i e f l yr e v i e w s t h es i m p l eh i s t o r ya n dt h ed e f i n i t i o n so ff r a c t i o n a ld e r i v a t i v e sa n di n t e g r a l s i nc h a p t e r2 ，w ef i r s tr e c a l lt h ep r e d i c t o r c o r r e c ta p p r o a c hf o rn u m e r i c a l l ys o l v i n g f r a c t i o n a ld i f f e r e n t i a le q u a t i o n ，t h e nw ea p p l yt h en u m e r i c a lt r a c k i n gt e c h n i q u et of i n d i n g t h e “e 伍c i e n td i m e n s i o n s o ff r a c t i o n a lb r u s s e l a t o r i nt h ef o l l o w i n gc h a p t e rt h ea d o m i a n sm e t h o df o rp a r t i a lf r a c t i o n a ld i f f e r e n t i a l e q u a t i o ni si n t r o d u c e d a n dt h ep r o c e d u r ef o rs e e k i n ga n a l y t i c a la s y m p t o t i c a ls o l u t i o n s o ff r a c t i o n a lp a r t i a ld i f f e r e n t i a le q u a t i o n sa r ea l s ol i s t e d f i n a l l y , t h et i m e - f r a c t i o n a l n a v i e r s t o k e se q u a t i o ni st a k e na sa ne x a m p l e i nt h el a s tc h a p t e r ，t h ei d e ao fp r e d i c t o r - c o r r e c ta n dt h ea d o m i a n sm e t h o df o r t h en u m e r i c a ls o l u t i o nt of r a c t i o n a ld i f f e r e n t i a le q u a t i o na r ec o n s t r u c t e d a p p l i e dt h e n o v e la l g o r i t h mt oo r d i n a r yn o n l i n e a rf r a c t i o n a ld i f f e r e n t i a le q u a t i o n se s p e c i a l l yp a r t i a l n o n l i n e a rf r a c t i o n a ld i f f e r e n t i a le q u a t i o n ，t h er e s u l t sr e v e a lt h a tt h ep r o p o s e dm e t h o di s v e r ye f f e c t i v e t h en o v e la l g o r i t h mi sc o m b i n e dw i t ht h en u m e r i c a lt r a c k i n gt e c h n i q u et o f i n dt h e “e f f i c i e n td i m e n s i o n s o ff r a c t i o n a ll o r e n zs y s t e m t h en u m e r i c a ls c h e m e sf o rt h e t i m e - f r a c t i o n a lb u r g e r se q u a t i o na r ed e s i g n e da n dv e r i f i e d t 豳c h a p t e ri st h ec e n t r a l p a r to ft h i sp a p e r k e y w o r d sf r a c t i o n a ld i f f e r e n t i a le q u a t i o n ；p r e d i c t o r - c o r r e c t o r ；n u m e r i c a lt r a c k i n g t e c h n i q u e ；a d o m i a n sm e t h o d 1 1 原创性声明 本人声明：所呈交的论文是本人在导师指导下进行的研究工作除了文中特别 加以标注和致谢的地方外，论文中不包含其他人已发表和撰写过的研究成果参与 同一工作的其他同志对本研究所做的任何贡献均已在论文中作了明确的说明并表 示了谢意 占导b 身 本论文使用授权说明 本人完全了解上海大学有关保留、使用学位论文的规定，即：学校有权保留论 文及送交论文复印件，允许论文被查阅和借阅，学校可以公布论文的全部或部分内 容 ( 保密的论文在解密后应遵守此规定) 签名：址导师签名：乏茎笙日期：占扫 第一章前言 本章在引言部分简要阐述分数阶微积分的历史及研究现状，在接下来的部分 回顾分数阶微积分的定义 1 1 引言 早在1 6 9 5 年，l e i b n i z 给l h o s p i t a l 写了一封信，询问：整数阶导数的概念能 否自然推广到非整数阶导数l h o s p i t a l 对这个问题感到很新奇，作为回信也反问 了一个简单的问题；“如果求导的阶数为二分之一，那么将会是怎样的情况呢? 在这一年的9 月3 0 号，l e i b n i z 给l h o s p i t a l 回信写到：这将会导致悖论，不过 总有一天会得到有用的结果”这个特殊的f i 子1 6 9 5 年9 月3 0 号因此被认为是 分数阶微积分的确切诞生日早期对分数阶微积分有贡献的数学家包括l i o u v i l l e 、 r i e m a n n 等 在近三个世纪里，对分数阶微积分理论的研究主要在数学的纯理论领域里进 行，似乎它只对数学家们有用自从1 9 7 4 年第一届国际分数阶微积分会议之后， 它在各个领域都找到了自己的应用3 0 年来，它在理论和应用方面都取得了很大 的发展特别是近十几年来，分数阶的微分方程在科学的不同领域已得到广泛的 应用，许多学者指出分数阶微积分非常适合刻画具有记忆和遗传性质的材料和过 程，在经典模型中这些性质常常是被忽略的现在它已经成功应用到控制，金融， 物理，化学，水文地质等领域国外还成立了专门的分数阶微分方程建模组织 我们可以把分数阶积分”理解为“礼重积分”的自然推广，同样也可以把分 数阶微分”理解为饥阶导数 的自然推广不过，分数阶积分和微分的定义都有几 种，现在都发现了各自的应用领域但是，简单说来，分数阶微积分是关于任意阶 微分和积分的理论近些年来，一些关于分数阶微积分的奠基工作相续出版，其中 包括分别由下列作者写的书：k b o h l h a m 和j s p a n i e r 1 ，s g s a m k o ，a a k i l b a s 和o i m a r i c h e v 2 ，k s m i l l e r 和b r o s s 3 】，i p o d l u b n y 4 尽管分数阶微积分与经 典的微积分有相同的历史，且已经有了相当的发展，但是远没有经典的微积分理论 那样完善随着分数阶微分方程在越来越多的学科领域里出现，无论对分数阶微分 方程的理论还是计算方法的研究都显得尤为迫切然而由于分数阶导数是拟微分 l 2 0 0 8 上海大学硕士学位论文 2 算子，它的非局部性在对现实问题进行了优美刻画的同时，也给理论分析和数值计 算带来困难本文旨在数值计算方面上作出努力，并应用于分数阶偏微分方程的计 算 1 2 分数阶微分算子和积分算子的定义 1 2 1g a m m a 函数和b e t a 函数 g a m m a 函数是阶乘概念的推广，定义为 r ( z ) ：= e 一铲d t ，( 驼( z ) o ) ，u 同时，g a m m a 函数可用极限形式表示为 它具有如下关系式： r ( z ) ：= n l 。i m 豸耳1 _ ( 跪( z ) 。) r ( z ) r ( 1 一z ) 2 五丽1 1 ，o o ) b e t a 函数也可用g a m m a 函数表示 即川叫邺) = 躺 1 2 2 分数阶积分与导数 由微积分的知识我们知道，对于一个函数求n ( n n ) 重积分可简化 a 矿u = 芒可瓜叫铲1 缸( 丁) 打 2 0 0 8 上海大学硕士学位论文 3 将上式推广到非整数阶的情形，并使用g a m m a 函数可给出如下r i e m a n n - l i o u v i l l e 分数阶积分定义 定义1 2 1r i e m a n n - l i o u v i l l e 分数阶积分令u ( t ) 定义在区间( a ，b ) 上，p 0 ，则 阶数为肛的r i e m a n n l i o u v i l l e 分数阶积分定义为 r l 脚) ：2 高上( 一丁) p _ 1 饥( 丁) d r ( 1 2 1 ) 将经典的整数阶导数与分数阶积分算子作复合运算便可以给出如下r i e m a n n - l i o u v i l l e 分数阶导数的定义 定义1 2 2r i e m a n n - l i o u v i l l e 分数阶导数令u ( t ) 定义在区间( a ，b ) 上，p 0 ，n 是大于p 的最小整数，及盯= 佗一p ，则阶数为p 的r i e m a n n - l i o u v i l l e 分数阶导数义 为 舰d ：, t u ：= d 剐t ) = 高杀( z 。( 亡一矿_ 1 缸( r ) 蚍 ( 1 2 2 ) c a p u t o 分数阶导数和r i e m a n n - l i o u v i l l e 分数阶导数用几乎同样的方式推广了 经典导数，只是它们的微分和积分顺序不同 定义1 2 3c a p u t o 分数阶导数令u ( t ) 定义在区间( a ，b ) 上，p 0 ，扎是大于p 的最小整数( n 一1 p n ) ，及盯= n p ，则阶数为p 的c a p u t o 分数阶导数义为 c 咣t u ( t ) 2 cd 三? d n u ( 亡) 2f 高上( t - - 1 ) a - l u ( n ) ( 丁) d t ( 1 2 3 ) r i e m a n n l i o u v i u e 、c a p u t o 分数阶导数都对经典导数进行了比较完美的推广， 从数学意义角度来看r i e m a n n - l i o u v i l l e 分数阶导数与经典导数吻合，c a p u t o 分数 阶导数在区域的左端点和经典导数有可能不吻合尽管c a p u t o 分数阶导数在这方 面不够完美，但是它在工程中却被广泛地应用这是因为对于经典或者分数阶微分 方程来说，一个显著的特征就是必须给定初值来保证解的唯一性如果我们研究用 c a p u t o 分数阶导数定义的分数阶微分方程时，给定的初值与经典微分方程给定的 初值保持一致，在应用中有很清楚的物理意义而如果用r i e m a n n - l i o u v i l l e 分数阶 导数定义的分数阶微分方程，尽管它在数学上是严格的、优美的，但是在分数阶导 数没有很好的物理或者几何解释的情形下，给出它的初值是困难的本文以下研究 的分数阶微分方程均运用c a p u t o 分数阶导数定义 本论文共分四章，第一章为前言；第二章第一节简单讨论了分数阶微分方程的 预估校正法，第二节在前一节的基础上，数值跟踪分数阶b r u s s e l a t o r 系统有极限环 2 0 0 8 上海大学硕士学位论文 4 的最小有效维数；第三章将介绍分数阶偏微分方程的a d o m i a n 方法，并求非线性时 间分数阶偏微分方程的解析近似解；第四章在前两章的基础上提出了基于a d o m i a n 技巧的分数阶微分方程的数值计算新方法，并运用到分数阶微分方程特别是分数阶 偏微分方程的计算，这部分内容是本文的核心内容 l i l - - 章 数值跟踪分数阶b r u s s e l a t o r 系统有极限环的最小 有效维数 本章的第一节简单讨论分数阶微分方程的预估校正法；第二节是数值算例，数 值跟踪分数阶b r u s s e l a t o r 系统有极限环的最小有效维数 2 1 预估校正格式的描述 本节简要地回顾和改进数值求解如下分数阶微分方程的预估校正格式 5 - 1 0 ， c 磁z ( 亡) = y ( t ，z ( ) ) ，x ( k ( o ) = 鞯，k = 0 州1 一，陋卜1 ， ( 2 1 1 ) 其中r o t l 是大于o t 的最小整数该方程等价于第二类非线性v o l t e r r a 积分方程， 邢) = 薹1 丽t k + 击肛丁广加肌 ( 2 1 2 ) 对于问题( 2 1 2 ) 的计算，应用复合梯形公式来计算积分，其中以( z n 十l 一) 口一1 为 权函数取节点岛0 = 0 ，1 ，2 ，几+ 1 ) ，其中h = 叫n ，t n = n h ，n = 0 ，1 ，n z + 则对( 2 1 1 ) 进行离散，我们得到a d a m s - m o u l t o n 公式 其中 ，n + l2 它等价于 州) = 薹1 两t k + + r - - 砩博_ _ o 啪胞一钏， ( 2 1 3 ) 1 1 , 。+ 1 一( n q ) ( 礼+ 1 ) ， ( 佗一j + 2 ) 。+ 1 + ( n j ) a + 1 1 ， z c 如川= 。薹1 丽t k + 币，c 沁州饥+ 高妻+ 。，c 勺一岛执c 2 工4 ， 5 m lv i + 0 n | | 0 ，k 2 = 0 ，1 ，q 2 1 1 ， ( 2 2 1 ) ； c d o 一, t x 。( ) = l ( t ，z ( ) ) ，z s ( o ) = z 磐，q 。 0 ，k 。= 0 ，1 ，0 1 月一1 其中x ( t ) = ( z 1 ( ) ，z 2 ( 亡) ，x s ( t ) ) 定义“有效维数”为：iqi = o i l + q 2 + + o t 。， c 瑶t z ( t ) 如果满足( q 1 ，a 2 q 8 ) = ( 1 ，1 1 ) ，则系统( 2 2 1 ) 就是一般的整数阶微分 方程组 2 0 0 8 上海大学硕士学位论文7 文献 1 6 】中通过讨论分数阶l o r e n z 系统存在混沌吸引子的有效维数，提出是不 是对于一般的分数阶混沌系统存在一个通用的有效维数而使得它们仍然有混沌吸 引子文献 17 】运用频域的方法寻求分数阶混沌系统有混沌吸引子的有效维数，而 且有效维数能取得较小的值但那里使用的方法存在一定的缺陷，详细评述见文献 1 8 】从计算数学的角度来说一般直接从方程本身出发( 时域方法) 进行数值求解，能 得到准确的数值结果我们知道对于一般的系统有极限环，空间维数至少为2 ，相应 地，对于分数阶系统是否能找到有极限环的有效维数呢? 这里，我们应用数值跟踪 原理和上节预估校正算法数值跟踪分数阶b r u s s e l a t o r 系统有极限环的有效维数 分数阶b r u s s e l a t o r 系统的一般形式 c d 器t x l ( t ) = a 一( p + 1 ) x l + z z 2 c d 鬯j z 2 ( 亡) = # x a z ；z 2 ， 当p 1 = p 2 = 1 ，它是整数阶b r u s s e l a t o r 系统由已知的结论我们知道当p a 2 + 1 ， 此古典系统有唯一的极限环，而且当( a 一1 ) 2 0 ) q ，o t r 很明显函数空间以是线性空间，且满足如果p q 则仅c 岛 定义3 2 2 函数f ( t ) ( t 0 ) 满足，( m ) g ，则称f ( t ) ( t 0 ) 叼，仇z + u _ o ) 引理3 2 3 如果m 一1 o rnl詹 七= o 接下来我们考虑如下形式的时间分数阶偏微分方程的求解过程 c d o n , t u ( x ，t ) = f ( u ， 钍( 南( z ，0 ) = c 七( z ) ， 方程( 3 2 1 ) 等价表示为 c d o n , u ( x ，t ) = l u ( x ， 让( 七( z ，0 ) = c 七( z ) ， t k ，u z z ) + g ( x ，t ) ， m 一1 o l m ， 后= 0 ，1 ，m 一1 t ) + n u ( x ，t ) + g ( x ，t ) ，m 一1 q o ，盼o 、尸， 如果q = 1 ，g = 0 ，则解析解为 u ( x ，秒，t ) = - s i n ( x + y ) e 一2 以，v ( x ，y ，t ) = s i n ( x - t - y ) e 一2 以， 这与一般的n a v i e r - s t o k e s 方程的一致上述解的情形我们在图四和图五中描述 如果我们将初值条件给定为 u ( x ，y ，0 ) = 一e ( 舛圳，v ( x ，y ，0 ) = e x + 掣 ( 3 2 8 b ) 同样的计算过程我们得到 o o u ( x ，y ，t ) = 让n 时虬禹一z ：c 褊+ 揣揣+ 一训( 2 + 志， 。 n u 脚 = 幻 y zu 2 0 0 8 上海大学硕士学位论文 2 0 v ( x ，秒，t ) = n = o 剐：二南叫：c 焉 = e 忸切玩“2 吮q ) 一采五 + 揣i - - - i - 荒川+ 蔬i 南葡i 而- 一 当q = 1 ，g = 0 得到与一般的n a v i e r s t o k e s 方程一致的解 u ( x ，y ，t ) = 一e ( 蚪v + 2 扪， v ( x ，y ，亡) = e ( 蚪什2 川 图六为初始值为( 3 2 8 b ) q = 0 5 ，t = 0 0 5 时的情形 2 0 0 8 上海大学硕士学位论文 2 1 a q = 1 ，u ( x ，y ，t ) b n = 0 5 ，u ( x ，y ，t ) e q = o 1 ，u ( x ，y ，t ) 图四分数阶n a v i e r s t o k e s 方程( 3 2 8 ) 初值为( 3 2 8 a ) ，( a ) ，( b ) ，( c ) 中参数t = 3 ，l ，= 0 5 ，9 = 0 2 0 0 8 上海大学硕士学位论文 2 2 a q = 1 ，v ( x ，! ，) b a = 0 5 ， ( z ，y ，t ) c 口= 0 1 ，t ，( z ，y ，t ) 图五分数阶n a v i e r - s t o k e s 方程( 3 2 8 ) 初值为( 3 2 8 a ) ，( a ) ，( b ) ，( c ) 中参数t = 3 ，= 0 5 ，夕= 0 2 0 0 8 上海大学硕士学位论文2 3 u v y a 口= 0 5 ，u ( z ，y ，) y b n = 0 5 ，t ，( z ，t ) 图六分数阶n a v i e r - s t o k e s 方程( 3 2 8 ) 初值为( 3 2 8 b ) ，( a ) ，( b ) 中参数t = 0 0 5 ，i ，= 0 5 ，g = 0 第四章基于a d o m i a n 技巧的分数阶微分方程的数值计算 本章第一节将用a d o m i a n 分裂方法应用于a d a m s - m o u l t o n 公式得到计算分数 阶微分方程的一种新的行之有效的格式，利用数值跟踪原理和新的数值算法寻求分 数阶微分方程产生复杂动力学行为的有效维数第二节给出基于a d o m i a n 技巧的 分数阶偏微分方程的数值计算 4 1 基于a d o m i a n 技巧的数值计算方法 4 1 1a d o m i a n 技巧的数值计算方法介绍 基于第二章的预估校正思想，我们重新考虑算法的结构，观察( 2 1 4 ) ，如果 f ( t ，z ( ) ) 是非线性函数，则对应的格式是隐格式，第二章构造的格式是预估校正格 式，当然还有一种直接的方法，就是用牛顿迭代方法，但是这种方法的稳定性并没有 预期的好，原因是初值要选取得非常的好；另一方面，若讨论的是方程组，则涉及 到逆矩阵的计算，计算量会大大增加这里我们将用a d o m i a n 方法来处理( 2 1 4 ) ， 具体如下：首先( 2 1 4 ) 可以等价表示为 础州) = 揣啾l x 州) + 邢州) ) + i 轰1 篆+ 高静州地一硼 ( 4 1 1 ) 其中 ln a + 1 一( 佗一q ) ( 佗+ 1 ) q ， 歹= 0 ， a j , n + l 。1n - 歹+ 2 ) q + 1 + ( 佗一j ) q + l 一2 ( n - 歹+ 1 ) q + 1 ，1 歹仡 l x ( t n + 1 ) 是，( 亡时1 ，z ( 厶+ 1 ) ) 关于z ( 亡n 十1 ) 的线性项， n x ( t 竹+ 1 ) 表示，( 芒n + 1 ，x ( t n + 1 ) ) 关于z ( 如+ 1 ) 的非线性项根据a d o m i a n 方法的基本原理，z ( 如+ 1 ) 可以用无穷级数 形式表示为 z ( t 卅1 ) = 研( 如+ 1 ) ， ( 4 1 2 ) 其中 n 1 - 1 ，+ 知厶q ” 州机) 2 三搿磊+ 可若酉卅，f ( t j 一如) ) - 2 0 0 8 上海大学硕士学位论文 2 5 问题归结到求解( t n + 1 ) 将( 3 1 2 ) ，( 3 1 3 ) ，( 3 1 4 ) 代入( 4 1 1 ) ，我们得到 z r ( t n + 1 ) = x o ( t n + 1 ) +蒜羔2 ( 三( 量 r ( a + ) r 、台 o 。 ( 如+ 1 ) ) + ( 由此，我们可以将x r ( n + 1 ) ) ，7 = 0 ，1 ，2 表示为 其中 f q l 1 x o ( t n + 1 ) = x l ( t n + 1 ) = x 2 ( t t l + 1 ) = 研+ 1 ( t n + 1 ) a o = n ( x o ( t n + 1 ) ) ， a t 砉嘉 ( p 1 + 2 p 2 + - r p e = - r r = l ，2 ， k = o ( 南) t 七 粕。丽+ 九q r ( q + 2 ) 九a r ( q + 2 ) h d a r ( q + 2 ) n j = o ( l x o ( t n + 1 ) + a o ) ， ( l x l ( t n + 1 ) + a a ) ， r ( q + 2 ) a z i ( t n + 1 ) ) 】a ：o r ；0 研( 如+ 1 ) ) ) ，n + t 厂( 如，z ( 巧) ) ， ( l x ，( t n + 1 ) + a r ) ， ( z 1 ( k + 1 ) ) p 1 ( x 2 ( t n + 1 ) ) p 2 p l ! ( x r ( 如+ 1 ) ) 阶 肼! ( 4 1 3 ) ( p l + p 2 + 一p r ( z o ( ，l + 1 ) ) ， r 在计算中我们根据实际要求取有限项的和z ( 如+ 1 ) = x r ( 计1 ) 来近似表示 r = o z ( t n + 1 ) 由( 4 1 3 ) 和( 4 1 2 ) ，我们得到一种新的迭代算法来计算( 4 1 1 ) 这个迭代算 法是收敛的a d o m i a n 方法应用到a d a m s - m o u l t o n 方法的误差分析，我们可以由文 献 9 和定理( 4 1 ) 直接得到在这个过程中误差由两部分组成，一部分是将( 2 1 2 ) 应用复合梯形公式来计算积分得到a d a m s - m o u l t o n 公式( 4 1 1 ) ，这部分的误差在文 咖 铷 2 0 0 8 上海大学硕士学位论文 2 6 献 9 】中给出，第二部分的误差来自孑( t n + 1 ) = r = 0 截断误差，这部分的误差在定理( 3 1 ) 中给出 x r ( 亡n + 1 ) 来近似表示x ( t n + 1 ) 产生的 若我们求解方程组( 2 2 1 ) ，相应的计算过程如下；首先方程组( 2 2 1 ) 用梯形公 式离散得 其中 x l ( t t l + 1 ) = + x 2 ( t n + 1 ) = + z 。( 亡n + 1 ) = + 危0 1 f ( a l + 2 ) h q l f ( a a + 2 ) 危q 2 f ( a 2 + 2 ) h n 2 f ( a 2 + 2 ) h q 。 r ( a 。+ 2 ) 九 r ( q 。+ 2 ) r o , 1 1 ( f l ( t n + 1 z ( t 1 ) ) ) + j = o k l = o a l , j ，竹+ 1 ( 如，z ( 如) ) ， ( 厶( 亡n + l ，x ( t n + 1 ) ) ) + j = o 口2 1 1 k 2 = o a 2 , j ，n f + 1 ，2 ( 巧，z ( 岛) ) ， z 驴1 蠢 舻篙 f q 。卜l ( 厶( n + 1 z 。( t n + 1 ) ) ) + n j = 0 n 。曩n + l f s ( t j ，z ( 岛) ) k 。= o ，( k ) t k 可 ( 4 1 4 ) i 竹啦+ 1 一( 佗一q i ) ( 佗+ 1 ) 似， 歹=