一维双边空间常系数对流扩散方程的二阶数值算法

王亚力

摘要：本文致力于一维空间常系数对流扩散方程的数值解法，基于有限差分法，提出了一种二阶精度的加权Crank-Nicholson离散差分格式，并讨论了其收敛性，最后给出了数值算例。

关键词：有限差分法 加权Crank-Nicholson格式

一、引言

（一）分数阶微分简介

分数阶微分或导数即是阶数为分数阶的微分或者导数.我们平时熟知的微分或者导数都是整数阶的，而分数阶微分或导数就是对其阶数进行一个从整数到分数的推广而得到的.虽然公众对分数阶微分或导数比较陌生，但其实分数阶微分或导数也有很长的历史.分数阶微分或者导数最早是由LHospital提出的，1695年在其写给Leibniz的一封信中提到：在n阶微分中，n通常为整数，而当n=1/2时，微分dny/dxn有着什么样的意义，其表达式又为什么.但是这个问题一直没有得到解答，直到一个世纪之后，在1812年Laplace给出分数阶导数的定义，在其定义中主要用到了积分.而到了1819年，Lacroix利用Laplace定义给出了常见函数y=x的阶数为1/2的分数阶导数：，这个结果与现在常用Reimann-Liouville定义下的分数阶微分是一样的.在此之后，Fourier通過一个变换（现在称之为Fourier变换）给出了分数阶导数的定义.Abel也为分数阶微分或导数的发展做出了重要贡献，他在研究Abel积分方程时发现：一般情况下一个常数的分数阶微分不为零.到了19世纪30年代，Liouville在前人特别是Abel和Fourier的启发下，并结合Gamma函数成功的将分数阶导数应用到位势理论中，最终给出了自己对于分数阶导数的定义.

在这一百多年的历史中，分数阶微分或导数发展极其缓慢，究其原因，主要为：一，计算复杂，即随着微积分的阶数从整数推广到分数，其计算复杂程度有着一个爆炸性的增长;二，缺乏应用背景，由于分数阶微分或导数的发展与当时科学技术的发展所需的数学理论的脱节，导致分数阶微分或导数不受科学界主流所重视，在很长一段时间内只有少数的几位纯理论数学家对其进行研究和探索.

（二）分数阶微分方程的研究意义

随着科技的发展分数阶微分方程逐渐涌入主流学界的视野中，并在实际应用背景下越来越受到关注，很多实际问题在进行数学建模的过程中可以抽象为分数阶微分方程.分数阶微分方程同时也是微分方程理论中的重要组成部分，丰富了微分方程理论研究.分数阶微分方程主要包括分数阶常微分方程与分数阶偏微分方程，分数阶常微分方程相较于分数阶偏微分方程来说较为简单，分数阶偏微分方程较为复杂，基本不可求得解析解.分数阶偏微分方程分为：空间分数阶偏微分方程、时间分数阶偏微分方程和空间和时间导数均为分数阶的偏微分导数。

分数阶偏微分算子由于其定义而自带全局性、遗传性及记忆性等特性，所以分数阶偏微分方程比较适合解决比较复杂的实际问题，特别是模拟一些带有时序记忆特性的复杂变化过程.由于分数阶偏微分方程特性，其被广泛应用于各个学科领域的数学建模：例如记忆材料学、流体力学、物理学、生物学等等学科领域.除此之外，分数阶偏微分还在现代工程计算中有着广发的应用，由此可见分数阶偏微分方程在理论与应用中都举足轻重。

虽然分数阶微分方程的前景一片光明，但是现在对其研究还有着很多的不足，主要表现在：一数值算法不完善，主要数值算法为有限差分和有限元法;二虽然有了一些数值算法，但没有将其封装为较为成熟的数值计算软件，不能很好满足现在工程计算的需求;三一些具有挑战性的数值计算问题尚未解决，例如长时间段的数值计算以及大空间域的数值计算。

（三）分数阶偏微分方程数值解法研究现状

对于解分数阶偏微分方程来说，最理想的状况就是能求得其解析解，但是绝大多数是无法求得解析解的.虽然对于分数阶常微分方程已经有了多种求解析解的解法，如分离变量法、Fourier变换法、Laplace变化法、镜像法、Mellin变换法等，但是由于分数阶偏微分的特性，将这些方法往分数阶偏微分方程求解的推广并没有取得良好的效果，因此只能对绝大多数分数阶偏微分方程求数值解.

关于对分数阶偏微分方程数值解的研究，最近几十年取得了一些成就，最主要的数值解法包括：有限差分法、有限元法、有限体积法、谱方法、变分迭代法、同伦扰动法等等.在这些方法中有限差分法起步最早，发展的最完善，总体来说较为成熟.谱方法是最新的研究成果，效果最好，但是计算相当复杂，发展还不够完善.本文主要采用有限差分法来解决特定的分数阶偏微分方程.有限差分法主要包括L1，L2、L2C、经典Grunwald公式和移位的Grunwald公式等方法，而且大多数时候还会利用这几种方法之间的结合即加权平均来形成新的解法。

（四）本文研究内容

分数阶微分方程有很强的应用背景，是从一系列的物理应用场景中抽象提取出来的一类微分方程.对于解决与之前所有时间段都相关的问题有着非常重要的意义，特别是在物理、化学、材料力学等方面都有着非常重要的应用[1].

二、差分离散格式

令h为分数阶对流扩散方程的空间步长，，为时间步长，，，。令表示在网格剖分点的值.则表示源项在该网格点出的值。

对式（1a）中的进行二阶中心时间差商离散，并选取作为离散的网格点.而对于式（1a）中的采用相邻网格点和处的二阶中心空间差商，并进行加权平均离散，这样就可以获得这两项的二阶精度离散格式。

而对于式（1a）中的双边分数阶导数和则采用经典的Grunwald公式、向前移位的和向后移位的Grunwald公式进行离散，并取其加权平均而得到空间上具有二阶精度的离散格式，该离散格式如下：

三、数值算例

考虑常系数双边分数阶对流扩散方程（1），且其边初值条件如下：

其中，，则知此时该方程的精确解为：

为了验证该数值解法的精度为二阶，取，最大误差：，精度计算公式：，下表给出了此方程在本文提出的差分格式下的数值解的误差和精度.由表1可知，本文提出的差分离散格式在空间步长减小的情况下，误差也随之减小，且该数值解法具有二阶精度.经典或者移位的Grunwald公式只能获得一阶精度的数值解法，而利用加权Crank-Nicholson格式将两者结合起来则会获得二阶精度的数值解法，所以考虑将经典Grunwald公式和多个移位Grunwald公式进行加权平均获得更高精度的数值解法将是一个很好的研究方向.

参考文献：

[1]郭柏灵，蒲学科，黄凤辉.分数阶偏微分方程及其数值解[M].北京：科学出版2011：1-3，91-92.

[2]Hao Z P ， Sun Z Z ， Cao W R . A fourth-order approximation of fractionalderivatives with its applications[J]. Journal of Computational Physics， 2015， 281：787-805.

[3]Podlubny I.. Fractional differential equations[M]. San Diego：Academic Press，1999：1-4，88-89.

（作者单位：华南理工大学数学学院）