多项分数阶常微分方程的数值积分法.docx

多项分数阶常微分方程的数值积分法摘要 近年来，分数阶微积分在解决各种具有遗传和记忆性质的物理、材料和力学、信息领域问题中大放异彩，在建立各种数学模型的场合最为普遍使用。本文主要是通过Riemann-Liouville 分数阶积分来研究并解决多项分数阶常微分方程初值问题，比较基于一阶矩形公式的显式方法和基于二阶卷积权方法对于具体算例解的误差分析，对求解多项分数阶常微分方程的数值积分方法进行探讨。在解初值问题之前，本文首先引入了Riemann-Liouville分数阶积分的定义及相关性质以便读者理解，给出了基于一阶左右矩形公式的显式方法和基于二阶卷积权逼近法两种求解多项分数阶常微分方程初值问题的常规方法；还演示了一例基础的三项分数阶微分方程初值问题如何用以上两种方法求解，并根据不同的步长进行多次迭代比较，对比迭代过程中的收敛阶数和近似解与精确解进行误差分析，评价了两种方法的适用环境和相关优劣性。希望找到新的理论方法，以打破现有的限制条件，力求构建一套完善的分数阶微分方程理论。关键词 多项分数阶常微分方程 基于一阶矩形公式的显式方法 二阶卷积权方法 Numerical integration method for multiple fractional ordinary differential equationsAbstract In recent years, fractional calculus has been used to solve various problems with genetic and memory characteristics in physics, materials, mechanics, information and other fields, as well as to establish various mathematical models. In this paper, Riemann-Liouville fractional integration is used to study and solve the initial value problem of multiple fractional ordinary differential equations. The error analysis of the specific solution based on the explicit method based on the first-order rectangular formula and the method based on the second-order convolution weight is compared. The numerical integration method for solving multi-terms fractional-order ordinary differential equation is discussed. Before solving the initial value problem, this paper first introduces the definition and related properties of Riemann Liouville fractional integral for readers to understand. Two general methods for solving initial value problems of polynomial fractional ordinary differential equations are given, which are based on the first order left and right rectangle formula and the second order convolution weight approximation method. In addition, how to use the above two methods to solve the initial value problem of a basic three term fractional ordinary differential equation is demonstrated. The iterative process is compared according to different steps. The error analysis between the order and the approximate solution and the exact solution is carried out, and the applicable environment and the relative advantages and disadvantages of the two methods are evaluated. Hope to find a new theoretical method to break the existing restrictions and strive to build a complete set of fractional differential equation theory.Key words Multi-terms fractional-order ordinary differential equation, An explicit method based on the first order rectangle formula, Second order convolution weight method目 录引 言11 预备知识21.1 Riemann -Liouville分数阶积分的定义及相关性质21.2 左右矩形公式21.3 二阶卷积权逼近法32解分数阶常微分方程的数值方法32.1 Newton迭代法32.2基于一阶矩形公式的显式方法52.3二阶卷积权逼近法63数值算例63.1例题与解63.2基于一阶矩形公式显式方法的误差分析73.3基于二阶卷积权逼近法的误差分析10结论14参考文献15致谢16附 录A17附 录B19引 言多项分数阶常微分方程是常微分方程理论中的重要一环，它基于现实的各种实际问题，关于它的讨论经久不衰。1695年，在德国数学家Leibniz和法国数学家L'Hopital通信中分数阶微积分的思想萌芽了，它几乎是和经典微积分同时出现。像GL定义、RL定义和Caputo定义，众多学者从不同的角度切入问题，得到了不同的分数阶微积分定义。因为它的独特性质，以及在力学、物理学、材料力学、水文地理学、聚合体流变学、电子信息学等领域的成功应用，一些自然科学以及工程应用领域的非经典现象能够被刻画、再现出来。自新世纪以来，分数阶微积分的应用越来越得到人们的重视，是当前建模和数学领域中一个热点。分数阶微分算子具有非局部性，故可以用来描述具有历史记忆效应、遗传性质和空间全域相关性等动力学和物理过程。而且相比于整数阶导数构造的数学模型时不时就需要构造新模型，只因为一个材料上或外在因素导致的小小变动，分数阶导数建模十分轻松、参数物理意义简明清晰，现已成为复杂力学与物理过程数学建模的重要工具。现在分数阶微分方程在学术界还处于理论探索阶段，现今获得的成果大多是对经典微积分方程理论的简单推广，大多数分数阶微分方程还处于混沌状态，通过现有的部分理论无法求解。而且，目前该领域的研究主要集中在分数阶偏微分方程方向，反而对最经典的初值问题研究的不多。我有意针对这方面的一部分知识进行探讨，希望能从中找到新的理论方法，以打破现有的限制条件，力求构建一套完善的分数阶微分方程理论。本文的第一部分列出了本篇文章所要用到的相关知识，便于读者理解；第二部分计算出了三种求解多项分数阶常微分方程初值问题的常规方法；第三部分则演示了一例基础的三项分数阶微分方程初值问题如何用基于一阶左右矩形公式的显式方法和二阶卷积权逼近法求解，并根据不同的步长进行多次迭代实验，对比迭代过程中的收敛阶数和近似解与精确解进行误差分析，评价了两种方法的适用环境和相关优劣性。211 预备知识在这一节中，我列出了本文中所涉及的一些必要知识，为了便于读者理解。本文以最简单的多项分数阶常微分方程为基础： 1.1 Riemann -Liouville 分数阶积分的定义及相关性质定义 (Riemann -Liouville 分数阶积分) 对于一个正实数，那么，一个定义在上的函数的阶Riemann Liouville分数阶积分定义为 其中Gamma函数的定义为它具有如下两个性质： Riemann Liouville分数阶积分具有如下性质：性质1： 性质2（积分可加性）： 1.2 左右矩形公式设 且 ，则有 和 其中 。1.3 二阶卷积权逼近法对，如果 且，则 其中为常数，是由生成函数产生的系数。生成函数具体计算方法如下：2 解分数阶常微分方程的数值方法在这一节中，我给出了三种求解多项分数阶常微分方程的数值方法：1.利用基础的Riemann -Liouville 分数阶积分及其性质，求出了一般的可用如Nowton迭代法等方法求解的差项迭代公式；2.对Riemann -Liouville 分数阶积分后的方程使用左右矩形公式，求出了该法特定的迭代公式；3.对Riemann -Liouville 分数阶积分后的方程进行二阶卷积权逼近，求出了该法特定的迭代公式。2.1 Newton迭代法 将式展开得到 在式左右分别进行Riemann -Liouville 分数阶积分 由性质1得到 由性质2得到 整理得 由式可得 将t分割，令 时， 令 ，由，式可化为 将求和项重新分布得 由上式可推出 两式相减整理后得到迭代公式当 时，当 时， 令，运用牛顿迭代法（Newton's method），即可解出 的近似解。由式继续迭代，即可解出 2.2 基于一阶矩形公式的显式方法问题 可以等价地转化为下面的积分方程： 用和分别离散等号左右的Riemann Liouville积分，则有 合并同类项整理后得迭代公式 基于一阶矩形公式的显式方法不需要求解非线性方程，但精度较低。2.3 二阶卷积权逼近法问题 可以等价地转化为下面的积分方程： 用二阶卷积权逼近，得 整理得到迭代公式 基于二阶卷积权方法需要求解非线性方程，但精度较高。3 数值算例在本节中演示了一例基础的三项分数阶微分方程初值问题如何用左右矩形公式和二阶卷积权逼近法求解的过程，并根据不同的步长进行多次迭代比较，对比迭代过程中的阶数和近似解与精确解进行误差分析。3.1例题与解 其中令, ,精确解为3.2 基于一阶矩形公式显式方法的误差分析问题 可以等价地转化为下面的积分方程：用和分别离散等号左右的Riemann Liouville积分，则有 合并同类项整理后得到迭代公式利用matlab程序进行求解，相关计算文件为附录A。在条件下，不同步长得到的数值解与精确解的误差和收敛阶数如下表：表1 矩形公式解得数值解与精确解的误差和误差阶数表步长误差收敛阶0.2560822751087210.1244419684707411.0411342781031170.0617755557819011.0103651309143630.0309887634247360.9952908984552180.0156094544471050.989325071008480不同步长下得到的数值解与精确解折线图如下：图1 h=1/10矩形公式解得数值解与精确解图图2 h=1/20矩形公式解得数值解与精确解图图3 h=1/40矩形公式解得数值解与精确解图图4 h=1/80矩形公式解得数值解与精确解图图5 h=1/160矩形公式解得数值解与精确解图可以看出，步长越小，左右矩形公式解得的数值解与精确解的误差越小，收敛阶大致为1，说明基于左右矩形公式的显示方法是一阶的计算方法。3.3 基于二阶卷积权逼近法的误差分析问题 可以等价地转化为下面的积分方程： 用二阶卷积权逼近，得 整理得到迭代公式 利用matlab程序进行求解，相关计算文件为附录B。在条件下，不同步长得到的数值解与精确解的误差和收敛阶数如下表：表2 二阶卷积权逼近法解得数值解与精确解的误差和误差阶数表步长误差收敛阶0.0084041562911410.0021757507614671.9495896907263495.526388329040177e-041.9771044548450491.392163516228395e-041.9890082686771523.493377475027870e-051.994634226628701不同步长下得到的数值解与精确解折线图如下：图6 h=1/10由二阶卷积权逼近法得到的数值解与精确解图图7 h=1/20由二阶卷积权逼近法得到的数值解与精确解图图8 h=1/40由二阶卷积权逼近法得到的数值解与精确解图图9 h=1/80由二阶卷积权逼近法得到的数值解与精确解图图10 h=1/160由二阶卷积权逼近法得到的数值解与精确解图可以看出，步长越小，二阶卷积权逼近法解得的数值解与精确解的误差越小，收敛阶大致为2，说明二阶卷积权逼近法是二阶的计算方法。结 论本文主要研究了多项分数阶常微分方程的数值积分法，首先列出了本篇文章为了求解初值问题所用到的一些预备知识，即Riemann Liouville积分、左右矩形公式和二阶卷积权逼近的相关知识。随后计算出了多项分数阶常微分方程一般初值问题的三种不同的迭代公式，对比基于左右矩形公式的显式方法和二阶卷积权逼近法对于具体算例解的误差分析，可以看出基于左右矩形公式的显式方法不需要求解非线性方程，但精度低，为一阶的计算方法；二阶卷积权逼近法需要求解非线性方程，但精度高，为二阶的计算方法。所以通过本文论证，对于具体的多项分数阶常微分方程初值问题可以选用不同的数值解法来进行求解，线性方程可以选用基于左右矩形公式的显式方法求解，非线性方程可以选用二阶卷积权逼近法求解来确保计算结果的精确性。参考文献1 杨晨航.分数阶常微分方程和修正反常次扩散方程D.厦门大学博士学位论文,2008：11-122 王小东.RiemannLiouville 分数阶微积分及其性质证明D.太原理工大学硕士学位论文,2008：3-43 华东师范大学数学系.数学分析M.北京：高等教育出版社,2010：284-286.4 白占兵.分数阶微分方程边值问题理论及应用M.北京：科学出版社,2013：243-244.5 刘发旺，庄平辉，刘青霞.分数阶偏微分方程数值方法及其应用M.北京：科学出版社，2013：106-110.6 郭柏灵，蒲学科，黄凤辉.分数阶偏微分方程及其数值解M.北京：科学出版社,2013：29-35.7 吴迎，黄健飞，赵维加.分数阶常微分方程的数值解法C.青岛大学学报,2013,31(3):8 龚晟.简明微积分M.北京：高等教育出版社,2017： 114-118，192-195.9 林然,刘发旺.分数阶常微分方程初值问题的高阶近似J .厦门大学学报,2004,43(1): 25-30.10 孙志忠，高广花.分数阶微分方程的有限差分方法M.北京：科学出版社,2018:196-201.11 苏新卫, 刘兰冬.分数阶微分方程解的存在性J.山西大学学报(自然科学版),2007, 30(4): 12-1612 DELBSCO D, RODINO L. Existence and uniqueness for a nonlinear fractional differential equationJ.Journal of Mathematical Analysis and Applications,1996, 204(2): 609-62513 Lakshmikantham V, Vatsala A S. Basic theory of fractional differential equationsJ. Nonlinear Analysis,2008,69 (8): 2677-268214 Liu F, Anh, V and Turner, I and P Zhuang, Numerical simulation for solute transport in fractal porous mediaJ.ANZIAM J,2004, 45( E):461-473.15 Deng J Q, Ma L F. Existence and uniqueness of solutions of initial value problems for nonlinear fractional differential equationsJ. Applied Mathematics Letters, 2010, 23(6): 676-680致 谢本毕业论文是在我的指导老师黄健飞黄老师的帮助下完成的。今年伊始就多灾多难，由于疫情原因，和导师的交流受到了极大的阻碍，黄老师平日里工作繁忙，但还是每周抽空给我们小组成员在微信上开指导会，督促进度检阅成果，协助我们解决难题，推荐相关参考资料书籍。在我完成毕业论文的每个阶段，黄老师从选题到推荐资料，论文大纲，论文的中期检查，后期格式的调整等各个环节中都给予了我悉心的指导。这几个月以来，黄老师不仅在学业上给我以精心指导，同时还在思想给我以无微不至的关怀，在此谨向黄老师致以诚挚的谢意和崇高的敬意。同时，感谢我的小组成员赵胜伟、王达仁、沈飞翱，你们群策群力的头脑风暴极大地启发了本篇论文的写作。感谢在整个毕业论文完成期间和我密切交流的同学，和曾经在各个方面给予过我帮助的伙伴们，在此，我再一次真诚地向帮助过我的老师和同学表示感谢！另外，我必须感谢我的父母，和其他所有关心我爱护我的亲人。感谢长辈对我的教诲，感谢他们对我学业上的关心和支持，也因此我有足够的信心和勇气战胜前进路上的艰难险阻。 最后，衷心感谢论文答辩委员会的各位专家老师能够在百忙之中审阅我的论文，祝各位老师工作顺利。附 录AA1 based_on_juxinggongshi.mclc;clear all;format long%global alpha1 alpha2 alpha3 alpha1 = 0.9; alpha2 = 0.6; alpha3 = 0.3; % 这三个参数可以调整，但必须保证： 1 > alpha1 > alpha2 > alpha3 > 0%T = 2; % 区间长度, 可以调整%tau = 1/10; % 步长，可以调整，减小步长，可以减小误差%TT = 0:tau:T;N = length(TT)-1;%w_alpha1 = zeros(1,N+1); w_alpha1_2 = zeros(1,N+1); w_alpha1_3 = zeros(1,N+1); % 初始化系数for k=0:N w_alpha1(k+1) = (k+1)(alpha1)-k(alpha1)/gamma(alpha1+1); w_alpha1_2(k+1)= (k+1)(alpha1-alpha2)-k(alpha1-alpha2)/gamma(alpha1-alpha2+1); w_alpha1_3(k+1)= (k+1)(alpha1-alpha3)-k(alpha1-alpha3)/gamma(alpha1-alpha3+1);end%coe_0 = 1 + tau(alpha1-alpha2)*w_alpha1_2(0+1) + tau(alpha1-alpha3)*w_alpha1_3(0+1);U = zeros(1,N+1); % 初始化解向量for n = 1: N temp = tau(alpha1)*right_fun(TT(n-1+1),U(n-1+1); for k = 0 : (n-2) temp = temp - tau(alpha1-alpha2)*w_alpha1_2(n-k-1+1)*U(k+1+1) - tau(alpha1-alpha3)*w_alpha1_3(n-k-1+1)*U(k+1+1); temp = temp + tau(alpha1)*w_alpha1(n-k-1+1)*right_fun(TT(k+1),U(k+1); end U(n+1) = temp/coe_0; %end % 画图 True_sol = true_fun(TT);plot(TT,U,'-')hold onplot(TT,True_sol,'r*')legend('数值解','精确解','Location','North')% 计算误差err = max(abs(U-True_sol);errA2 right_fun.mfunction f_val = right_fun(t,u)global alpha1 alpha2 alpha3 temp1= gamma(3+alpha1)/gamma(3)*t2 + gamma(3+alpha1)/gamma(3+alpha1-alpha2)*t(2+alpha1-alpha2) + gamma(3+alpha1)/gamma(3+alpha1-alpha3)*t(2+alpha1-alpha3);temp2= (t(2+alpha1)2- u2;f_val = temp1 + temp2;endA3 true_fun.mfunction u = true_fun(t)global alpha1 u = t.(2+alpha1);end附 录BB1 based_on_juanjiquan.mclc;clear all;format long%global alpha1 alpha2 alpha3 alpha1 = 0.9; alpha2 = 0.6; alpha3 = 0.3; % 这三个参数可以调整，但必须保证： 1 > alpha1 > alpha2 > alpha3 > 0%T = 2; % 区间长度, 可以调整%tau = 1/10; % 步长，可以调整，减小步长，可以减小误差%TT = 0:tau:T;N = length(TT)-1;%BW1=zeros(1,N+1); BW1_2=zeros(1,N+1); BW1_3=zeros(1,N+1);BW1(1)=1;BW1_2(1)=1;BW1_3(1)=1;for k=1:N BW1(k+1)=(1-(1-(alpha1)/k)*BW1(k); BW1_2(k+1)=(1-(1-(alpha1-alpha2)/k)*BW1_2(k); BW1_3(k+1)=(1-(1-(alpha1-alpha3)/k)*BW1_3(k);end%w_alpha1 = zeros(1,N); w_alpha1_2 = zeros(1,N); w_alpha1_3 = zeros(1,N); % 初始化系数for k=0:N temp_w1=0;temp_w1_2=0;temp_w1_3=0; for i=0:k temp_w1=temp_w1+BW1(i+1)*BW1(k-i+1)/3(k-i); temp_w1_2 = temp_w1_2 + BW1_2(i+1)*BW1_2(k-i+1)/3(k-i); temp_w1_3 = temp_w1_3 + BW1_3(i+1)*BW1_3(k-i+1)/3(k-i); end w_alpha1(k+1)=(3/2)(-alpha1)*temp_w1; w_alpha1_2(k+1)=(3/2)(alpha2-alpha1)*temp_w1_2 ; w_alpha1_3(k+1)=(3/2)(alpha3-alpha1)*temp_w1_3;end%coe_0 = 1 + tau(alpha1-alpha2)*w_alpha1_2(0+1) + tau(alpha1-alpha3)*w_alpha1_3(0+1);U = zeros(1,N+1); % 初始化解向量for n = 1: N temp = 0; for k = 0 : (n-1) temp = temp - tau(alpha1-alpha2)*w_alpha1_2(n-k+1)*U(k+1) - tau(alpha1-alpha3)*w_alpha1_3(n-k+1)*U(k+1); temp = temp + tau(alpha1)*w_alpha1(n-k+1)*right_fun(TT(k+1),U(k+1); end temp0 = U(n); % % 求解非线性方程 while 1 temp1 = temp/coe_0 + tau(alpha1)*w_alpha1(0+1)*right_fun(TT(n+1),temp0)/coe_0; if abs(temp0-temp1) < 10(-7) U(n+1) = temp1; break; else temp0 = temp1; end end %end % 画图 True_sol = true_fun(TT);plot(TT,U,'-')hold onplot(TT,True_sol,'r*')legend('数值解','精确解','Location','North')%