一类分数阶常微分方程的数值解
常微分方程的数值解
f ( x, y1 ) f ( x, y2 ) L y1 y2
（其中 L 为 Lipschitz 常数）则初值问题（ 1 ）存 在唯一的连续解。
求问题（1）的数值解，就是要寻找解函数在一 系列离散节点x1 < x2 <……< xn < xn+1 上的近似 值y1, y 2,…,yn 。 为了计算方便，可取 xn=x0+nh，(n=0,1,2,…)， h称为步长。
(1)，(2)式称为初值问题，(3)式称为边值问题。 在实际应用中还经常需要求解常微分方程组：
f1 ( x, y1 , y2 ) y1 ( x0 ) y10 y1 (4) f 2 ( x, y1 , y2 ) y2 ( x0 ) y20 y2
本章主要研究问题(1)的数值解法，对(2)~(4)只 作简单介绍。
得 yn1 yn hf ( xn1 , yn1 )
上式称后退的Euler方法，又称隐式Euler方法。 可用迭代法求解
二、梯形方法 由
y( xn1 ) y( xn )
xn1 xn
f ( x, y( x))dx
利用梯形求积公式： x h x f ( x, y( x))dx 2 f ( xn , y( xn )) f ( xn1 , y( xn1 ))
常微分方程的数言 简单的数值方法 Runge-Kutta方法 一阶常微分方程组和高阶方程
引言
在高等数学中我们见过以下常微分方程：
y f ( x, y, y) a x b y f ( x, y ) a x b (2) (1) (1) y ( x ) y , y ( x ) y 0 0 0 0 y ( x0 ) y0 y f ( x, y, y) a x b (3) y(a) y0 , y(b) yn
数值分析常微分方程数值解法
第8页/共105页
➢ 数值积分方法(Euler公式)
设将方程 y=f (x, y)的两端从 xn 到xn+1 求积分, 得
y( xn1) y( xn )
xn1 f ( x, y( x))dx :
xn
xn1 F ( x)dx
xn
用不同的数值积分方法近似上式右端积分, 可以得到计算 y(xn+1)的不同的差分格 式.
h2 2
y''( )
Rn1
:
y( xn1)
yn1
h2 2
y''( )
h2 2
y''( xn ) O(h3 ).
局部截断误差主项
19
第20页/共105页
➢ 向后Euler法的局部截断误差
向后Euler法的计算公式
yn1 yn hf ( xn1, yn1 ), n 0, 1, 2,
定义其局部截断误差为
y 计算 的n递1 推公式，此类计算格式统称为差分格式.
3
第4页/共105页
数值求解一阶常微分方程初值问题
y' f ( x, y), a x b,
y(a)
y0
难点: 如何离散 y ?
➢ 常见离散方法
差商近似导数 数值积分方法 Taylor展开方法
4
第5页/共105页
➢ 差商近似导数(Euler公式)
(0 x 1)
y(0) 1.
解 计算公式为
yn1
yn
hfn
yn
h( yn
2xn ), yn
y0 1.0
n 0, 1, 2,
取步长h=0.1, 计算结果见下表
13
一阶常微分方程组数值解法
k12 2
,
y2n
k22 2
,...,ymn
km2 2
)
ki4 hfi (tn h, y1n k13, y2n k23, ,...,ymn km3 )
(i 1,2,...m) (n 0,1,2,...)
一阶常微分方程组的R-K算法
(1) 输入：m, h,t, yi (i 1,2,...m),te; (2) u1 0.5h;u2 0.5h;u3 h;u4 h;u5 0.5h;
;
(5) 输出：t, yi (i 1,2,...m);
(6) if t tethen 停机 else 返回（3）。
一阶常微分方程组的R-K方法
1.2 刚性方程组
设常系数线形微分方程组
dy Ay (t)
dt
y 其中, A Rmm; ( y1, y2,...,ym )T； (t) (1(t),2 (t),...,m (t))T。
h 2
,
yn
k3
hf
(tn
h 2
,
yn
k4 hf (tn h, yn
k1 ) 2
k2 ) 2
k3)
（n 0,1,2,...)
这里向量
y1n
y hf1(tn , n ) k11
y k f y y n
y2n ...
，
1h
(tn ,
n
)
hf2
(tn ,
n
)
k12
n
3
)
k
4m
其分量形式是
yin1 ki1 hfi ki2 hfi
yin (tn (tn
1 6
, y1n
h 2
(ki1 2ki2 ki3 ki4
计算方法：一阶常微分方程的数值解
改进Euler公式
LOGO
Runge-Kutta法的基本思想（2）
以上两组公式都使用函数f ( x , y )在某些点上的 值的线性组合来计算y( xn1 )的近似值yn1。 Euler公式：每步计算一次f ( x , y )的值，为一阶方法。 改进Euler公式：需计算两次f ( x , y )的值，二阶方法。
( i 2, 3 , p )
于是可考虑用函数f ( x , y )在若干点上的函数值的 线性组合来构造近似公式，构造是要求近似公式在 ( xn , yn )处的Taylor展开式与解y( x )在xn处的Taylor展开 式的前面几项重合，从而使近似公式达到所需要的阶 数。即避免求偏导，又提高了方法的精度，此为RK方 法的基本思想。
'
xn1
xn
f ( x, y )dx (n 0,1,)
用yn 1 , yn 代替y( xn 1 ), y( xn ), 对右端积分采用 取左端点的矩形公式
则有
xn1 xn
f ( x , y )dx h f ( xn , yn )
yn1 yn h f ( xn , yn ) (n 0,1,2, ... )
LOGO
忽略高阶项,取近似值可得到Euler公式
yn1 yn h f ( xn , yn ) (n 0,1,2, ... )
3. 数值积分法区间 将方程 y' f ( x, y ) 在区间 [ xn , xn1 ] 上积分
xn1
xn
y dx
LOGO
yn (c1 c2 ) f ( xn , yn )h c2 [a2 f ( xn , yn ) b21 f ( xn , yn ) f ( xn , yn )]h O ( h )
常微分方程中的数值方法
常微分方程中的数值方法常微分方程是数学中的一个重要分支。
它主要研究的对象是随时间变化的函数。
在实际应用中，我们需要求解这些函数的解析解，但通常情况下，解析解并不容易得到，甚至是不可能得到。
因此，我们需要使用数值方法来求解这些函数的数值近似解。
在本文中，我们将介绍常微分方程中的数值方法。
一、欧拉法欧拉法是常微分方程数值解法中最基本的一种方法。
它是根据欧拉公式推导而来的。
具体地，我们可以将一阶常微分方程dy/dt=f(t,y)写成如下形式：y(t+h)=y(t)+hf(t,y(t))其中，h是步长，f(t,y)是t时刻y的导数。
欧拉法就是通过上面的公式进行逐步逼近，然后得到最终的数值解。
欧拉法的计算过程非常简单，但所得到的解可能会出现误差。
这是因为欧拉法忽略了f(t+h,y(t+h))和f(t,y(t))之间的变化。
因此，我们需要使用更为精确的数值方法来解决这个问题。
二、改进欧拉法为了解决欧拉法中的误差问题，我们可以使用改进欧拉法。
改进欧拉法又称作四阶龙格-库塔法。
它的基本思想是对欧拉法公式进行改进，以提高计算精度。
具体地，根据龙格-库塔公式，可将改进欧拉法表示为：y(t+h)=y(t)+1/6(k1+2k2+2k3+k4)其中，k1=h*f(t,y)k2=h*f(t+h/2,y+k1/2)k3=h*f(t+h/2,y+k2/2)k4=h*f(t+h,y+k3)改进欧拉法的计算过程比欧拉法要复杂些，但所得到的数值解比欧拉法更精确。
这种方法适用于一些特殊的问题，但在求解一些更为复杂的问题时，还需要使用其他的数值方法。
三、龙格-库塔法龙格-库塔法是求解常微分方程中数值解的常用方法之一。
它最常用的是四阶龙格-库塔法。
这种方法的基本思想是使用四个不同的斜率来计算数值解。
具体地，我们可以将四阶龙格-库塔法表示为：y(t+h)=y(t)+1/6(k1+2k2+2k3+k4)其中，k1=h*f(t,y)k2=h*f(t+h/2,y+k1/2)k3=h*f(t+h/2,y+k2/2)k4=h*f(t+h,y+k3)与改进欧拉法相比，龙格-库塔法的计算复杂度更高，但所得到的数值解更为精确。
一阶常微分方程的数值求解
作业 利用Euler方法和R-K方法求解一个 常微分初值问题，并比较数值结 果，计算数值解和解析解的误差。
在 [ xk , xk 1 ] 内多取几个点，将它们的导数加权平均代 替 f ( x, y( x)) ，设法构造出精度更高的计算公式。
常用的是经典的 四阶R-K方法
y0 y( x0 ), xk 1 xk h yk 1 yk h (L1 2 L2 2 L3 L4 )/6
若 f 在 D {a x b,| y | } 内连续，且满足 Lip 条件：
dy f ( x , y) , y( x0 ) y0 , x [a, b] dx
L 0, s.t.| f ( x, y1 ) f ( x, y2 ) || y1 y2 | , 则上述问题的连续可
Matlab函数数值求解
[T,Y] = solver(odefun,tspan,y0)
其中 y0 为初值条件，tspan为求解区间；Matlab在数值求解 时自动对求解区间进行分割，T (向量) 中返回的是分割点的 值(自变量)，Y (向量) 中返回的是解函数在这些分割点上的函 数值。
solver 为Matlab的ODE求解器（可以是 ode45、ode23、
其中
L1 L2 L3 L4
f ( xk , yk ) f ( xk h / 2, yk hL1 / 2) f ( xk h / 2, yk hL2 / 2) f ( xk h, yk hL3 )
例3：利用四阶R-K方法求解例1与例2，并与Euler方法的 数值解进行比较。
y( xk 1 ) y( xk ) y( xk 1 ) y( xk ) dy O ( h) dx x h h k
微分方程的解析与数值解法
微分方程的解析与数值解法微分方程既是数学分析的重要分支，也是许多学科领域的基础。
在实际问题的求解中，我们常常需要寻找微分方程的解析解或者数值解。
本文将围绕微分方程的解析和数值解法展开讨论。
一、微分方程的解析解解析解指的是通过代数计算得到的方程的解。
对于某些简单的微分方程，我们可以通过分离变量、变量代换等方法得到解析解。
下面以一阶线性常微分方程为例，讨论解的求解过程。
考虑一阶线性常微分方程形式如下：$$\frac{dy}{dx} + P(x)y = Q(x)$$其中，$P(x)$和$Q(x)$为已知函数。
我们可以通过以下步骤求解该微分方程：1. 将方程改写为标准形式：$\frac{dy}{dx} + P(x)y - Q(x) = 0$2. 求解齐次线性微分方程：$\frac{dy}{dx} + P(x)y = 0$。
记其解为$y_h$，即$y_h = Ce^{-\int P(x)dx}$，其中$C$为常数。
3. 利用常数变易法，假设原方程的解为$y = u(x)y_h$，其中$u(x)$为待定函数。
4. 将$y = u(x)y_h$代入原方程，得到关于$u(x)$的方程。
5. 求解$u(x)$的方程，得到$u(x)$的表达式。
6. 将$u(x)$代入$y = u(x)y_h$，得到原方程的解析解。
上述过程就是一阶线性常微分方程求解的一般步骤。
对于其他类型的微分方程，也有相应的解析解求解方法。
但并非所有微分方程都存在解析解。
二、微分方程的数值解法对于一些复杂的微分方程，无法找到解析解，此时我们需要借助数值方法求解。
常见的数值解法包括欧拉法、改进的欧拉法、四阶龙格-库塔法等。
1. 欧拉法欧拉法是一种较为简单的数值解法，其基本思想是通过离散化微分方程，将微分方程转化为差分方程。
具体步骤如下：将求解区间$[a, b]$等分成$n$个小段，步长为$h = \frac{b-a}{n}$。
利用微分方程的导数定义，将微分方程转化为差分方程，即$y_{i+1} = y_i + h \cdot f(x_i, y_i)$，其中$f(x, y)$为微分方程右端的函数。
常微分方程数值解法
y ( xi ).
f ( xi , y i )
yi+1 = yi + hf (xi，yi). 即为Euler公式。
若记
y ( xi 1 ) y ( xi ) h
y ( xi 1 ) f ( xi 1 , y ( xi 1 ))
yi+1 = yi + hf (xi+1，yi+1).
（希望）
yi h( c1 c2 ) y ' ( xi ) c2 a2 h y" ( xi ) O ( h )
2 3
14
希望：ei+1 = y(xi+1) – yi+1 = O(h3). 则应：
c1 c2 1 1 c2 a 2 2 b21 1 a2
改进Euler
1.0959
1.1841 1.2662 1.3434 1.4164
y(tn)
1.0954
1.1832 1.2649 1.3416
tn
0.6
0.7 0.8 0.9
Euler法
1.5090
1.5803 1.6498 1.7178 1.7848
改进Euler
1.4860
1.5525 1.6165 1.6782 1.7379
特例：a2 = 1 c1 = c2 = 1/2，b21 = 1，得2阶R-K公式：
h y y i ( K1 K 2 ) i 1 2 K1 f ( xi , y i ) K f ( x h, y hK ) i i 1 2
改进欧拉公式。
y(tn)
1.4832
1.5492 1.6125 1.6733 1.7321
分数阶微分方程的数值求解算法
分数阶微分方程的数值求解算法分数阶微分方程是一类在科学和工程领域中广泛应用的数学模型。
与传统的整数阶微分方程不同，分数阶微分方程包含了非整数阶导数，具有更广泛的物理意义和应用背景。
然而，由于分数阶导数的非局部特性，分数阶微分方程的数值求解相对困难。
为了克服这一问题，人们提出了多种数值求解算法。
一种常用的分数阶微分方程数值求解算法是基于离散化方法的。
这种方法将分数阶微分方程表示为一个离散的差分方程，并使用适当的差分格式进行数值求解。
其中，最为常见的是基于格点和差分点的有限差分法。
该方法通过将域离散化为一组格点和差分点，将分数阶导数转化为近似导数，进而得到离散化的差分方程。
然后，通过迭代求解离散化的差分方程，得到分数阶微分方程的数值解。
除了基于离散化方法的数值求解算法，还有一种常用的方法是基于变换方法的。
这种方法基于分数阶微分方程的变换理论，将分数阶微分方程转化为整数阶微分方程或者整数阶积分方程。
然后，使用传统的整数阶微分方程或积分方程的数值求解算法来求解得到的转化方程。
常见的变换方法包括拉普拉斯变换法、傅里叶变换法和小波变换法等。
这些方法在降低了求解难度的同时，可能会引入一定的误差，需要对误差进行适当的控制和修正。
此外，还有一些其他的数值求解算法被用于分数阶微分方程的求解，如基于插值法的算法、基于迭代法的算法和基于逼近法的算法等。
这些算法在特定情况下可能具有较好的数值性能和求解效果。
但是，无论使用何种算法，都需要注意数值稳定性和精度控制的问题，以确保得到可靠的数值解。
综上所述，分数阶微分方程的数值求解算法是一个重要而具有挑战性的问题。
各种不同的算法都有其适用范围和特点，在具体应用中需根据实际情况选择合适的算法。
随着对分数阶微分方程理论和方法的深入研究，相信会有更多高效、精确的数值求解算法被提出，并得到广泛应用。
常微分方程初值问题的数值解法
1
1、基本概念： 基本概念：
•一阶常微分方程初值问题是： 一阶常微分方程初值问题是： 一阶常微分方程初值问题是 y ′ =f(x,y) (1.1) y(x0)=y0 (1.2) 其中f是已知的xy平面上某个区域D上连续函数，式（1.1）是微分 方程，有无穷多解，式（1.2）是确定解的初始条件。 如一元函数y(x)对一切a≤x≤b 满足 （1） （x,y(x))∈D （2） y(x0)=y0 （3） y′存在，且y′(x)=f(x,y(x)) 则称y(x)是初值问题（1.1）、（1.2）在[a,b]上的解。 定义1.1 如果存在正常数L>0,使得对一切x∈[a,b]及y1 、y2 ∈[c,d] 有 |f(x,y1)-f(x,y2)|≤ L|y1- y2| (1.3) 则称f(x,y)满足对y的lipschitz(李普希滋）条件，其中L称为 lipschitz常数。
~ = y + hf ( x , y ) y k +1 k k k y k +1
2011-12-12
h y = y k + ( f ( x k , y k ) + f ( x k +1 , ~k +1 )) 2
12
Euler公式的几何意义 公式的几何意义
Y=y(x)
a
b
2011-12-12
13
例题：用Euler的数值解(取步长h=0.1计算到y3)：
2011-12-12 4
求初值问题，是给出它的解在某些节点数值的近似值，这称为 数值离散方法，若要求出在[a,b} 上的离散解，引入点列{xK} ， 其中 xk= xk-1 + hk-1 ，k=1，2，3….. xk称为节点， hk 称为步长。 通常，步长h不变，取为等距步长 hk =h= （b-a）/N ，N为等份 区间[a,b]分割数,节点变为等距节点，此时有 kh, xk= xk-1 + h= x0 + kh 记式（1.1），（1.2）的准确解y(x)在节点xk 之值为y(xk) ，而 记y(xk)的近似值为yk ，又记fk=f（xk ， yk ），
分数阶微分方程数值解的一种逼近方法讲解
分数阶微分方程数值解的一种逼近方法By:Pankaj Kumar, Om Prakash Agrawal摘要本文提出了一类分数阶微分方程（FDEs）的数值解方案.在这种方法中,FDEs 被Caputo型分数阶导数所表现. Caputo型分数阶导数的属性可以让一个分数阶微分方程减弱为一个Volterra型积分方程. 这样做了之后,许多研究Volterra 型积分方程的数值方法也可以应用于寻找FDEs的数值解. 本文总时间被划分为一组小区间,在两个连续区间中，用二次多项式逼近未知函数. 这些近似被替换成转化的Volterra型积分方程由此获得一组方程. 这些方程的解提供了FDE的解. 这种方法被应用于解决两种类型的FDEs,线性和非线性. 用这里给出的方法得到的解能与解析解和其他方法的数值解较好的吻合. 同时结果说明这种数值方法是稳定的.1.引言本文讨论分数阶微分方程的数值解. 分数阶导数和分数阶积分近年来收到了广泛的关注. 在许多实际应用中，分数阶导数和分数阶积分为考虑的系统提供了更加精确地模型. 比如，分数阶导数已经被成功地运用到模拟许多粘性材料的依赖频率的阻尼行为.1980年之前，Bagley 和Torvik提出了这个领域已经被研究的工作的一个回顾，并且说明了半阶导数模型可以非常好地描述阻尼材料的频率以来. 另一些学者说明了分数阶导数和分数阶积分在电化学过程，电解质极化，有色噪声，粘性材料和混沌领域的应用. Mainardi，Rossikhin和Shitikova 提出了分数阶导数和分数阶积分在一般固体力学，特定粘弹性阻尼模型中的应用的调查. Magin提出了分数阶微积分在生物工程的三个关键部分的回顾. 分数阶导数和分数阶积分在其他领域的应用以及相关的数学工具和技巧还可以在许多其他文献上找到.系统模型中分数阶导数的引进大多会导致分数阶微分方程的出现. 对某些特定的分数阶微分方程在通常系统条件下的解，已经有几种方法被找到. 这些方法包括，拉普拉斯变换，傅里叶变换，模态综合法和特征向量展开法，数值法以及基于Laguerre积分公式的方法. 然而，这些方法中大多数不能被应用到非线性分数阶微分方程. 更进一步的，正如Diethelm等人指出的，这些方法很多只能应用到特定类型的分数阶微分方程，并且人们并不知道他们能否被推广. 并且，在很多作者的研究成果中，并没有出现系统性的收敛性分析.最近，对于能被应用到线性和非线性分数阶微分方程的数值稳定数值逼近技巧，人们的兴趣愈发浓厚. 这些方法技巧大多利用了分数阶微分方程可以被减弱为Volterra型积分方程的特性. 因此，Volterra型积分方程的数值解法也可以应用到分数阶微分方程的解当中. Diethelm等人提出了分数阶微分方程数值解的一种PECE方法，其中P,C,E分别代表预测，校正和估计. 这样一来很多学者又推广了应用于常微分方程和分数阶微分方程的Adams–Bashforth–Moulton 型预测-校正格式. 这种方法的提出也是利用分数阶微分方程可以被转化为Volterra型积分方程的特点. 这些作者同时提出了误差分析和用Richardson外推法改善数值精度的延伸. Ford和Simpson提出了一种阶数大于1的分数阶微分方程的数值解法. 在该公式中，阶大于1的分数阶微分方程被减弱为阶小于1的分数阶微分方程，然后用相应的数值解法解由此导出的系统. 在所有这些方法当中，节点之间的未知函数用线性函数逼近. Kumar and Agrawal提出了阶数大于1的分数阶微分方程的数值解法. 这种方法要求就y(t)和它的导数在时间节点上连续.本文基于古典分数阶微分方程可以转化为Volterra型积分方程的特点也提出了一种数值方法来逼近分数阶微分方程的解. 特别地，我们用二次逼近函数来建立这种算法，结果说明这种方法可以被应用到寻求分数阶微分方程的数值解. 我们还通过两个例子，线性和非线性问题的解决，说明了这种方法的高效和准确，并且这种数值方法是稳定的.2.数值算法关于分数阶导数的定义已经出现有好几种，它们包括Riemann–Liouville, Grun-wald–Letnikov, Weyl, Caputo, Marchaud,和Riesz分数阶导数. 这里，我们规定使用Caputo导数.其中，Caputo导数的定义是, (n-1<α<n),(1)其中，α是导数的阶数，n是比α大的最小的整数.式(1)早在19世纪就在Liouville的论文中被提出，在Caputo的论文发表前一年它被Rabotnov所用. 然而，在文献中，被（1）式所定义的分数阶导数作为Caputo导数被广泛认知.在接下来的讨论中，我们考虑含有Caputo导数的初值问题:(2)在初始条件：, k=0,1,...,n-1,(3)下的解，其中，f是任意函数，是y的k阶导数，，k=0,1,…,n-1是指定初值条件. 假设这个函数关于参数和积分区间都是连续的，并且对于它的第二个参数满足Lipschitz条件.在纯数学中，Riemann–Liouville导数比Caputo导数应用更加广泛. 然而，这里考虑Caputo导数是因为以Riemann–Liouville导数为基础的分数阶微分方程要求y(t)在t=0点的导数和积分为0.一般来说，这些条件的物理意义不是已知的，并且在实际应用中，他们是不可用的. Lorenzo and Hartley讨论了寻找在更一般的情况满足下初始条件的正确格式的问题. 在Diethelm and Ford的文章中，方程(2)和(3)被证明可以等价描述为:,(4)其中g(t)为:. (5)为了解释以二次多项式为基础的数值方法，我们假设我要求的是由(2)式定义的分数阶微分方程从0到T的积分. 为了达到这个目的，我们把时间T等分成N 份，令h=T/N，作为时间区间的每一个部分. 时间在网格点上被表示为. 同时假设y(t)的数值逼近值被网格点所决定. 该方法的基本思想是在相邻的两个时间节点和上数值地获取函数y(t)的值，然后重复这个过程来接近所求积分直到取到终点T.为了便于接下来的讨论，我们规定如下记号:, ,这里的方法需要对方程(4)每一步求两次积分值. 这里有两种方法来达成目的.第一种,用一些近似函数逼近y(t)，然后用一种数值方法确定式(4)的积分值.这里需要在未知积分的情况下对和作初始的估计. 第二种, 都用近似函数来显式地逼近y(t)和f(t,y(t))以及确定式(4)的积分. 注意在这种情况下，和会作为参数出现在函数f当中. 本文利用的是第二种逼近方法.现在，我们给出算法的详细思路. 首先我们需要确定y(t)，，的值. 用二次插值函数可以在区间[0，]上逼近y(t)和f(t,y(t)):(6)以及(7)其中，是函数在第k个时间节点的值，，k=0,1,2是二次插值函数，其中下标（j,k）代表在第j+1，j=0,…,m步的第k个近似函数.我们首先确定y(t)在和处的值. 把(7)式带入(4)式，并积分得到:(8)其中,,(9)可以精确计算得到. 注意式(8)需要知道f在和的值(或者间接地说，y的值). 为了得到，在[，]上把f(t,y(t))近似为:，(10)其中，,是另一个二次插值函数. 函数,k=1,2,3由下式给出，函数由相似的办法定义.把(10)代到(4)中，积分得到:, (11)其中，, (12)可以(9)中一样被精确计算出来. 由(7)可以得到的值为:(13) 这里，我们充分利用了二次多项式的性质. 在非二次多项式的情况下，将会有不同的参数.把(13)代到(11)得到,+, (14)注意到(8)和(14)是关于两个未知量和的方程，可以用Newton–Raphson法，不动点迭代或者其他非线性方法求解. 这里，我们用Newton–Raphson法求解这些非线性方程. 这个方法需要对和作一个初始的估计. 当α大于1, 由t=0处的斜率可以得到关于和的更好的估计. 然而，在这里，对于α>1和α<1我们对把这些变量的初值估计为. 注意在每一次迭代式，时间步长取2h.现在我们假定在处，y的值是已知的，我们要求的是和处的值. 根据以上的逼近方法，和可以被表示为：+, (15)++(16)其中，,k=2m,2m+2,2m+2,,k=2m,2m+1/2,2m+1是和,用同样方法确定的系数. 注意(15),(16)的积分可以被数值确定因为y(t)在,处的值是已知的. 这些方程含有两个未知量和,而他们可以通过Newton–Raphson法得到. 本文中，我们把作为和的初值估计. 这样一来，方程(1)就可以在需要的区间上求积分.作为特殊情况，考虑如下非线性系统:.这种条件下，，式(16)和(15)减弱为:, (17) 其中=, (18)=, (19)=-, (20)=1+(21)=-, (22)=- . (23)(17)是一组线性联立方程，可以用线性方法求解.请注意以下两个补充说明. 第一，方程(1)只含有一个y(t). 当y(t)是一个向量函数时，算法同样成立.不过，所有的y(t)和f(t,y(t))必须换成相应的近似向量函数. 第二，算法需要保存所有算过的的y的值. 这是很多分数阶微分方程的典型特征. 这将会导致一些问题，特别是当y的维数和分数阶有限元系统一样大时. 这种情况下，系统有临近的短期记忆，y(t)在过去一段时间长度的值可以忽略不计,以此来改善对存储空间的需求和计算效率.3. 数值结果为了说明这种方法的效率,我们分别考虑一个线性和一个非线性的算例. 讨论这些例子是因为他们解的逼近格式是可靠的，并且可以用其他数值方法求解. 这样我们就可以把用这种方法得到的结果和解析解以及其他数值方法的结果相比较.３.１例１线性方程在第一个例子中，我们考虑如下给出的线性方程:，０＜α＜２，(24) 且. (25) 初始条件仅当α>1是成立. (24),(25)的解是：，(26) 其中,. (27) 是Mittag–Leffler函数的阶.图1.α=0.75时y(t)的比较（O:解析解，X:本文数值方法）表1.α=0.75时h在不同值下函数y(t)的误差对不同的α和h可以得到很多结果，这里给出其中一些. 在各种情况下，我们另T=6.4s.考虑这个区间是因为它接近α=2的系统的时间. 这里图(1)比较了α=0.75时的解析解和二次数值方法. 在这种情况下，我们令h=0.1s. 注意到这两个结果几乎完全重合. 为了强调收敛性，表(1)列出了当α=0.75,h分别等于0.4，0.2，0.1，0.05和0.025时的结果. 注意到随着步长的缩小，误差也如期望一样的缩小了. 在大部分节点误差比R=E(2h)/E(h)都非常接近3.37，这表明误差阶为1.75(或E(h)=O()).图2. α不同时y(t)的比较(O:α=0.25，X:α=0.5，+:α=0.75，Δ:α=0.95，*:α=1.)图3.α=1.5时y(t)的比较(O:解析解，X:本文数值方法)表2.α=1.5时h在不同值下函数y(t)的误差图(2)展示了h=0.025,α分别等于0.25，0.5，0.75，0.95和1时y(t)的数值结果.因为解析解和数值结果完全一致，因此图中没用画出解析解. 当α=1时，精确解为y(t)=. 注意到随着α越来越接近1，数值解越收敛到解析解y(t)=,即在极限情况下，分数阶微分方程的解接近整数阶导数的解.更进一步地，我们给出了α>1的一系列结果.α<1和α>1的结果是分离的，因为y(t)的斜率在α=1处有一个跃迁.图３比较了y(t)在α=1.5，h=0.4时的解析解和数值解. 两个结果又一次几乎完全重合.为了突出收敛性，表2给出了α=1.5，h分别等于0.4，0.2，0.1，0.05和0.025时的数值解. 正如之前观察到的一样，在这种情况下，随着步长的减小，误差也随之减小. 在这种情况下，误差比接近5.5，这表明误差阶为2.5(或E(h)=O()). 这样一来，通过观察α<1和α>1的收敛结果，可以知道误差的收敛阶为1+α(或E(h)=O()),即误差的阶不仅依赖于h，还依赖于导数的阶α.图4. α不同时y(t)的比较(O:α=1.25，X:α=1.5，+:α=1.75，Δ:α=1.95，*:α=2.)图5.α=0.25，0.75，1.25，1.75时y(t)的比较(O:解析解，X:本文数值方法)表3.本文数值方法和文献(35)中y(t)的误差的比较.图4展示了h=0.025,α分别等于1.25，1.5，1.75，1.95和2时y(t)的数值结果.又一次，因为解析解和数值结果完全一致，因此图中没用画出解析解. 当α=2时，精确解为y(t)=cos(t). 注意到随着α越来越接近2，数值解逐渐收敛到整数阶导数的解. 图2和图4展示的收敛结果十分重要，因为他们说明了在极限情况下，分数阶微分方程和他们的解逼近整数阶微分方程以及他们的解析解. 表3比较了t=1.0文献(35)的解的误差和用本文数值方法在α=0.5和1.25，h=0.1,0.05,0.025的解的误差. 注意到本文的方法得到了更低阶的误差. 这是因为，这里的方法是一种高阶方法. 当α和h取其他值时这种趋势也能显现出来.3.2 例2.非线性方程在第二个例子中，我们考虑一个如下定义的非线性方程:(28) 且. (29) 和之前一样，第二个初始条件仅适用于α>1. (28)(29)的精确解在文献(35)中已给出,(30)注意到当α<1, 方程的解在t=0处的斜率趋近于无穷. 因此，他可能导致在t=0附近出现一个巨大的数值误差.表4. α=0.75和1.5，h取不同值下y(t) 的误差.表5. 文献（35）中y(t)的误差和用本文数值方法得到的y(t)的误差的比较上面给出了一些在不同α和h下的数值结果. 图5表示了h=0.05,α分别等于0.25，0.75，1.25，1.75时解析解和数值解的结果. 由它可以得到(1).解析解和数值解基本重合，当α取其他值是，可以得到同样的结果. (2)尽管在t=0处斜率非常大，但是方法给出了非常精确的结果. （3）正如预期的那样, 在t=1处，对所有的α，y(t)的值收敛到0.25. 表4列出了当α=0.75和1.5，h=0.1,0.05和0.025的数值解的误差. 注意到误差随着步长的减小而减小. 同样的趋势在α取其他值时也能观察到. 在尝试过的α的值中，误差比R=E(2h)/E(h)表明没有特定的收敛阶. 然而，当α=0.75和1.5时，收敛的误差的平均值接近12和15.表5比较了文献(35)中在t=1.0处解的误差和用这种方法在α=0.25和1.25，h=0.05时的误差. 这里我们用的是文献(35)中用Richardson外推法得到的值. 观察得到，本文的方法又一次给出了更小的误差. 当α=0.25时，这种方法给出了比Richardson外推法小得多的误差. 这可能是因为，当α等于0.25时y(t)在t=0附近的斜率改变非常迅速,并且线性方法不能精确地获得结果. 需要指出的是，这种数值方法对于α和h取其他值时同样给出了更加精确的结果.4．结论本文给出了一种经典的分数阶微分方程的数值逼近方法. 这里的分数阶微分方程是依据Caputo型分数阶导数给出的, 这种导数的性质可以把分数阶微分方程减弱为Volterra型积分方程. 时间区间被分成一组网格，通过3个连续节点的二次插值多项式逼近未知的和已知的函数y(t)和f(t,y(t)). 把这些多项式带入Volterra型积分方程可以得到一组代数方程，这种数值方法的提出就是用来解这些方程以及获取y(t)的解. 通过一个线性一个非线性的例子的解决，说明了这种数值方法的作用. 用这种方法得到结果和解析解以及其他数值方法的结果是一致的. 还得到一个结论就是结果随着步长的减小而收敛. 在极限情况小，当α逼近整数值，数值方法会得到一个整数阶系统的解. 结果还表明这种方法是数值稳定的.。
分数阶微分方程的数值解法及其MATLAB实现
分数阶微分方程的数值解法及其MATLAB实现1.分数阶微分方程的定义和性质分数阶导数是一般导数的推广，可以用分数阶微分方程来描述分数阶导数的性质。
分数阶微分方程一般形式如下：D^αy(t)=f(t,y(t)),0<α≤1其中，α是分数阶指数，D^α是分数阶导数符号，y(t)是未知函数，f(t,y(t))是已知的函数。
2.分数阶微分方程的数值解法由于分数阶导数的定义较为复杂，常见的数值解法有两种：格点差分法和Laplace变换法。
2.1格点差分法格点差分法是将连续的分数阶导数问题离散化为离散的整数阶微分方程问题，然后利用传统的整数阶微分方程数值解法求解。
具体步骤如下：(1)将时间区间[0,T]平均分为N段，使得Δt=T/N。
(2)将分数阶导数D^α近似为整数阶导数D_m^β，其中m>α>m-1，β=m-α。
(3)将方程D^αy(t)=f(t,y(t))离散为差分方程D_m^βy(t_k)≈f(t_k,y(t_k))，其中t_k=kΔt，k=0,1,2,...,N。
(4)解差分方程，得到数值解y_k，k=0,1,2,...,N。
2.2 Laplace变换法Laplace变换法是将分数阶微分方程转化为Laplace变换形式的整数阶微分方程，然后利用Laplace变换的性质和经典的整数阶微分方程数值解法求解。
具体步骤如下：(1) 对分数阶微分方程D^αy(t) = f(t, y(t)) 进行Laplace变换，得到整数阶微分方程s^αY(s) - s^αy(0) + s^α-1Y(s) = F(s)，其中Y(s) 和 F(s) 分别为 y(t) 和 f(t,y(t)) 的Laplace变换。
(2)将Y(s)和F(s)用有限差分近似替代，并将方程离散化为差分方程，得到Y_k(s)和F_k(s)，其中k是离散时间步长。
(3) 解差分方程，得到 Y_k。
将 Y_k 用逆Laplace变换转换为 y_k。
常微分方程的数值解法分析解析
第八章 常微分方程的数值解法一．内容要点考虑一阶常微分方程初值问题:⎪⎩⎪⎨⎧==00)(),(y x y y x f dx dy微分方程的数值解：设微分方程的解y (x )的存在区间是[a,b ]，在[a,b ]内取一系列节点a= x 0< x 1<…< x n =b ,其中h k =x k+1-x k ；(一般采用等距节点，h=(b-a)/n 称为步长)。
在每个节点x k 求解函数y(x)的近似值:y k ≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。
用数值方法，求得f(x k )的近似值y k ，再用插值或拟合方法就求得y(x)的近似函数。
(一)常微分方程处置问题解得存在唯一性定理对于常微分方程初值问题:⎪⎩⎪⎨⎧==00)(),(y x y y x f dx dy如果：(1) 在B y y A x x 00≤-≤≤,的矩形内),(y x f 是一个二元连续函数。
(2) ),(y x f 对于y 满足利普希茨条件，即2121y y L y x f y x f -≤-),(),(则在C x x 0≤≤上方程⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy的解存在且唯一，这里C=min((A-x 0),x 0+B/L),L 是利普希茨常数。
定义：任何一个一步方法可以写为),,(h y x h y y k k k 1k Φ+=+，其中),,(h y x k k Φ称为算法的增量函数。
收敛性定理：若一步方法满足： (1)是p 解的.(2) 增量函数),,(h y x k k Φ对于y 满足利普希茨条件.(3) 初始值y 0是精确的。
则),()()(p h O x y kh y =-kh =x -x 0,也就是有0x y y lim k x x kh 0h 0=--=→)((一)、主要算法 1．局部截断误差局部截断误差：当y(x k )是精确解时，由y(x k )按照数值方法计算出来的1~+k y 的误差y (x k+1)- 1~+k y 称为局部截断误差。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
验证 了算 法的有 效性． 学 院 ｕ 学 关键 词 ： 分数 阶 常微分 方程 ； ａ ｕｏ分数 阶导数 ； Ｃｐ ｔ 降阶 法 ； 数值 解
Ｏ １５ 报 中图分 类号 ： ７
文 献标 识码 ： Ａ
文章编 号 ：６ ３—２ １ （ ０ １ ０ —０ ２ —０ １７ ６８ ２１ ）６ ０ ７ ４
０ 引 言
最近 几年 ， 数 阶微积 分在 许多 学科 和现 代工 程计 算 中得 以广泛关 注 和应 用 ． 分 在使 用分数 阶导 数 的模 型中， 大部 分 情况 下会 出现 一 系列 的分数 阶微 分方 程． 数 阶微分 方程 的数值 求解 方法 成为 了近年研 究 的 分 热点 课题 之一 ． ｌ ｒ Ｒ ｓ 在文 献［ ］ Ｍｉｅ 和 ｏ ｓ ｌ １ 中给 出了一种 将分 数 阶常系 数线性 常微 分方 程
Ｈ
摘 要 ： 对一 类分数 阶 常 系数 线性 常微 分 方程 ， 于降 阶 的思 想 ， 针 基 通过 转 换将 其 转化 为低 滨 阶 的分数 阶 方程 组 的形式 ， 构造 了一 种新 的数 值 解 法 ， 出 了具 体 的 计算 格 式 ， 给 并通 过数 值 算 例
州 Ｂ
ａ
（） １
收 稿 日期 ： ０ １—１ ２１ Ｏ—Ｏ ７
基 金 项 目 ： 家 自然 科 学 基 金 项 目（ ０ ７ ０ ８ ， 州 学 院青 年 人 才 创 新 工 程 科 研 基 金 项 目（ Ｚ ＹＱＮ Ｇ ０ ００ 国 １９ １ １ ）滨 ＢＸ Ｌ ２ １１ ） 作者简介 ： 王
１６ ｃ ｒ ． ３． ｏ ｎ
ａＤ￣（ ＋ｎ１ ” （ ＋…＋ｎＤｙ ） ０ ｏ（ 一０ 一÷，∈Ｎ ．Ｔ Ｄ ￡ Ｃｙ） 一 ｙ） １ Ｔ（ ＋ａ ￡ ， ｃ ） ｑ
转化 为分 数 阶常 微分 方 程组 的方 法 ， 在此 基础 之上 求 出了这类 方程 的解 析解 ． ｉｔｅｍ 和 Ｆ ｒ 并 Ｄ ｅｈ ｌ ｏｄ在 文献 ［ ］中将 这 种 降阶算 法应 用 到 Ｂ ｇｅ ２ ａ ｌｙ—Ｔ ｒ ｉ ｏ ｖｋ方程
—
Ｏ
其 中 ０＜
＜
＜ … ＜
≤ １ ｍ 一 ［ ，。 ， 口］ ＋ ＋ … ＋ 一 ．
利用 如下代 换把 式 （ ）改写 ， ３
Ｙ — Ｙ， Ｙ ： Ｄ ， ， ，０ ｌ — Ｏｙ ， ／ — Ｄ。 一 。 ｌ … ｎ＋ Ｙ 。 … ， Ｙ － － Ｄ Ｙ．
等 价 的 矩 阵 的 形 式 为
第 ２ 卷 ７
（ ２）
对 于方 程 （）不妨 设 ａ １， ＞
ａ ｏ ｏ ，
＞ … ＞ａ ＞ ａ ≥ ０ 将 方 程 （ ）改 写 成 ： 。 ， １
（ ３）
（ ＋∑ （， Ｊ ｉ（ ＋ｂ Ｄ ｄ （ ＋…＋６ Ｄ ｙ ） ＝ｆ ） ６Ｄ＋ ２ Ｊ２ １ ％ ） ， ＋， ） Ｊ ｊ ｊ Ｊ （） ，
证 明参 见文 献 Ｅ ］ ４． 由引理 １和引理 ２易 知上 述转 化是等 价 的．
Ｄ ［Ｄ ＿（） ［ ｆ（） ｅ ￡． ］一 Ｄｆ ＤＴ ］一 ＤＦ ｆ（） 厂 证 明 ： 见文 献［ ］ 参 ３．
引理 ２ 如 果 （）∈ Ｃ Ｅ ， ， ￡ ｏ 丁］ 这里 Ｔ＞ ０ 并且 是∈ Ｎ， ， 如果 ａ Ｎ ， ａ ｋ 则 Ｄ Ｙ ０ ０＜ ＜ ， （ ）一 。
Ａ Ｄ Ｙ＋ＢＤＰ ｙ＋Ｃ ｙ一 厂 ． （）
本 文 将这 种降 阶 的转化 方法 应用 到任 意 阶的常 系数 线性 常微 分方程
ａＣ ｙ ￡ ＋ ａ一 Ｄ ｙ ￡ ＋ … ＋ ａ Ｙ ￡ ＋ ａ Ｙ ． Ｄ￣ （） １ （） １ Ｄ （） ０ Ｔ （）一 厂 ￡ ， Ｄｏ （）
得到 这类 分数 阶常系 数线 性常 微分 方程 的数 值解 ， 以 Ｂ ｇｅ —Ｔｏｖｋ方 程 为例验 证算法 的有 效性 ． 并 ａ ｌｙ ｒｉ
分数 阶导数 的定 义有 很多 种 ， 中选 用 的是 Ｃ ｐ ｔ 数 阶导数 ． 文 ａ ｕｏ分
定 义 （ ａ ｕ ｏ分 数 阶 导 数 定 义 ［ ） 设 Ｙ ∈ Ｃ ， Ｃ ｐｔ 。 ］ ” 一 １ 口 ｎ ∈ Ｎ， ＞ ０则 称 ＜ ≤ ， ｔ
磊 （ ９ Ｏ ） 男 ， 东 滨 州 人 ， 师 ， 士 ， 要 从 事 微 分 方 程 的 数 值 解 法 研 究 ， － ａｌｋｌ １ｎ １８ 一 ， 山 讲 硕 主 Ｅ ｍ ｉ ｅ２ ｃ ＠ ： ｉ
２ ８
滨州 学 院学报
¨ （ ）一 ｃ ， ｋ一 ０ １ … ，口 ］ ０ ＾ ， ， ［ ．
Ｄ１ ０
，— —
０ Ｄ
ｉ
０ ０
］厂ｙ 一 Ｌｙ 『∑ － Ｐ＿ ｙ０ ］ ＿ｎ ｌ —
…
（） ４
其中
为 ｎ 的矩阵 ×
０
●
：
Ｄ＾ 一
： ● …来自０ 引理 １ Ｃ ｐｔ ａ ｕｏ分数 阶微分 算子 满足 交换 律 ， 且满 足叠加 关 系 并
㈤一 ＿ 』 而ｌ
１ 算 法 构 造
对 于分 数 阶线性 常 系数多 项 常微分 方程
ａ ｒ
为 ＿ 厂 ）在 Ｃ ｐ ｔ （ ａ ｕｏ意义 下 的 ａ阶导 数． 了方便起 见 ， 为 下文 中均 记为 Ｄ ．
口ｃ Ｄ Ｙ（）＋ ａ－ ￡ ｎ１Ｄ Ｙ（）４ … ＋ 口 ＇ － １ Ｙ（）＋ ａ ｏ ￡ ＤＴ ｏＤＴ Ｃ ｙ（）一 ，（） ，
第 ２ ７卷 第 ６期
Ｖｏ１ ２ Ｎｏ ． ７， ．６
２１ ０ １年 １ ２月
Ｄｅ ．， ０１ ｃ ２ １
【 分 方 程 与 动 力 系 统研 究】 微
一
类 分 数 阶 常 微 分 方 程 的数 值 解
王 磊
Ｔ Ｊ
Ｏ
（ 滨州 学 院 数 学与 信息科 学 系 ， 山东 滨州 ２ ６ ０ ） ５ ６ ３