王 穎，樊孝仁

(太原工業(yè)學(xué)院 理學(xué)系，山西 太原 030008)

0 引言

分數(shù)階導(dǎo)數(shù)在描述不同物質(zhì)的記憶與遺傳性質(zhì)方面提供了有力的工具，在科學(xué)和工程的不同領(lǐng)域都用分數(shù)階方程(組)來描述動力系統(tǒng).分數(shù)階微積分在數(shù)學(xué)、物理、化學(xué)、控制等領(lǐng)域都有廣泛的應(yīng)用，已經(jīng)吸引了許多科學(xué)家和工程學(xué)家的關(guān)注和研究[1].同整數(shù)階微分方程一樣，許多分數(shù)階微分方程的解都是理論上存在，其解析解并不能明確給出.然而，除了模型，對于方程出現(xiàn)突然發(fā)散，收斂或分歧的情況下，對于抓住其關(guān)鍵點是至關(guān)重要的.因此，高精度數(shù)值解仍然是需要的[2].因此，許多分數(shù)階常微分方程的數(shù)值解法被提出，最常見的技巧是區(qū)域分解法[3-5]，變分法[6-7]，分數(shù)階微分變換法[8-12]，分數(shù)階差分法[13]及冪級數(shù)解法[14].同時，整數(shù)階微分方程的經(jīng)典方法也經(jīng)常用到.如 ：Laplace變換法[15].然而，常微分方程初值問題的數(shù)值解法的經(jīng)典方法之一Runge-Kutta方法卻很少被人們所提到.

本文將把Runge-Kutta方法應(yīng)用到分數(shù)階常微分方程的初值問題中，給出其算法格式.

1 定義

以下是該文可能用到的基本定義.

定義1[16]對t>0，函數(shù)x(t)的階數(shù)為α∈R+的分數(shù)階積分定義為

其中，Γ(·)是伽馬函數(shù).

定義2[16]對t>0，函數(shù)x(t)的階數(shù)為α∈[n-1，n)，n∈Z+的Riemann-Liouville分數(shù)階導(dǎo)數(shù)定義為

定義3[16]如果f(t，x(t))是連續(xù)函數(shù)，那么分數(shù)階微分方程

(1)

等價于以下的第二類非線性Volterra積分方程

(2)

2 算法

2.1 Runge-Kutta型方法簡介

Runge-Kutta型方法為單步方法，針對方程

(3)

給出其算法.一般而言，若xn已知，則

xn+1-xn=hΦ(tn，xn，h)，n=0，1，….

(4)

算法的核心在于構(gòu)造增量函數(shù)Φ. Taylor級數(shù)展開法是用函數(shù)f在同一點(tn，xn)的高階導(dǎo)數(shù)表示Φ(tn，xn，h),這不便于數(shù)值計算.Runge-Kutta型方法是用函數(shù)f在一些點上的值來表示Φ(tn，xn，h),使單步法的局部截斷誤差與Taylor級數(shù)展開法相等.我們針對方程(3)，在[0，h]上取m個點0=t1≤t2≤t3≤…≤tm≤h，若知道ki=f(ti，x(ti))，i=1，2，…，m,則用它們的一次組合去近似f，

(5)

問題是如何計算ki,(因x(ti)未知).一個直觀的想法是：設(shè)已知(t1，k1)=(t1，f(t1，x(t1))，由Eurler方法，x(t2)=x(t1)+(t2-t1)f(t1，x(t1))=x(t1)+(t2-t1)k1，于是

k2=f(t2，x(t1)+(t2-t1)k1).

(6)

同樣，利用Eurler方法，又可以從(t2，k2)算出，

k3=f(t3，x(t1)+(t2-t1)k1+(t3-t2)k2).

(7)

如此可以繼續(xù)進行下去，節(jié)點{ti}和系數(shù)ci可如此選擇，使得近似式(5)有盡可能的逼近階.

為便于推導(dǎo)，先引進若干記號.首先令ti=t1+aih，i=1，2，…，m.其中，ai與h無關(guān).再引進如下的三角系數(shù)陣：

其中，bij與h無關(guān)，

(8)

又ci≥0，

(9)

{ai}，{bij}和{ci}已給定，則Runge-Kutta法的計算過程如下：

un+1=un+hΦ(tn，un，h).

(10)

其中，

(11)

(12)

系數(shù){ai}，{bij}和{ci}按如下原則確定：將ki關(guān)于h展開，以之代到(11)，使得l次冪hl(l=0，1，…，p-1)的系數(shù)和

(13)

中的同次冪的系數(shù)相等.如此得到的算法(10)稱為m級p階Runge-Kutta算法.

常用的Runge-Kutta算法為如下的4級4階經(jīng)典的Runge-Kutta算法.

(14)

2.2 Runge-Kutta型方法的穩(wěn)定性與收斂性

定理1[17]設(shè)Φ(t，x，h)(t0≤t≤T，0≤h≤h0，x∈(-，+))關(guān)于u滿足Lipschitz條件(Lipschitz常數(shù)為L)，則單步法(10)穩(wěn)定.

定理2[17]設(shè)Φ(t，x，h)滿足定理1的條件，又單步法相容，則當(dāng)h→0時，它的數(shù)值解xn→x(t)，只要t0+nh→t，初值x0→x(t0).更設(shè)單步法(10)是p階單步法，則還有收斂估計|en·≤eLT[(1+Lh)|e0·+cThp]. 設(shè)f(t，x)(t0≤t≤T，x∈(-，+))上連續(xù)，且關(guān)于u滿足Lipschitz條件.由Φ(t，x(t)，h)的表達式(12)知Φ(t，x，0)=f(t，x)，即Runge-Kutta法相容，所以定理1和定理2對Runge-Kutta法恒成立.

2.3 分數(shù)階微分方程的Runge-Kutta法構(gòu)造

考慮如下的分數(shù)階微分方程：

(15)

的Runge-Kutta方法的構(gòu)造.其中，f(t，x(t))是[0，1]×R上的連續(xù)函數(shù).

由引理3，問題(15)可以寫成如下的積分形式：

(16)

步驟2：計算x1，

設(shè)xn-1已經(jīng)計算出，下一步xn.

這樣，我們就可以計算出方程(15)的數(shù)值解{xn}.

3 說明