一类分数阶常微分方程的数值解
常微分方程的数值解
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、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
验证 了算 法的有 效性． 学 院 ｕ 学 关键 词 ： 分数 阶 常微分 方程 ； ａ ｕｏ分数 阶导数 ； Ｃｐ ｔ 降阶 法 ； 数值 解
Ｏ １５ 报 中图分 类号 ： ７
文 献标 识码 ： Ａ
文章编 号 ：６ ３—２ １ （ ０ １ ０ —０ ２ —０ １７ ６８ ２１ ）６ ０ ７ ４
０ 引 言
最近 几年 ， 数 阶微积 分在 许多 学科 和现 代工 程计 算 中得 以广泛关 注 和应 用 ． 分 在使 用分数 阶导 数 的模 型中， 大部 分 情况 下会 出现 一 系列 的分数 阶微 分方 程． 数 阶微分 方程 的数值 求解 方法 成为 了近年研 究 的 分 热点 课题 之一 ． ｌ ｒ Ｒ ｓ 在文 献［ ］ Ｍｉｅ 和 ｏ ｓ ｌ １ 中给 出了一种 将分 数 阶常系 数线性 常微 分方 程
Ｈ
摘 要 ： 对一 类分数 阶 常 系数 线性 常微 分 方程 ， 于降 阶 的思 想 ， 针 基 通过 转 换将 其 转化 为低 滨 阶 的分数 阶 方程 组 的形式 ， 构造 了一 种新 的数 值 解 法 ， 出 了具 体 的 计算 格 式 ， 给 并通 过数 值 算 例
州 Ｂ
ａ
（） １
收 稿 日期 ： ０ １—１ ２１ Ｏ—Ｏ ７
基 金 项 目 ： 家 自然 科 学 基 金 项 目（ ０ ７ ０ ８ ， 州 学 院青 年 人 才 创 新 工 程 科 研 基 金 项 目（ Ｚ ＹＱＮ Ｇ ０ ００ 国 １９ １ １ ）滨 ＢＸ Ｌ ２ １１ ） 作者简介 ： 王
１６ ｃ ｒ ． ３． ｏ ｎ
ａＤ￣（ ＋ｎ１ ” （ ＋…＋ｎＤｙ ） ０ ｏ（ 一０ 一÷，∈Ｎ ．Ｔ Ｄ ￡ Ｃｙ） 一 ｙ） １ Ｔ（ ＋ａ ￡ ， ｃ ） ｑ
转化 为分 数 阶常 微分 方 程组 的方 法 ， 在此 基础 之上 求 出了这类 方程 的解 析解 ． ｉｔｅｍ 和 Ｆ ｒ 并 Ｄ ｅｈ ｌ ｏｄ在 文献 ［ ］中将 这 种 降阶算 法应 用 到 Ｂ ｇｅ ２ ａ ｌｙ—Ｔ ｒ ｉ ｏ ｖｋ方程
—
Ｏ
其 中 ０＜
＜
＜ … ＜
≤ １ ｍ 一 ［ ，。 ， 口］ ＋ ＋ … ＋ 一 ．
利用 如下代 换把 式 （ ）改写 ， ３
Ｙ — Ｙ， Ｙ ： Ｄ ， ， ，０ ｌ — Ｏｙ ， ／ — Ｄ。 一 。 ｌ … ｎ＋ Ｙ 。 … ， Ｙ － － Ｄ Ｙ．
等 价 的 矩 阵 的 形 式 为
第 ２ 卷 ７
（ ２）
对 于方 程 （）不妨 设 ａ １， ＞
ａ ｏ ｏ ，
＞ … ＞ａ ＞ ａ ≥ ０ 将 方 程 （ ）改 写 成 ： 。 ， １
（ ３）
（ ＋∑ （， Ｊ ｉ（ ＋ｂ Ｄ ｄ （ ＋…＋６ Ｄ ｙ ） ＝ｆ ） ６Ｄ＋ ２ Ｊ２ １ ％ ） ， ＋， ） Ｊ ｊ ｊ Ｊ （） ，
证 明参 见文 献 Ｅ ］ ４． 由引理 １和引理 ２易 知上 述转 化是等 价 的．
Ｄ ［Ｄ ＿（） ［ ｆ（） ｅ ￡． ］一 Ｄｆ ＤＴ ］一 ＤＦ ｆ（） 厂 证 明 ： 见文 献［ ］ 参 ３．
引理 ２ 如 果 （）∈ Ｃ Ｅ ， ， ￡ ｏ 丁］ 这里 Ｔ＞ ０ 并且 是∈ Ｎ， ， 如果 ａ Ｎ ， ａ ｋ 则 Ｄ Ｙ ０ ０＜ ＜ ， （ ）一 。
Ａ Ｄ Ｙ＋ＢＤＰ ｙ＋Ｃ ｙ一 厂 ． （）
本 文 将这 种降 阶 的转化 方法 应用 到任 意 阶的常 系数 线性 常微 分方程
ａＣ ｙ ￡ ＋ ａ一 Ｄ ｙ ￡ ＋ … ＋ ａ Ｙ ￡ ＋ ａ Ｙ ． Ｄ￣ （） １ （） １ Ｄ （） ０ Ｔ （）一 厂 ￡ ， Ｄｏ （）
得到 这类 分数 阶常系 数线 性常 微分 方程 的数 值解 ， 以 Ｂ ｇｅ —Ｔｏｖｋ方 程 为例验 证算法 的有 效性 ． 并 ａ ｌｙ ｒｉ
分数 阶导数 的定 义有 很多 种 ， 中选 用 的是 Ｃ ｐ ｔ 数 阶导数 ． 文 ａ ｕｏ分
定 义 （ ａ ｕ ｏ分 数 阶 导 数 定 义 ［ ） 设 Ｙ ∈ Ｃ ， Ｃ ｐｔ 。 ］ ” 一 １ 口 ｎ ∈ Ｎ， ＞ ０则 称 ＜ ≤ ， ｔ
磊 （ ９ Ｏ ） 男 ， 东 滨 州 人 ， 师 ， 士 ， 要 从 事 微 分 方 程 的 数 值 解 法 研 究 ， － ａｌｋｌ １ｎ １８ 一 ， 山 讲 硕 主 Ｅ ｍ ｉ ｅ２ ｃ ＠ ： ｉ
２ ８
滨州 学 院学报
¨ （ ）一 ｃ ， ｋ一 ０ １ … ，口 ］ ０ ＾ ， ， ［ ．
Ｄ１ ０
，— —
０ Ｄ
ｉ
０ ０
］厂ｙ 一 Ｌｙ 『∑ － Ｐ＿ ｙ０ ］ ＿ｎ ｌ —
…
（） ４
其中
为 ｎ 的矩阵 ×
０
●
：
Ｄ＾ 一
： ● …来自０ 引理 １ Ｃ ｐｔ ａ ｕｏ分数 阶微分 算子 满足 交换 律 ， 且满 足叠加 关 系 并
㈤一 ＿ 』 而ｌ
１ 算 法 构 造
对 于分 数 阶线性 常 系数多 项 常微分 方程
ａ ｒ
为 ＿ 厂 ）在 Ｃ ｐ ｔ （ ａ ｕｏ意义 下 的 ａ阶导 数． 了方便起 见 ， 为 下文 中均 记为 Ｄ ．
口ｃ Ｄ Ｙ（）＋ ａ－ ￡ ｎ１Ｄ Ｙ（）４ … ＋ 口 ＇ － １ Ｙ（）＋ ａ ｏ ￡ ＤＴ ｏＤＴ Ｃ ｙ（）一 ，（） ，
第 ２ ７卷 第 ６期
Ｖｏ１ ２ Ｎｏ ． ７， ．６
２１ ０ １年 １ ２月
Ｄｅ ．， ０１ ｃ ２ １
【 分 方 程 与 动 力 系 统研 究】 微
一
类 分 数 阶 常 微 分 方 程 的数 值 解
王 磊
Ｔ Ｊ
Ｏ
（ 滨州 学 院 数 学与 信息科 学 系 ， 山东 滨州 ２ ６ ０ ） ５ ６ ３