（计算数学专业论文）分数阶常微分方程的高阶多步法和变分数阶扩散方程的数值方法.pdf

摘要 摘要 分数微积分的出现已有3 0 0 多年的历史，它的应用领域很广，包含各种材料的 记忆、力学和电特性描述、地震分析、电力分形网络、分数阶正弦振荡器、机器 人、电子电路、电解化学、分数电容理论、电极电解质接口描述、分形理论，特 别是描述自相似和多孔结构的动态过程、分数阶控制器设计、弹粘性系统和柔 软构造物体的振动控制、分数阶生物神经元和概率论等。分数阶微分方程的特 点是含有非整数阶导数，能非常有效的描述各种各样的物质的记忆和遗传性质， 在工程，物理，金融，水文等领域发挥越来越重要的作用。 这篇文章主要由下面几个部分组成。 绪言部分介绍了关于分数阶微积分的一些预备知识，给出了分数阶微积分 一些基本定义和性质。 接下来的第二章中，首先从基本的分数阶常微分方程出发，对l u b i c h 提出的 一个关于分数阶导数的高阶近似，将其应用于分数阶微分方程，构造高阶数值差 分格式来进行分数阶微分方程的数值求解，并在理论上给出这一算法的误差分 析，证明了它的相容性，收敛性和稳定性。 第三章对于一个推广到分数阶的松驰方程，直接利用g r i i n w a l d - l e t n i k o v 分 数阶导数定义进行离散，得到分数阶松驰方程一个数值方法，并给出了相容性， 收敛性和稳定性的证明。 在第四章中，进一步的考虑更复杂的非线性分数阶常微分方程，同样利用的 是l u b i c h 提出分数阶导数的高阶近似，构造相应的数值格式，并给出这一算法的 误差分析，即相容性，收敛性和稳定性的证明。 第五章考虑变分数阶的微分方程，在近来提出的一些模型中，分数阶导数的 阶数会随着时间或空间的变化而变化，因此在最后一章中我们讨论基于r i e s z 分 数阶导数的一类变分数阶扩散方程，给出求解这样一个方程的一个数值方法，并 对其相容性，收敛性和稳定性进行了证明。 关键词：分数阶方程；高阶多步法；变分数阶； a b s t r a c t a b s t r a c t f r a c t i o n a lc a l c u l u sh a sah i s t o r yo f3 0 0y e a r s ，i t s a p p l i c a t i o n i s v e r y b r o a d ，i n c l u d i n gt h em e m o r yo fm a n yk i n d so fm a t e r i a l ，c h a r a c t e r i z a t i o no fm e c h a n - i c sa n de l e c t r i c i t y , e a r t h q u a k ea n a l y s i s ，r o b o t s ，e l e c t r i cf r a c t a ln e t w o r k ，f r a c t i o n a ls i n e o s c i l l a t o r , r o b o t , e l e c t r o n i cc i r c u i t s ，e l e c t r o l y s i sc h e m i c a l ，f r a c t i o n a lc a p a c i t a n c et h e - o r y , e l e c t r o d ee l e c t r o l y t ei n t e r f a c ed e s c r i p t i o n ，f r a c t a lt h e o r y , e s p e c i a l l yi nt h ed y n a l i l i c p r o c e s sd e s c r i p t i o no fp o r o u ss t r u c t u r e ，f r a c t i o n a lc o n t r o l l e rd e s i g n ，v i b r a t i o nc o n t r o l o fv i s c o e l a s t i cs y s t e ma n dp l i a b l es t r u c t u r eo b j e c t s ，f r a c t i o n a lb i o l o g i c a ln e u r o n sa n d p r o b a b i l i t yt h e o r y ，e t c t h ec h a r a c t e r i s t i co ff r a c t i o n a lo r d e rd i f f e r e n t i a le q u a t i o ni s c o n t a i n i n gt h en o n - i n t e g e ro r d e rd e r i v a t i v e i tc a l le f f e c t i v e l yd e s c r i b et h em e m o r y a n d t r a n s m i s s i b i l i t yo fm a n yk i n d so fm a t e r i a l ，a n dp l a y sa ni n c r e a s i n g l yi m p o r t a n tr o l ei n e n g i n e e r i n g ，p h y s i c s ，f i n a n c e ，h y d r o l o g ya n do t h e rf i e l d s t h i st h e s i sc o n s i s t so ft h ef o u rc h a p t e r s i n t r o d u c t i o ng i v e ss o m ec o n c e r n i n gf r a c t i o n a lc a l c u l u st op r e p a r et h ek n o w l e d g e a n dp r e s e n tb a s i cd e f i n i t i o n sa n dp r o p e r t i e so ff r a c t i o n a lc a l c u l u s i nc h a p t e r2 ，s t a r t i n gf r o mt h eb a s i cf r a c t i o n a lo 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 ，w e a p p l yah i g ho r d e ra p p r o x i m a t i o no ff r a c t i o n a ld e r i v a t i v ea d v a n c e db yl u b i c ht of r a e t i o n a ld i f f e r e n t i a le q u a t i o n ，c o n s t r u c tah i g hn u m e r i c a ld i f f e r e n c es c h e m et os o l v et h e 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 ，p r e s e n te r r o ra n a l y s i so ft h ea l g o r i t h m st h e o r e t i c a l l y , a n d p r o v et h ec o n s i s t e n c y , c o n v e r g e n c ya n ds t a b i l i t y i nc h a p t e r3 ，c o n s i d e r i n gf r a c t i o n a lr e l a x a t i o ne q u a t i o n ，w em a k eu s eo fd i r a y t h eg r u n w a l d - l e t n i k o vd e f i n i t i o nt od i s c r e t ef r a c t i o n a ld e r i v a t i v e ，o b t a i nan u m e r i c a l m e t h o do ff r a c t i o n a lr e l a x a t i o ne q u a t i o n ，a n dg i v et h ep r o o fo fc o n s i s t e n c y , c o n v e r - g e n c ya n ds t a b i l i t y i nc h a p t e r4 ，w ec o n s i d e rm o r ec o m p l e xf r a c t i o n a ln o n l i n e a rd i f f e r e n t i a le q u a t i o n ，a l s ou s i n gt h eh i g l lo r d e ra p p r o x i m a t i o np r e s e n t e db yl u b i c ht oc o n s t r u c tc o r r e s p o n d - i n gn u m e r i c a ls c h e m ea n dg i v i n gt h ee l t o ra n a l y s i so f t h ea l g o r i t h m s i nc h a p t e r5 ，n o t i c i n gt h a tr e c e n t l yi ns o m em o d e l st h eo r d e ro ff r a c t i o n a ld e r i v a - t i v ec a nv a r i e t yw i t ht i m eo rs p a c e ，w ed i s c u s s e dv a r i a b l eo r d e rf r a c t i o n a ld i f f u s i o n e q u a t i o nb a s e do nr i e s zf r a c t i o n a ld e r i v a t i v e w ep r e s e n t e dan u m e r i c a lm e t h o dt o s o l v et h i sk i n do fe q u a t i o n ，a n dp r o v ei t sc o n s i s t e n c y , c o n v e r g e n c ya n d s t a b i l i t y k e yw o r d s ：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 ；h i g ho r d e rm u l t i p l em e t h o d ；v a r i a b l e o r d e rf r a c t i o n a ld i f f u s i o ne q u a t i o n ； l i i 厦门大学学位论文原创性声明 兹呈交的学位论文，是本人在导师指导下独立完成的研究成果。 本人在论文写作中参考的其他个人或集体的研究成果，均在文中以 明确方式标明。本人依法享有和承担由此论文产生的权利和责任。 声明人( 签名) ：撇 w 矸年 6 月勺e l f 厦门大学学位论文著作权使用声明 本人完全了解厦门大学有关保留、使用学位论文的规定。厦门大 学有权保留并向国家主管部门或其指定机构送交论文的纸质版和电 子版，有权将学位论文用于非赢利目的的少量复制并允许论文进入 学校图书馆被查阅，有权将学位论文的内容编入有关数据库进行检 索，有权将学位论文的标题和摘要汇编出版。保密的学位论文在解密 后适用本规定。 本学位论文属于 1 、保密( ) ，在年解密后适用本授权书。 2 、不保密( 作者签名：刊影陶。 导师签名：卅亿协丢 日期：砷年6 月1 e 1 日期：岬年6 , 97 日 第一章绪论 第一章绪论 1 1引言 分数微积分理论是数学分析的一个分支，是专门研究任意阶积分和微分的 数学性质及其应用的领域。最早提出这一思想的是数学家l e i b n i z ( 1 6 9 5 年) ，在 它之后有e u l e r ( 1 7 3 0 年) 、l a p l a c e ( 1 8 1 2 年) 、f o u r i e r ( 1 8 2 2 年) 、a b e l ( 1 8 2 3 1 8 2 6 年) 、l i o u v i l l e ( 1 8 3 2 1 8 7 3 年) 、r i e m a n n ( 1 8 4 7 年) 等先后在该领域做出了杰出的贡 献。到本世纪初n u t t i n g 1 1 ( 1 9 2 1 年) 、g e m a n t 2 】f 3 1 ( 1 9 3 6 年，1 9 3 8 年) 、s c o t t b l a i r 4 1 ( 1 9 4 9 年) 将分数微积分理论引入粘弹性材料的本构关系；同时，前苏联著名科 学家r a b o t n o v s 1 6 ( 1 9 4 8 年，1 9 8 0 年) 在介绍它的具有弱奇核遗传效应固体力学 时也暗示了需要分数微积分理论；另外，在早期的研究工作中c a p u t o 7 ( 1 9 6 6 年) 、m a i n a r d i 8 1 9 10 9 7 1 年) 等将分数微积分理论在有关地震、冶金等方面也得到 很好的应用。 近三百年以来，分数阶微积分或分数阶演算这一重要的纯数学分支已 渐成体系。但是，对于绝大多数工程技术界学者而言，它还鲜为人知。直 至m a n d e l b r o t 1 0 提出分形学说，将r i e m a n n - l i o u v i l l e 分数阶微积分用以分析和 研究分形媒介中的布朗运动以来，分数阶微积分才在许多学科的现代工程计算 中得以广泛关注和应用。在近十几年来，这一理论得到了相当的关注【l 5 1 ，像在 黏弹性材料 1 9 - 2 2 1 ，电气化学过程c 2 3 1 ，电介质极化现像 2 4 1 ，色噪声 2 5 1 ，不规则 扩散，信号处理【凋，控制理论【1 5 】，自然界中多孔或裂断介质中溶质的反应和扩 散 2 7 , 2 8 以及混沌【2 9 】等现像的描述中得到了诸多应用。主要是因为在这些物理 化学现像中，相对于原有的整数阶模型，用新的分数阶模型描述的结果会更精 确。同时，对于整数阶微分方程，相关数值算法理论已比较成熟，而对于这些分 第一章绪论 数阶模型中的分数阶微分方程，数值算法研究起步不久，特别是理论分析方面目 前还比较有限 3 0 - 3 5 。 这篇论文主要分为以下几个部分，引言部分是关于分数阶微积分的一些预 备知识，给出了分数阶微积分一些基本定义和性质。 第二章讨论分数阶常微分方程。由于许多分数阶方程可以化为等价的阶数 在0 ，l 之间的分数阶方程组来求解，因而首先从基本的分数阶常微分方程出发， l l o o , a y ( t ) = 厂( f ) ，0 o b l ， r 【o ，丁】， i ) ，( o ) = 0 ， 利用l u b i c h 提出的一个关于分数阶导数的高阶近似， np o d t a y ( t ) h 川o f l a _ ；y ( t j ) + h 川w j y ( t j ) j = 0j = l 其v p h = t n ，0 = j h ，p 表示方法的阶数，即近似阶数为d ( 胪) ，同时也等于需要的 初始点的个数。可以看到，上面的数值格式是由系数( 武口) ：七= 0 ，l ， ， 咐： j = l ，2 ，p ；k = l ，2 ，) 确定的。这些系数同导数阶数口和方法的阶数p 有关， 它们的定义和性质在后面详细给出。共给出t p = 2 ，3 ，4 ，5 ，6 时近似格式系数的 求法。然后将其应用于分数阶微分方程，构造高阶数值差分格式来进行分数阶 微分方程的数值求解，并利用线性代数的克莱姆( c r a m e r ) 法则，通过求矩阵行列 式的方法，得出数值解和精确解之间的误差表达式，在理论上给出了这一算法的 误差分析，证明了它的相容性，收敛性和稳定性。 第三章中，对于一个推广到分数阶的松驰方程【1 5 】， f lo d f y ( t ) + a o ) y o ) = 厂o ) ，0 0 。直接利用g f i i n w a l d l c m i k o v 分数阶导数定义进行离散，得到分数 2 第一章绪论 器。一扩。 掣= 毗f ) ) 如卅m ，砒 1 o ，胗o 。 关于分数阶微积分的定义，目前有许多不同的表达形式，主要有下面这么几 类。 1 2 分数阶积分 口是一个正实数，令r l 一1 o t n ，那么，一个定义在【口，6 】区间上的函 数，( f ) 的a 阶分数阶积分的定义是 儿) = 溉警棚m ) ， 口町a 儿) = 溉 a 量引饨一曲) ， 4 第一章绪论 它还有另一个表达彤式 口d 7 口九) = 丽1r ( ，叫弘1 m ) 把 其中r ( z ) 表示通常的g 锄m a 函数，即 r ( z ) 2f o e 一t z - 1 出。 已经证明了这两个形式是等价的【1 5 1 。 1 3r i e m a n n l i o u v i l l e 分数阶导数 对于一个正实数a ，令n 一1 口n ，一个定义在陋，翻区间上的函 数，( f ) 的口阶慰锄锄n - “o u v i l l e 分数阶导数的定义是； 口d t a f ( 归：a n ( 志r 喾尚) 。 1 4g r o n w a l d l e t n i k o v 分数阶导数 同样的，对于一个正实数a ，令n 一1 口刀，一个定义在妇，纠区间上的函 数，( f ) 的口阶g 釉n w a l d k ：t n ：蛔分数阶导数的定义是 口钟饨) = l i 枷r ai l 一口簖，( f ) = 2 觋 一a ( t - 毛a ) h ( 一1 ) ，( ；) ，( a + t - j h ) ， 其中( 擘) 表示二项式系数，即 a 、a ( a 一1 ) ( a j + 1 ) t j ) 5 了一。 戤e m 锄- l i o u v i l l e 分数阶导数和g r i i n w a l d - l c m i k o v 分数阶导数之间有这样 的一个等价关系，对于正实数饯，n 一1 a n ，如果定义在a ，纠区间上的函 数，( f ) 有直到万一1 次连续导数，并且厂( 矗) ( f ) 在【口，纠是可积的，那么这时r i e m 锄n 。 l i o u v i l l e 分数阶导数和g m n w a l d l e t n i k o v 分数阶导数是等价的【1 5 】。 5 第一章绪论 1 5 其他的一些定义 c a p u t o 分数阶导数定义。对于正实数a ，r - - l 口n ，一个定义在 口，纠区 间上的函数厂( f ) 的o t 阶c a p u t o 分数阶导数的定义是 口d a y ( r ) = 币l ， 厂( 一) ( f ) d f 瓦丽雨丽。 有时函数需要定义在无穷区间上，所以有这样的分数阶导数的定义， 一。砰饨) = 币1 丽 厂( n ) ( f ) d f 瓦可雨丽。 c a p u t o 分数阶导数和r i c m a n n l i o u v i l l e 分数阶导数也有一个关系。对于正实 数口，n l 口 l0 ，0 f 口 口班( ) 巾) = 科“f ) ， 口 1 1 0 , 0 f 口 其它一些函数的r i e m a n n l i o u v i l l e 分数阶导数， 口d 尹厂口e l ，l a ( 九f ) 7 第一章绪论 口d 尹c o s h ( 瓜) = 厂a e 2 ，l a ( 九f 2 ) n as i n h 亓( v 佤一 t ) = t l - a 助，2 一口( 九f 2 ) 、l t 1 6 2 分数阶算子的复合运算 分数阶算子在一定条件下可以进行相互的复合，一个定义在【口，纠区间上 的函数，( f ) ，对于分数阶积分来说，分数阶积分之间的阶数是可以直接相加的， 令口，卢 0 ， 0 0 7 口( 口矿饨) ) = 口町州饨) 。 而对于分数阶导数，这里使用r i e m a n n 1 i o u v i i l e 分数阶导数定义，有这样一 些关系，令n l 0 ，则先求积分 再求导数时 口钟( 口d 卢巾) ) = 口群邛饨) 先求导数再求积分， 口玎砰饨) ) = 口群- 卢弛) - ，_ a d t a - j f ( f ) ，剐高差尝击 接下来，当整数阶导数和分数阶导数复合运算时， 历d k ( 口d t a f ( f ) ) = 口矽a 竹) 8 第一章绪论 口砰( 掣) = 硝州一k 岛- , 可# j ) ( a 丽) ( t - a 可) j - 口- k 如果是分数阶导数之间的复合运算，令n 一1 0 ，那么有 4 钟( 口钟巾) ) = 口群邯弛) _ 削 a d t a - j f ( r ) ，剐而( t - f a ) 两- # - j y 1 8 3 分数阶导数的积分变换 下面介绍分数阶导数的积分变换。通常积分变换是可以用来解析求解整数 阶微分方程，对于分数阶导数来说，各种积分变换也是解析求解分数阶微分方程 的重要手段之一。 首先是分数阶导数的l a p l a c e 变换【1 4 1 1 3 】。一个函数的l a p l a c e 变换是 f ( s ) = 厂( f ) ；j ) = e - s t f ( t ) d t 对于一个函数的r i e m a n n 1 i o u v i l l e 分数阶导数，它f l 勺l a p l a c e 变换是 l 口砰饨) ；j ) = ，f ( 墨) 一n 毛- i 口群一饨) ，剐 其o e n - l = e i d x g ( t ) d t 可以看到，对一个函数g ( f ) 进行f o u r i e r 变换，要求g ( f ) 在整个实轴上有定义，所以 采用无穷区间上的分数阶导数定义是合理的。因此分数阶导数定义是 一捌比) = 币丽1 g ( 一) ( f ) d f f 丽雨丽 其中万一1 口n 。这时分数阶导数的f o u r i e r 变换就是 b 一。z 垆g ( f ) ；) = ( - i o ) a g ( ) 。 9 第二章分数阶常微分方程的高阶数值方法 第二章分数阶常微分方程的高阶数值方法 对于整数阶常微分方程的数值解法，如欧拉法、线性多步法等都已有较完 善的理论。而对于分数阶微分方程，数值方法和误差估计的理论研究相对较少。 在这章中，我们考虑最基本的分数阶常微分方程，这一章方程的分数阶导数的阶 数是0 a l 的，许多更高阶的方程常常可以化为阶数小于或等于l 的方程组来 处理，因此首先研究阶数在0 ，1 2 间的分数阶方程。对这一方程，我们引进分数 阶的线性多步法，导出了分数阶常微分方程初值问题的高阶近似，证明了其方法 的相容性和收敛性，并且给出稳定性分析。最后给出一些数值例子，证实了这个 分数阶线性多步法是解分数阶常微分方程的一个有效方法 2 1 分数阶常微分方程 我们考虑如下的分数阶常微分方程 t；od竺ta)(yi(t)三-f。(t，)，。口1 c 2 l ， 定义2 i ： 科巾) = 历d t l ( 南r 若尚) , n - - l a n 还有一个与之等价的c 黼n w a l d k t n i l 【o v 定义 定义2 2 ： 1 0 第二章分数阶常微分方程的高阶数值方法 。钟饨) = l i 川m a 川簖巾) = l i m 。 川笺( 一) ，( ；) 厂( f 一弘) ， 2 2 分数阶导数的离散 首先我们考虑对一个已知函数求分数阶导数。对一个正实数a ，( m 一1 a m ) ，当函数y ( f ) 在区间【o ，t 】中有。至m 一1 阶连续导数，且) ，( 坍) ( f ) 在【o ，t 】中 可积时，已经证明了r i e m a n n l i o u v i l l e 定义的分数阶导数和g r i i n w a l d l c t n i k o v 定 义的分数阶导数两者是一致的。这就相当于给出了对分数阶导数的一个近似 方法，因为很自然的，将g r i i n w a l d l e t n i k o v 定义中的极限号去掉，就可以得 到由有限个点表示的分数阶导数的离散形式。这实际上就给出了一个1 阶方 法，l u b i c h 又进一步提出y 2 6 阶形如下面形式的分数阶导数的近似方法 4 2 1 。 科y ( f ) 。蹲y ( f ) = 一口童! ( 巧) ( 2 2 ) i - - o 其中j l = t n ，t j = j h ，系数 ；口) ；j = o ，i ，2 ) 是一个生成函数的泰勒展开系 数，即 矗口) + d 口k + 吐a ) p + = w 参口( 曲 阶数p 不同，生成函数也不同，下面这些函数分别对应2 至6 阶方法的生成函数。 疃哟= ( 主一及+ 三p ) 口， 叫哟= ( 芸一缸+ 三一三，) 口， 以o ) = ( 篙一缸+ 辞一三p + 丢，) a ， 蜘( 等一缸+ 好一萼，+ 三一扣， 第二章分数阶常微分方程的高阶数值方法 咙口) ) = ( 等一缸+ 萼p 一了2 0 工1 十百1 5 工4 一亏6 ，+ 丢) a 。 事实上，上面提到的i 扫g r i i n w a l d l e t n i k o v 定义得到的近似方法就是 由w f 口g ) = ( 1 一工) 口的泰勒展开系数 ；口) = ( 一1 ) ( 了) ；j = 0 ，1 ，2 ) 得到的近 似方法。另外，a o 时，由上述函数展开系数得到的公式( 2 2 ) 实际上是分数阶 积分的近似公式。本文只考虑分数阶导数，但仍有用到口 0 ，p 对应方法的阶数( 2 6 ) 。我们可以发现对于固定的j 9 l ，即使p 提高， 误差阶也只有d ( 胪) 。所以，为了有更高阶的格式，我们利用l u b i c h ( 1 9 8 6 ) 提出的 技巧 4 2 1 ，用下式来近似分数阶导数 o d t a y ( t ) h 川t o ( n a _ ；y ( t j ) + h 川w ( o ) ( 2 3 ) 其思想就加上一定的修正项使之能去掉d ( 胪) 这一项，这样，总的误差阶就只 有o ( h p ) t 。 ( 2 3 ) 中的f 是由下列线性方程组得到的 善p 咐尸= 揣矿口一善p ! 删，一1 ) 这一格式的近似阶数在后面由引理给出。 2 3 分数阶常微分方程的数值方法 程 我们把上节中分数阶导数的近似方法应用于下面最基本的分数阶常微分方 o d t a y ( t ) = ，( f ) ，y ( o ) = 0 1 2 第二章分数阶常微分方程的高阶数值方法 其中f ( 0 ，列，0 p 时，氪由下式求出 j l a 西口员+ j l 一口k - i 畦竺；y ( 巧) + 一a pm 拶( 。) = ， 可得 j l 一口毋孤一j l 一口毋y 做) + j l a k 竺( 。) + j l 一哇帆( 。) = 厶， j - - - oj = t h - a 矗恹一) ，( 块) ) + o d ? y ( t k ) + o ( h p ) = f k ， 再由o d t y ( t k ) = f ( t k ) = f k ，有 h - a 武口慨一y ( 玖) ) = o ( h p ) ， 即 一y(tk)=o(hp+yk y ( t k 口1 。 一 ) =“) o 口 现在我们给出总体截断误差。同样的，分数阶线性多步法的总体截断误差 定义茭jf h y o ，y l ，强一1 求得的y 七与) ，陬) 之间的误差。 1 4 第二章分数阶常微分方程的高阶数值方法 定理2 2 - 对于分数阶微分方程( 2 1 ) ，由分数阶线性多步法( 2 4 ) 求得的强的总体 截断误差为o ( h p ) 。 h 川趟竺加+ j l 叫t - m = f k 。 鼻 一 j = o，= 1 一口k4 竺；y ( 。) + 一噬( 。) + d ( 矿) ： 。 j l 一口k 趟- 口j 勺= o ( h p ) 。 曰：坩_ 二篓置、， 1 5 第二章分数阶常微分方程的高阶数值方法 f f e k i 一 一l p 二：、， 其中言：( d ( 胪) ，o ( h p ) ) r 。对这一线性方程组，我们只须解出纨。易求b 的 行列式吲= o 一口武口) 七- p ，由克莱姆法则，还需求下列矩阵的行列式 二篓； 为计算i i ，记a i = h 按最后一行展开，可以得到 i = ( 一1 ) 扣p + 1 l d la 2 a oa l 0o + ( 一1 ) 七一p + 七一p 口0 1 6 h - a t o a _ ) ( p + 1 ) 俨摧1 ) 一 j l a 武a ) 鲰一( p + 1 ) a ( k 1 ) 一( p + 1 ) c 礓一( p + 2 ) a c k 1 ) 一q + 2 ) 惰 州 忡 知 肭 啡 廿 知 柚 p l 砚 内 0 k 坟 口 ： 雄 ，= 气 l 口， 一 一 一 铆 如 0 l 2 坟 卜 一 外 翻 盼 缈 第二章分数阶常微分方程的高阶数值方法 由引理2 1 ，设 则可推出 ( a o + a l x + a 2 x 2 + ) ( c 1 d + c 1 工+ c 2 + ) = 1 。 = b p + l a 七。一p q 一1 ) + 口o = b p + l a ：- p q 一1 ) + 口o ( 一1 ) 七一1 ) + l 2 + 碚 圾口l 巩一1 a o 圾口l 鲰一( 升2 ) 坟一1 a o a ( k 1 ) 一( p + 2 ) 2 0 口l 口2 a k 一0 + 2 ) 口。口l a ( k 1 ) 一( p + 2 ) 00 a k 一( p + 3 ) a ( k 1 ) 一( p + 3 ) + 3 0 a o = b p + l 奇p c 七一1 ) + 2 古p q 一( p + 2 ) + 碚 如此继续则有 a l 圾口1 坟一1 a o a k 一0 + 3 ) 口( 七一1 ) 一0 + 3 ) b p + 3 0 知 l = b p + l a k o - p c k 一1 ) + 2 口o k - p c k 一( 舛2 ) + 1 7 第二章分数阶常微分方程的高阶数值方法 + 圾一2 罐一p c 2 + 古p + 2 ) 坟口1 坟一1a o = b p + l a ：- p 铅一( 肘1 ) + b + 2 口：一p q 一( 舛2 ) + + 玩一2 罐一p c 2 + a ：一( 舛2 ( 一圾一l 口l + a o b k ) = b p + j a i l - p c 七一( p + 1 ) + 易p + 2 口：一p c 七一( p + 2 ) + + b k 一2 口o k c 2 + b k 一1 k - p ( 一口l a 0 2 ) + b k a k o - p a 0 1 = 古p ( 坼l q 一1 ) + + 2 c k 一( p + 2 ) + + b k l c l + b k c o ) 。 所以再由吲= ( 一口西a ) 卜p = a k o - p 推出 推得 e k = i b i i b i = ( 匆舛l c 七一( p + 1 ) + b p + 2 c k 一( p + 2 ) + + 玩一i c l + b k c o ) ， 又因为 所以 l e k l _ m 攀 i b i l ( 1 0 ，l e v i c l i 口一，由此得到蚓= i h 口i c l h a i a 。因 此再f h n h = t ，有 e k l m a x n l b i l c i t a n - ah a - ! + 0 一1 ) 口一1 + + 2 口一1 + l a 一1 ) 。 1 8 第二章分数阶常微分方程的高阶数值方法 注意到 n a ( n a 一1 + 一1 ) 口一+ + 2 a 一+ 1 一1 ) =l(n)a-1+n ( 孚) 弘1 + + ( 署) a - i + ( 圹) 而 ，l i m 妻1 ( ) a - i f o -= 1 ，一l 出= 否1 。 因而易得存在常数c 2 ，使得 ：1 ( ( n ) a - , n - i - ，口一l + + ( 三) a - i + ( 三) 口一1 ) p 时，y k ，畦由下式求得， h 川竺加+ 川pm = ， 七 j = oj = t 和 川善趟竺物+ j i 川暑p 巧= 露。 七 1 9 第二章分数阶常微分方程的高阶数值方法 两式相减得到 ；句：觚。川；句= 觚。 j = p + 1 同样的由上式得到下列线性方程组 j l a 矗口) o | l a “口) 一口罐a ) o 0 h - a c o 丘( 一a ) 叶2 ) 一鹾1 ) 一( p + 3 ) 记为& = _ b ，记口f = 一口西引，设 类似的可以得到 所以 j l 一口矗a ) o h - a o o 丘( 弋a ) p + 1 ) 一1 ) 七+ 2 ) h - a o j f a ) h - a o a ) ( 知+ 口l x + 哦，+ ) ( c o + c l 工+ c 2 ，+ ) = 1 。 七= ( l 一( 舛1 ) - t - b p + 2 c k 一( 肘2 ) + + 圾一l c l + b k c o ) ， i 鼠l 1 黔 i 协i ) ( i q 一( 舛1 ) i + i q 一( 肘2 ) i + + i c l l + i c 0 1 ) 曙删) 喈 | 蝴i ) _ c i l f ( t ) l l 一。 2 0 靠 岛一i 聃2 研1 口 o o 蛳 蛳蛳 第二章分数阶常微分方程的高阶数值方法 2 5 数值例子 ，o d ? y ( f ) = ，( f ) ，te 【0 ，州， 1 ) ，( o ) = 如，o 口l 。 f ( t ) = ( 0 5 0 c + 0 5 ) 2 ，r ( 4 + o 5 a ) ，3 0 5 口3 f ( 3 + o 5 a ) 。2 0 5 口。2 f ( 2 + 0 5 a ) 1 0 5 口、 , r ( 4 - o 5 a ) 1一r ( 3 一o 5 a ) 十了i f = 丽。 ) ，( f ) = ( o 5 a + o 5 ) 2 t 3 刖口一3 t 2 + o 5 a + 力l 枷a ) 口= 0 3 时，应用分数阶线性多步法，分别取方法的阶数p = l 和p = 3 ，求得 口= 0 8 时，同样分别取方法的阶数p = l 和p = 3 ，求得的数值解和精确解的 yod。ta，y：(t)，。=，f。(t),口t_，0，丁】 2 l 第二章分数阶常微分方程的高阶数值方法 表2 1 例2 1 的数值解和解析解的表( 口= 0 3 ，p = 1 ) 最大误差是1 6 2 4 1 7 5 3 9 e - 0 0 2 ( h = 0 0 5 ) 图2 1 例2 1 的图形a = 0 3 ，p = 1 第二章分数阶常微分方程的高阶数值方法 表2 2 例2 1 的数值解和解析解的表( a = 0 3 ，p = 3 ) 最大误差是2 8 2 6 9 9 9 0 1 e 一0 5 ( h = 0 0 5 ) 图2 2 例2 1 的图形a = 0 3 ，p = 3 第二章分数阶常微分方程的高阶数值方法 表2 3 例2 1 的数值解和解析解的表( 口= 0 8 ，p = 1 ) 最大误差是3 8 1 7 9 7 0 2 0 e - 0 0 2 ( h = 0 0 5 ) 图2 3 例2 1 的图形a = 0 8 ，p = 1 第二章分数阶常微分方程的高阶数值方法 表2 4 例2 。l 数值解和解析解的表( a = 0 8 ，p = 3 ) 最大误差是6 0 6 11 0 8 3 2 e - 0 0 4 ( h = 0 0 5 ) 图2 4 例2 i 的图形口= 0 8 ，p = 3 第二章分数阶常微分方程的高阶数值方法 t 图2 50 口l 时例2 1 的图形 o123 4 567 t 图2 60 口l 时例2 2 的图形 第二章分数阶常微分方程的高阶数值方法 图2 70 口l 时例2 3 的图形 例2 3 ： 厂册( f ) ：饨【0 ，孔 ly ( o ) = ) i o ，o a 1 ( t ) = s i n ( t ) ，t = 2 石，y ( 0 ) = o 时，给出0 口l 时，解在随着阶数变化时的图 形，见图2 7 ，从中可以看到三角函数在分数阶算子作用下的性态，正弦函数在分 数阶算子作用下过渡到它的导数余弦函数。 第三章分数阶松驰方程的数值方法 3 1 引言 第三章分数阶松驰方程的数值方法 进一步的，当分数阶常微分方程为如下形式时，我们称其为分数阶松驰方 程1 1 5 1 ， o p 芦y o ) + a p 沙o ) = f ( t ) ，0 0 。这一章我们给出这个方程的一个数值方法及其误差分析。方程中 分数阶导数的定义采用r i e m a n n l i o u v i l l e 定义 。珊，= 志嘉拦尚) , n - - l a 如 并利用等价的( h u n n w a l d k 炳i l 【0 v 定义【1 5 】， 科如) = l i 川m h 一口i 毛t h i l i m h ( - 1 ) 也) 北刊科) ，( 归j i 川川毛( 。1 ) 吣沙一弘) 得到分数阶导数的一个近似公式， o d t a y ( t ) h - a 壹( _ 1 ) 、) y ( t - j h ) ，n h j = o j = r ( 一1 ) j ( 1 ) ，= i 其近似阶数为一阶 1 5 1 。 j i l l3 1 假定g ( f ) 在【o ，引上充分可微，并令4 a = ( 一1 ) j ( 了) ，则 o d t a y ( t ) = d ( 7 1 ) + j l a 砖a y o j h ) ，n h = t ，f 【o ，卅。 ( 3 2 ) 第三章分数阶松驰方程的数值方法 实际上弓口= ( 一1 ) ，( 了) 是生成函数w ( 石) = ( 1 一工) a 的泰勒展开系数，即 ( i - x ) 口= 武a + 4 口) 工+ 吐a ) + 。 因此对于方程( 3 1 ) ，令j l = r n ，0 = y h ，= 0 ，l