分数阶微积分_高等应用数学问题的MATLAB求解chapter10

384 第 10章 数学问题的非传统解法 10.6 分数阶微积分学及其应用 第 3章详细介绍了微积分学的相关 内容 财务内部控制制度的内容财务内部控制制度的内容人员招聘与配置的内容项目成本控制的内容消防安全演练内容 ，后面还陆续介绍了基于MATLAB语言的微 积分问题的计算方法。我们都熟知，dny/dxn 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示 y 对 x的 n阶导数，但若 n = 1/2时 是什么含义呢？这是 300多年以前法国著名数学家 Guillaume Franc¸ois Antoine L’Hoˆpital 问过微积分学创造者之一 Gattfried Wilhelm Leibnitz的一个问题[47]。分数阶微积分理论 建立至今已经有 300年的历史了，但早期主要侧重于理论研究，近年来在很多领域都已 经开始应用分数阶微积分学理论，例如在自动控制领域出现了分数阶控制理论等新的分 支。本节将介绍分数阶微积分的定义及各种计算方法，并介绍分数阶线性及非线性微分 方程的求解方法。 10.6.1 分数阶微积分的定义 10.6.1.1 分数阶微积分的各种定义及性质 在分数阶微积分理论发展过程中，出现了很多种函数的分数阶微积分的定义，如 由整数阶微积分直接扩展而来的 Cauchy积分公式，Gru¨nwald-Letnikov分数阶微积分定 义、Riemann-Liouville分数阶微积分定义以及 Capotu定义等，本节将先介绍这些定义及 其等效关系，再给出分数阶微积分的各种性质： ①分数阶 Cauchy积分公式 从简单整数阶积分直接扩展出来 Dγf(t) = Γ(γ + 1) 2pj ∫ C f(τ) (τ − t)γ+1dτ (10-6-1) 其中 C为包围 f(t)单值与解析开区域的光滑曲线。 ② Gru¨nwald-Letnikov分数阶微积分定义 该定义为 aD α t f(t) = lim h→0 1 hα [(t−a)/h]∑ j=0 (−1)j ( α j ) f(t− jh) (10-6-2) 其中， ( α j ) 为二项式系数。 ③ Riemann-Liouville分数阶微积分公式 其分数阶积分的定义为 aD −α t f(t) = 1 Γ(α) ∫ t a (t− τ)α−1f(τ)dτ (10-6-3) 其中 0 < α < 1，且 a为初值，一般可以假设零初始条件，即令 a = 0，这时微分记号 可以简写成 D−αt f(t)。Riemann-Liouville定义为目前最常用的分数阶微积分定义。特别 地，D 左右侧的下标分别表示积分式的下界和上界[14]。 由这样的积分还可以定义出分数阶微分，假设分数阶 n − 1 < β ≤ n，则定义其分数 阶微分为 aD β t f(t) = dn dtn [ aD −(n−β) t f(t) ] = 1 Γ(n− β) dn dtn [∫ t a f(τ) (t− τ)β−n+1dτ ] (10-6-4) 10.6 分数阶微积分学及其应用 385 ④ Caputo分数阶微分定义 Caputo分数阶微分的定义为 0D α t y(t) = 1 Γ(1− γ) ∫ t 0 y(m+1)(τ) (t− τ)γ dτ (10-6-5) 其中，α = m+ γ, m为整数，0 < γ ≤ 1。类似地，Caputo分数阶积分定义为 0D γ t = 1 Γ(−γ) ∫ t 0 y(τ) (t− τ)1+γ dτ, γ < 0 (10-6-6) 可以 证明 住所证明下载场所使用证明下载诊断证明下载住所证明下载爱问住所证明下载爱问 [38]，对很广一类实际函数来说，Gru¨nwald-Letnikov分数阶微积分定义及 Riemann-Liouville 分数阶微积分定义是完全等效的，Caputo 定义和 Riemann-Liouville 定义的区别主要表现在对常数求导的定义上，前者对常数的求导是有界的 (为 0)，而后 者求导是无界的，Caputo定义更适用于分数阶微分方程初值问题的描述。本节主要研究 Riemann-Liouville分数阶微积分问题，Cauchy积分定义和其他几种定义有一些区别，这 将在后面内容中通过例子来演示。 分数阶微积分有如下各条性质[37]： ①解析函数 f(t)的分数阶导数 0Dαt f(t)对 t和 α都是解析的。 ② α = n为整数时，分数阶微分与整数阶微分的值完全一致，且 0D0t f(t) = f(t)。 ③分数阶微积分算子为线性的，亦即对任意常数 a, b，有 0D α t [af(t) + bg(t)] = a 0D α t f(t) + b 0D α t g(t) (10-6-7) ④分数阶微积分算子满足交换律，并满足如下的叠加关系 0D α t [ 0D β t f(t) ] = 0D β t [ 0D α t f(t) ] = 0D α+β t f(t) (10-6-8) 10.6.1.2 分数阶微积分的积分变换 函数的分数阶积分表达式的 Laplace变换为 L [ D−γt f(t) ] = s−γL [f(t)] (10-6-9) 函数分数阶微分的 Laplace变换为 L [ 0D α t f(t) ] = sαL [f(t)]− n−1∑ k=1 sk [ 0D α−k−1 t f(t) ] t=0 (10-6-10) 特别地，若函数 f(t)及其各阶导数的初值均为零，则L [0Dαt f(t)] = sαL [f(t)]。 函数 f(t)的分数阶微积分的 Fourier变换可以统一写成 F [ −∞Dαt f(t) ] = (jω)αF [f(t)] (10-6-11) 其中 α既可以为正数，表示分数阶微分，也可以为负值，表示积分。 386 第 10章 数学问题的非传统解法 10.6.2 分数阶微积分的计算 10.6.2.1 利用 Fourier级数计算周期函数的分数阶微积分 在介绍基于 Fourier级数的分数阶微积分计算之前，有必要先回顾一下 Fourier级数 及其求解方法。对周期为 2L的函数，可以利用第 3.2.2节介绍的 Fourier展开方法，将 t ∈ [−L,L]区间内信号 f(t)写成三角函数的级数形式 f(t) = a0 2 + ∞∑ n=1 ( an cos np L t+ bn sin np L t ) (10-6-12) 其中  an = 1 L ∫ L −L f(t) cos npt L dt, n = 0, 1, 2, · · · bn = 1 L ∫ L −L f(t) sin npt L dt, n = 1, 2, 3, · · · (10-6-13) 若 t ∈ (a, b)，则可以计算出 L = (b−a)/2，引入新变量 tˆ，使得 t = tˆ+L+a，则可以 将 f(tˆ)映射成 (−L,L)区间上的函数，可以对之进行 Fourier级数展开，再将 yˆ = t−L−a 转换成 t的函数即可。用该节编写的函数 fseries()可以得出 Fourier级数及 Fourier系 数 an, bn。 假设已知正弦、余弦函数，其整数阶微分表达式为 dk dtk [ sin at ] = ak sin ( at+ kp 2 ) , dk dtk [ cos at ] = ak cos ( at+ kp 2 ) (10-6-14) 由 Cauchy积分公式可以证明，对分数阶的微分来说，亦即 k为分数时，上述公式仍 然成立[43]。所以，可以考虑用整数阶 Fourier级数展开去逼近已知的 f(t)函数，利用式 (10-6-14)给出的式子，就可以将 f(t)函数的 γ 阶导数用下面的公式直接计算出来 aD γ t f(t)= a0 Γ(1−γ) t −γ + ∞∑ n=1 (np L )γ [ an cos (np L t+ γp 2 ) +bn sin (np L t+ γp 2 )] (10-6-15) 若 γ 为大于 1的分数，则可以将其分解成m + γ∗，其中m为正整数，0 < γ∗ < 1，这时 可以先对原 f(t)函数先求取m阶导数，再对结果求 γ∗阶微分。 若 γ 为负数，则立即计算出 |γ| 阶积分的值。根据上述的算法，可以编写出基于 Fourier级数展开的分数阶微分函数，其中 0 < γ < 1时建议采用此方法，否则用户可以 自己决定是否需要将 γ 分解成整数与小数之和的形式再进行微积分运算。 function F=fdiffint(A,B,t,gam,a,b) a0=A(1); A=A(2:end); n=length(B); L=(b-a)/2; if gam>=0, F=0; else, F=a0*t^gam/gam; end for i=1:n an=i*pi/L; bn=gam*pi/2; 10.6 分数阶微积分学及其应用 387 F=F+an^gam*(A(i)*cos(an*t+bn)+B(i)*sin(an*t+bn)); end if a+b, F=subs(F,t,t-L-a); end 【例 10-31】考虑例 3-19中给出的函数 y = t(t − p)(t − 2p), t ∈ (0, 2p)，试用 Fourier级数展开的 方法绘制出该函数的一些分数阶微分与积分函数曲线。 【求解】利用 Fourier级数展开方法，可以得出如图 10-43a所示的拟合效果，可见，对本例来 说，直接采用 Fourier级数，只用 12项就能得出较好的拟合效果。 >> syms t; f=t*(t-pi)*(t-2*pi); [A,B,F]=fseries(f,t,12,0,2*pi); ezplot(f,[0,2*pi]); hold on; ezplot(F,[0,2*pi]) 0 1 2 3 4 5 6 −10 −5 0 5 10 a) Fourier级数的拟合效果 0 0.5 1 1.5 2 2.5 3 −15 −10 −5 0 5 10 15 20 γ=0.75 γ=0.5 γ=0.25 γ=−0.75 b) 分数阶微积分函数曲线 图 10-43 基于 Fourier级数的方法求出的函数分数阶微积分 利用有限项的 Fourier级数，就可以得出各个阶次下的微积分函数曲线，如图 10-43b所示。 >> F1=fdiffint(A,B,t,0.25,0,2*pi); ezplot(F1,[0,pi]); hold on F1=fdiffint(A,B,t,0.5,0,2*pi); ezplot(F1,[0,pi]); F1=fdiffint(A,B,t,0.75,0,2*pi); ezplot(F1,[0,pi]); F1=fdiffint(A,B,t,-0.75,0,2*pi); ezplot(F1,[0,pi]); 用下面语句还可以绘制出分数阶微积分曲面，如图 10-44所示。 >> tt=0:0.1:2*pi; aa=-1:0.2:1; T=[]; for r=aa F1=fdiffint(A,B,t,r,0,2*pi); T=[T; double(subs(F1,t,tt))]; end surf(tt,aa,T) 【例 10-32】假设给定函数 f(t) = e−t sin(3t+ 1)，t ∈ (0, p)，试求出其分数阶导数。 【求解】对给定的区间进行 Fourier级数展开，用不同的项数进行拟合，则将得出如图 10-45a所 示的拟合效果，显然，这些拟合效果都不好，尤其表现在曲线的边界上。 >> syms t; f=exp(-t)*sin(3*t+1); ezplot(f,[0,pi]); hold on [A1,B1,F1]=fseries(f,t,10,0,pi); ezplot(F1,[0,pi]); 388 第 10章 数学问题的非传统解法 0 2 4 6 8 −1 −0.5 0 0.5 −10 0 10 20 γ t 图 10-44 不同阶次下的微积分曲面 [A2,B2,F2]=fseries(f,t,15,0,pi); ezplot(F2,[0,pi]); [A3,B3,F3]=fseries(f,t,20,0,pi); ezplot(F3,[0,pi]); 0 0.5 1 1.5 2 2.5 3 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 a) (0,p)区间 0 0.5 1 1.5 2 2.5 3 −0.4 −0.2 0 0.2 0.4 0.6 0.8 b) (−0.5,p+ 0.5)区间 图 10-45 不同指定区间的 Fourier级数展开拟合效果 为了解决这样的问题，可以适当放大函数的定义域，在区间 (−0.5, p + 0.5)上重新拟合原函 数，这样在 (0, p)区间段内原函数的拟合效果如图 10-45b所示，显然在这一区间内拟合效果明显 改善了。 >> syms t; f=exp(-t)*sin(3*t+1); ezplot(f,[0,pi]); hold on [A1,B1,F1]=fseries(f,t,10,-0.5,0.5+pi); ezplot(F1,[0,pi]); [A2,B2,F2]=fseries(f,t,20,-0.5,0.5+pi); ezplot(F2,[0,pi]); 利用这样得出的 Fourier级数展开，可以求出其分数阶导数，如图 10-46a所示，然而，对这 样给出的函数来说，得出的导数曲线不平滑，因为 Fourier级数不能平滑地拟合原函数。 >> F1=fdiffint(A2,B2,t,0.25,-0.5,pi+0.5); ezplot(F1,[0,pi]); hold on F1=fdiffint(A2,B2,t,0.5,-0.5,pi+0.5); ezplot(F1,[0,pi]) F1=fdiffint(A2,B2,t,0.75,-0.5,pi+0.5); ezplot(F1,[0,pi]) 对原始函数进行 60项 Fourier级数展开，则会发现函数曲线拟合效果有显著改善，图 10-46a 中还给出了分数阶微分函数曲线，低阶微分曲线很光滑，但高阶微分曲线仍然不光滑。 >> [A1,B1,F1]=fseries(f,t,60,-0.5,0.5+pi); 10.6 分数阶微积分学及其应用 389 0 0.5 1 1.5 2 2.5 3 −1 −0.5 0 0.5 1 a) 20项展开 0 0.5 1 1.5 2 2.5 3 −1 −0.5 0 0.5 1 b) 60项展开 图 10-46 函数的分数阶微分曲线 F1=fdiffint(A1,B1,t,0.25,-0.5,pi+0.5); ezplot(F1,[0,pi]); hold on F1=fdiffint(A1,B1,t,0.5,-0.5,pi+0.5); ezplot(F1,[0,pi]) F1=fdiffint(A1,B1,t,0.75,-0.5,pi+0.5); ezplot(F1,[0,pi]) 基于 Fourier级数的分数阶微积分计算方法完全取决于该函数是否能由 Fourier级数 精确近似，另外需要已知该函数。 10.6.2.2 用 Gru¨nwald-Letnikov定义求解分数阶微分 其实，求解分数阶微积分的最直接数值方法是利用 Gru¨nwald-Letnikov定义的方法 aD α t f(t) = lim h→0 1 hα [(t−a)/h]∑ j=0 (−1)j ( α j ) f(t− jh) ' 1 hα [(t−a)/h]∑ j=0 w (α) j f(t− jh) (10-6-16) 其中，w(α)j = (−1)j ( α j ) 为函数 (1 − z)α 的多项式系数，该系数还可以更简单地由下面 的递推公式直接求出 w (α) 0 = 1, w (α) j = ( 1− α+ 1 j ) w (α) j−1, j = 1, 2, · · · (10-6-17) 假设步长 h足够小，则可以用式 (10-6-16)直接求出函数数值微分的近似值，并可以 证明[38]，该公式的精度为 o(h)。所以，利用 Gru¨nwald-Letnikov定义可以立即编写出下 面的函数来求取给定函数的分数阶微分函数。 function dy=glfdiff(y,t,gam) h=t(2)-t(1); dy(1)=0; y=y(:); t=t(:); w=1; for j=2:length(t), w(j)=w(j-1)*(1-(gam+1)/(j-1)); end for i=2:length(t) dy(i)=w(1:i)*[y(i:-1:1)]/h^gam; end 390 第 10章 数学问题的非传统解法 其中 y, t分别为给定函数的采样值与时刻值构成的向量，要求 t为等间距向量，且 gam为 分数阶的阶次，这样得出的 dy 向量为函数的分数阶导数。 【例 10-33】考虑例 10-32中给出的函数 f(t) = e−t sin(3t+ 1)，t ∈ (0, p)，试求出其分数阶导数。 【求解】分别选择采样周期 T = 0.01和 0.001，则可以得出在这两个采样周期下计算出的 0.5阶导 数函数曲线，如图 10-47a所示，可见，在 t接近 0的区域外二者是很接近的，对本函数而言，选 择 T = 0.01可以足够精确地求出函数的分数阶微分。 0 0.5 1 1.5 2 2.5 3 −1 −0.5 0 0.5 1 1.5 2 2.5 a) 不同采样周期的结果比较 0 1 2 3 0 0.5 1 0 1 2 3 4 5 6 t γ b) 分数阶微分曲面 图 10-47 函数的分数阶微分 >> t=0:0.001:pi; y=exp(-t).*sin(3*t+1); dy=glfdiff(y,t,0.5); plot(t,dy); t=0:0.01:pi; y=exp(-t).*sin(3*t+1); dy=glfdiff(y,t,0.5); line(t,dy) 对不同的 γ 选值，还可以调用下面的语句绘制出分数阶导函数的三维图，如图 10-47b所示。 >> Z=[]; t=0:0.001:pi; y=exp(-t).*sin(3*t+1); for gam=0:0.1:1, Z=[Z; glfdiff(y,t,gam)]; end surf(0:0.01:pi,0:0.1:1,Z); axis([0,pi,0,1,-1.2,6]) 【例 10-34】试用不同定义求取函数为 f(t) = sin(3t+ 1)的 0.75阶微分，并比较得出的结果。 【求解】 由 Cauchy 定义，按式 (10-6-14) 给出的分数阶微分公式可以立即求出 0.75 阶微分为 0D 0.75 t f(t) = 3 0.75 sin ( 3t+1+ 0.75p 2 ) ，而用Gru¨nwald-Letnikov定义的微分可以由 glfdiff()函 数得出，这样由下面语句即可以绘制出由这两个定义得出的分数阶微分曲线，如图 10-48所示。 >> t=0:0.01:pi; y=sin(3*t+1); y1=3^0.75*sin(3*t+1+0.75*pi/2); y2=glfdiff(y,t,0.75); plot(t,y1,t,y2) 比较两种定义得出的结果，可见用 Cauchy积分公式定义的算法没有初始突变过程，在 t = 0 区域以外二者几乎是一致的。从得出的结果很难断定哪个定义是正确的，因为这完全取决于 t ≤ 0时 y 的值是什么，若 y = 0，则显然在 t = 0+ 时刻 y 的值从 0突变到 sin 1，所以这时分 数阶微分值应该为∞，突变的影响亦将持续一段时间，故 Gru¨nwald-Letnikov定义能正确反映 该函数的分数阶微分变化，而不适合用 Cauchy 积分公式；若在 t ≤ 0 时原函数仍然满足函数 y(t) = sin(3t+ 1)，则在 t = 0+时函数没有突变，则更适合使用 Cauchy积分公式。 10.6 分数阶微积分学及其应用 391 0 0.5 1 1.5 2 2.5 3 −3 −2 −1 0 1 2 3 4 5 6 图 10-48 不同定义下的分数阶微分曲线 10.6.2.3 分数阶微积分的 Fourier变换算法 由前面介绍的 Fourier变换理论可以将给定函数 f(t)的 Fourier变换及反变换为 F (ω) = ∫ ∞ −∞ f(t)e−jωtdt, f(t) = 1 2p ∫ ∞ −∞ F (ω)ejωtdω (10-6-18) 双边取分数阶微分，则可以得出 Dγf(t) = Dγ [ 1 2p ∫ ∞ −∞ F (ω)ejωtdω ] = 1 2p ∫ ∞ −∞ F (ω)Dγ [ ejωt ] dω = 1 2p ∫ ∞ −∞ (jω)γF (ω)ejωtdω (10-6-19) 亦即该函数的 γ 阶微分实际上就是 (jω)γF (ω)函数的 Fourier反变换。这些公式还有离散 版本，对给定序列 f(k)，可以得出该序列的离散 Fourier变换 F (ω) = ∞∑ k=−∞ f(k)e−jkω (10-6-20) 这样，函数分数阶导数可以表示成 Dγf(k) = 1 2p ∞∑ k=−∞ (jω)γF (ω)ejkω (10-6-21) 10.6.2.4 分数阶微积分的滤波算法 前面介绍的各种分数阶微分运算的前题是被微分函数 f(t)为已知函数，在实际应用 中该信号经常是无法预先知道的，所以应该采用其他形式来求取分数阶微分，例如通过 构造滤波器的方式来对信号进行数值微分处理。 信号的滤波器可以有连续和离散两种形式，分别用来拟合 Laplace 变换算子 sγ 和 Fourier变换算子 (jω)γ，从效果上看，函数的分数阶数值微分相当于原来信号需要通过这 样的滤波器得出的输出信号。 392 第 10章 数学问题的非传统解法 文献 [37] 中列出了多种连续滤波器的实现算法，这里只介绍其中的 Oustaloup 算 法[35]，假设选定的拟合频率段为 (ωb, ωh)，则可以构造出连续滤波器的传递函数模型 Gf (s) = K N∏ k=−N s+ ω′k s+ ωk (10-6-22) 其中滤波器零极点、增益可以由下式直接求出 ω′k=ωb ( ωh ωb )k+N+12 (1−γ) 2N+1 , ωk=ωb ( ωh ωb )k+N+12 (1+γ) 2N+1 , K= ( ωh ωb )−γ 2 N∏ k=−N ωk ω′k (10-6-23) 根据上述的算法，可以直接编写出如下的函数， 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 连续滤波器，这样若 y(t)信号 通过滤波器进行过滤，则可以认为输出的信号是 Dγt y(t)的近似。 function G=ousta_fod(r,N,w_L,w_H) mu=w_H/w_L; k=-N:N; w_kp=(mu).^((k+N+0.5-0.5*r)/(2*N+1) )*w_L; w_k=(mu).^((k+N+0.5+0.5*r)/(2*N+1) )*w_L; K=(mu)^(-r/2)*prod(w_k)/prod(w_kp); G=tf(zpk(-w_kp’,-w_k’,K)); 其中，r为分数阶的阶次，2N + 1为滤波器的阶次，w L和 w H分别为用户选定的拟合频 率下限 ωb和上限 ωh，一般在该区域内能较好地拟合分数阶微分算子，而其外的区域将和 微分算子相差很多。 【例 10-35】假设 ωb = 0.01, ωh = 100 rad/sec，试设计出连续滤波器，对 f(t) = e−t sin(3t+ 1)信 号计算 0.5阶微分。 【求解】可以用下面的语句设计出 >> G=ousta_fod(0.5,2,0.01,100), bode(G) Transfer function: 10 s^5 + 298.5 s^4 + 1218 s^3 + 768.5 s^2 + 74.97 s + 1 ------------------------------------------------------- s^5 + 74.97 s^4 + 768.5 s^3 + 1218 s^2 + 298.5 s + 10 >> t=0:0.001:pi; y=exp(-t).*sin(3*t+1); y1=lsim(G,y,t); y2=glfdiff(y,t,0.5); plot(t,y1,t,y2) 当然，用该算法还可以在更大的频率范围内拟合分数阶微分函数，这时需要适当增大拟合 的阶次，下面给出在 (10−4, 104)频段内的拟合效果，如图 10-50所示。可见，对这样大的频率范 围，不再适合 N = 2的近似，而应该采用更大的 N 值，如 N = 4。 >> G=ousta_fod(0.5,2,1e-4,1e4); G1=ousta_fod(0.5,3,1e-4,1e4); G2=ousta_fod(0.5,4,1e-4,1e4); G3=ousta_fod(0.5,5,1e-4,1e4); bode(G,’-’,G1,’--’,G2,’:’,G3,’-.’) 10.6 分数阶微积分学及其应用 393 −20 −10 0 10 20 M ag ni tu de (d B) 10−3 10−2 10−1 100 101 102 103 0 30 60 Ph as e (de g) Bode Diagram a) 不同采样周期的结果比较 0.5 1 1.5 2 2.5 3 0 1 2 3 4 b) 分数阶导数曲面 图 10-49 函数的分数阶导数 −50 0 50 M ag ni tu de (d B) 10−5 100 105 0 30 60 Ph as e (de g) Bode Diagram Frequency (rad/sec) 图 10-50 不同阶次下的滤波器近似效果 对离散系统来说，由于纯粹的分数阶微分器是不能直接物理实现的，所以需要用 FIR 滤波器 [43]或 IIR滤波器[4]的形式对其近似。利用滤波器设计工具箱中的 filt()函数，由 Ivo Petra´sˇ博士开发的分数阶微积分器的 FIR滤波器设计函数①，其核心部分 清单 安全隐患排查清单下载最新工程量清单计量规则下载程序清单下载家私清单下载送货清单下载 如下 function H=dfod2(n,T,r) if r>0 bc=cumprod([1,1-((r+1)./[1:n])]); H=filt(bc,[T^r],T); elseif r<0 bc=cumprod([1,1-((-r+1)./[1:n])]); H=filt([T^(-r)],bc,T); end 其中 n为期望的滤波器阶次，T 为滤波器的采样周期，r为所需的导数阶次，用该函数可 以直接设计出滤波器 H。若某信号经过滤波器 H，则可以得出该信号的 r阶导数，经过 滤波器得出的输出信号可以由 filter()函数计算出来。 【例 10-36】试用离散 FIR滤波器对例 10-35中给出的信号求取分数阶微分。 ①可以从MathWorks免费下载，见 http://www.matlabcentral.com下的信号处理滤波器设计栏目。 394 第 10章 数学问题的非传统解法 【求解】用下面的语句可以直接设计出滤波器，并绘制出 Bode图，如图 10-51a所示，由得出的 频率响应看，n = 100明显优于 n = 50时的滤波器，但显然远不如前面介绍的用 Oustaloup算法 设计的连续滤波器，在这些滤波器下还可以得出给定信号的分数阶微分，如图 10-51b所示，然而 计算精度不甚理想。 >> t=0:0.001:pi; y=exp(-t).*sin(3*t+1); y3=glfdiff(y,t,0.5); G1=dfod2(50,0.001,0.5); G2=dfod2(100,0.001,0.5); bode(G1,G2) y1=filter(G1.num{1},G1.den{1},y); y2=filter(G2.num{1},G2.den{1},y); figure; plot(t,y1,t,y2,t,y3) 0 10 20 30 M ag ni tu de (d B) 100 101 102 103 0 30 60 Ph as e (de g) Bode Diagram n=50 n=100 n=100 n=50 a) 滤波器的频域响应曲线 0.5 1 1.5 2 2.5 3 −1 0 1 2 3 4 n=50 n=100 b) 分数阶单数的滤波近似 图 10-51 不同 (n, a)组合下的数值微分器效果比较 FIR滤波器的阶次较高，所以可以考虑采用 IIR滤波器来生成分数阶微积分信号。从 连续系统角度看，分数阶微积分可以用 Laplace变换写成 s±γ，对其进行离散化，则可以 引入变换函数 s = w (z−1)来近似分数阶微分或积分运算。 文献 [4]给出了一种选择变换函数的新方法，该方法兼顾了 Simpson积分算法和梯形 积分算法，引入了加权系数 a，使得滤波器可以表示为H(z)=aH S (z)+(1−a)H T (s)，a ∈ [0, 1]。这样可以推导出 IIR滤波器 G ( z−1 ) = k±γ0 [ 1− z−2 (1 + bz−1)2 ]±γ (10-6-24) 其中，b = 3 + a− 2 √ 3a 3− a ，k0 = 6b T (3− a)，且 γ ∈ (0, 1)，T 为离散化的采样周期。对高 阶微积分来说，可以先进行整数阶微积分，对其结果再进行分数阶微积分。文献 [4]还给 出了基于连分式的滤波器设计方法和程序，对其修改后如下 function H=iir_frac(r0,a,n,T) syms x r b; aas=((1-x^2)/(1+b*x)^2)^r; n2=2*n; b0=(3+a-2*sqrt(3*a))/(3-a); aas=char(subs(aas,{r,b},{r0,b0})); maple(’with(numtheory)’); maple([’cfe:=cfrac(’ aas ’,x,n2)’]); 10.6 分数阶微积分学及其应用 395 p=collect(maple(’nthnumer’,’cfe’,n2),x); n=sym2poly(p); q=collect(maple(’nthdenom’,’cfe’,n2),x); d=sym2poly(q); H=tf(n(end:-1:1)*(6*b0/T/(3-a))^r0,d(end:-1:1),’Variable’,’z^-1’,’Ts’,T); 其中 r0为阶次，a为加权系数，n为滤波器阶次，T 为采样周期，这样可以由该函数直接 设计出离散的滤波器模型 H。下面将通过例子演示 IIR滤波器的设计及分数阶微积分近 似解法。 【例 10-37】试用不同的 n, a组合构造出 s0.5 阶滤波器，并比较其频域响应效果。 【求解】 用前面给出的程序可以立即设计出如下的滤波器，并绘制出滤波器的 Bode 图，如图 10-52a所示，从得出的频率响应曲线看，若选择 a = 0.5则相位曲线在给定的区域内不为恒值，所 以拟合效果将不理想，但高频处赋值将有下降趋势，对噪声将有抑制作用。 >> H1=iir_frac(0.5,0,3,0.01); H2=iir_frac(0.5,0.5,3,0.01); H3=iir_frac(0.5,0,4,0.01); H4=iir_frac(0.5,0,5,0.01); bode(H1,’-’,H2,’--’,H3,’:’,H4,’-.’) 用 Hn,a ( z−1 )来表示滤波器，则可以得出设计出的各类滤波器，例如 H3,0 ( z−1 ) = 113.1− 56.57z−1 − 56.57z−2 + 14.14z−3 8 + 4z−1 − 4z−2 − z−3 0 10 20 30 M ag ni tu de (d B) 10−1 100 101 102 −45 0 45 Ph as e (de g) a) 不同 (n, a) 0.5 1 1.5 2 2.5 3 0 1 2 3 4 5 b) 滤波效果比较 图 10-52 不同 (n, a)组合下的数值微分器效果比较 用下面的语句还可以绘制出各个滤波器的滤波效果，如图 10-52b所示，然而从得出的滤波效 果看，阶次太低则近似效果不佳，所以在设计中应该适当增加拟合阶次。 >> t=0:0.01:pi; y=exp(-t).*sin(3*t+1); y0=glfdiff(y,t,0.5); y1=filter(H1.num{1},H1.den{1},y); y2=filter(H2.num{1},H2.den{1},y); y3=filter(H3.num{1},H3.den{1},y); y4=filter(H4.num{1},H4.den{1},y); plot(t,y0,t,y1,’:’,t,y2,’--’,t,y3,’:’) 用Maple语言提供的连分式近似函数 cfrac()有局限性，当选择的 n ≥ 6不能得出 有意义的近似，所以可以考虑 Pade´近似方法，用其直接拟合理想滤波器的 Taylor幂级数 展开式，这样，改写原来的离散化函数，则可以写出如下的滤波器设计函数，该函数可 以用于高阶滤波器的设计。 396 第 10章 数学问题的非传统解法 function H=iir_pade(r0,a,n,T) syms x r b; aas=((1-x^2)/(1+b*x)^2)^r; b0=(3+a-2*sqrt(3*a))/(3-a); aas=subs(aas,{r,b},{r0,b0}); c=taylor(aas,x,2*n); c=sym2poly(c); c=c(end:-1:1); [n,d]=padefcn(c,n-1,n); H=tf(n(end:-1:1)*(6*b0/T/(3-a))^r0,d(end:-1:1),’Variable’,’z^-1’,’Ts’,T); 【例 10-38】仍考虑前面的例子，试提高滤波器的阶次，观察近似效果。 【求解】采用Maple语言提供的连分式方法不能处理这样的高阶问题，所以应该采用上面给出的 Pade´算法，可以设计出相应的滤波器，并绘制出各个滤波器的 Bode图，如图 10-53a所示。从实 际滤波器的增益看，和真值相差较远，虽然比低价滤波器有所改善。 >> H1=iir_pade(0.5,0,10,0.01); H2=iir_pade(0.5,0.75,10,0.01); H2=minreal(H2); bode(H1,’:’,H2,’--’) −20 −10 0 10 20 30 M ag ni tu de (d B) 10−2 10−1 100 101 102 −45 0 45 90 Ph as e (de g) a) 不同 a值的 Bode图 0.5 1 1.5 2 2.5 3 0 1 2 3 4 b) 实际滤波效果 图 10-53 n = 10时滤波器实现及效果比较 用下面的语句还可以绘制出滤波效果，如图 10-53b所示，其中 a = 0的滤波器产生的分数阶 微分信号被噪声淹没，所以这里只给出 n = 10, a = 0.5时滤波器的滤波效果，可见，采用高阶滤 波器时精度将明显改善，较好地逼近了实际的分数阶信号，由此可以看出 IIR滤波器比 FIR滤波 器表现出来的巨大优势。 >> t=0:0.01:pi; y=exp(-t).*sin(3*t+1); y0=glfdiff(y,t,0.5); y2=filter(H2.num{1},H2.den{1},y); plot(t,y0,t,y2,’:’) 10.6.3 分数阶微分方程的求解方法 10.6.3.1 分数阶线性微分方程的解法 分数阶线性微分方程的一般形式为[38] anD βn t y(t) + an−1D βn−1 t y(t) + · · ·+ a1Dβ1t y(t) + a0Dβ0t y(t) = u(t) (10-6-25) 其中 u(t)可以由某函数及其分数阶微分构成。为后面叙述方便，不妨对微分方程式作假 设 βn > βn−1 > · · · > β1 > β0 > 0。若需要研究的微分方程出现下面两种特殊情况，则需 要先对其变换，再进行求解。 10.6 分数阶微积分学及其应用 397 ①若研究的微分方程各个阶次不满足上述大小关系，则可以先进行排序。 ②若涉及到负的 βi 值，则原来的方程为分数阶微积分方程，需要选择引入新的变量 z(t) = Dβ0t y(t)，这样就可以将原来的微积分方程变换成关于 z(t)的微分方程了。 假设函数 y(t)具有零初始条件，则可以对该方程进行 Laplace变换，得出 G(s) = Y (s) U(s) = 1 ansβn + an−1sβn−1 + · · ·+ a1sβ1 + a0sβ0 (10-6-26) 这里 G(s)又称为分数阶传递函数。文献 [38]中给出了求解分数阶线性微分方程的精确解 法，然而该算法即使用计算机实现也有较大难度，所以这里讨论用其他数值算法。 考虑式 (10-6-16)中给出的 Grn¨wald-Letnikov定义，用离散方法可以将其改写成 aD βi t y(t) ' 1 hβi [(t−a)/h]∑ j=0 w (βi) j yt−jh = 1 hβi yt + [(t−a)/h]∑ j=1 w (βi) j yt−jh  (10-6-27) 其中 w(βi)0 可以由下面的递推公式得出 w (βi) 0 = 1, w (βi) j = ( 1− βi + 1 j ) w (βi) j−1, j = 1, 2, · · · (10-6-28) 代入式 (10-6-25)则可以直接推导出微分方程数值解为 yt = 1 n∑ i=0 ai hβi ut − n∑ i=0 ai hβi [(t−a)/h]∑ j=1 w (βi) j yt−jh  (10-6-29) 该算法可以直接由 MATLAB实现，假设将各项的系数和微分阶次用 a和 n向量表 示，并给出等间距时间向量 t和这些点上的输入函数值 u，则可以得出微分方程的数值解 y向量，所以用该函数可以容易地求解分数阶线性微分方程。 function y=glfode(u,t,a,n) h=t(2)-t(1); D=sum(a./[h.^n]); W=[]; y=zeros(length(t),1); w=y; for i=1:length(n) w=1; for j=2:length(t), w(j)=w(j-1)*(1-(n(i)+1)/(j-1)); end W=[W; w]; end for i=2:length(t) for j=1:length(n), A(j)=W(j,2:i)*[y(i-1:-1:1)]; end y(i)=(u(i)-sum(A.*a./[h.^n]))/D; end 【例 10-39】试用数值方法求解下面的分数阶线性微分方程。 D3.5t y(t) + 8D 3.1 t y(t) + 26D 2.3 t y(t) + 73D 1.2 t y(t) + 90D 0.5 t y(t) = 90 sin(t 2) 398 第 10章 数学问题的非传统解法 【求解】由给出的方程可以写出 a和 n向量，从而直接调用编写的 glfode()函数得出该微分方 程的解，用绘图语句可以绘制出输出和输入信号的曲线，如图 10-54所示。为提高得出数值解的 精度，通常需要选择较小的 h值，这里得出的结果精度较高，再进一步减小 h的值，例如选择 h = 0.001，得出的仿真结果与图中给出的结果看不出任何区别。 >> a=[1,8,26,73,90]; n=[3.5,3.1,2.3,1.2,0.5]; t=0:0.002:10; u=90*sin(t.^2); y=glfode(u,t,a,n); subplot(211), plot(t,y); subplot(212), plot(t,u) 0 0.2 0.4 0.6 0 2 4 6 8 10 −100 0 100 y(t) u(t) 图 10-54 方程解及输入信号 在实际计算中对大规模问题有时计算量较大，所以可以考虑采用“短时记忆”效应 的原理对求和式进行近似[38]。 10.6.3.2 非线性分数阶微分方程近似解法 由前面的叙述可见，对未知信号进行分数阶微分数值运算的一种有效途径是采用 Oustaloup算法设计连续滤波器对信号进行滤波处理，另外考虑到该滤波器分子和分母阶 次一致，可能导致在仿真过程中出现代数环，所以应该在其后面再接一个低通滤波器， 将其截止频率设置为 ωh，这样可以建立起如图 10-55a所示的分数阶微分器模块，通过适 当选择频段和阶次则可以较好近似分数阶微分的效果。注意，虽然 Oustaloup算法设计的 滤波器理论上可以求取任意阶次的分数阶微积分，但从数值微积分精度看，该滤波器更 适合求取 1阶以内的分数阶微积分，所以应该将高阶微积分先进行整数阶微积分运算，再 对结果进行滤波处理。 利用 Simulink的模块封装技术[52]，可以将该模型进行封装，得出如图 10-55b所示的 分数阶微分器模块，双击该模块则可以得出如图 10-55c所示的对话框，允许用户填写设 计 Oustaloup滤波器所需的参数，在模块封装初始化栏目应该填写下面的语句，以便在使 用模块前先自动设计出滤波器，并根据阶次正确显示图标。 wb=ww(1); wh=ww(2); G=ousta_fod(gam,n,wb,wh); num=G.num{1}; den=G.den{1}; T=1/wh; str=’Fractional\n’; if isnumeric(gam) if gam>0, str=[str, ’Der s^’ num2str(gam) ]; else, str=[str, ’Int s^{’ num2str(gam) ’}’]; end 10.6 分数阶微积分学及其应用 399 1 Out1T.s+1 1 Transfer Fcn1 num(s) den(s) Transfer Fcn 1 In1 a) 分数阶微分滤波器 c10mfodeFractionalDer s^0.9 b) 封装模块 c) 分数阶微分器参数对话框 图 10-55 分数阶数值微分设计模块 else, str=[str, ’Der s^gam’]; end 在实际仿真过程中，由于搭建起来的系统一般为刚性系统，所以在选择求解算法时 应该选择 ode15s或 ode23tb等，因为这些算法可以保证较高的计算效率和精度。下面将 通过例子演示该模块在分数阶微分方程近似求解中的应用。 【例 10-40】试用滤波器的思想求取例 10-39中分数阶线性微分方程的数值解，并与该例中所用方 法得出的结果进行比较。 【求解】求解分数阶线性微分方程问题不如例 10-39中给出的方法直观，在求解之前，需要引入 辅助变量 z(t) = D0.5t y(t)，这样原来的微分方程可以直接变换成下面的形式 z(t) = sin(t2)− 1 90 [ D3t z(t) + 8D 2.6 t z(t) + 26D 1.8 t z(t) + 73D 0.7 t z(t) ] 根据该方程可以搭建起如图 10-56所示的 Simulink仿真框图，对该框图进行仿真，则可以得出该 微分方程的数值解，将两种方法得出的数值解在同一坐标系下绘制，则得出如图 10-57所示的曲 线，从曲线上基本看不出两者的差别，表明此算法可以以很高精度求解线性方程。 c10mfode1.mdl 1 Out1 1/90 Gain4 73 Gain2 26 Gain1 8 Gain Fractional Int s^{−0.5} Fractional Der s^0.7 Fractional Der s^0.8 Fractional Der s^0.6 90*sin(u[1]^2) Fcn du/dt Derivative2 du/dt Derivative1 du/dt Derivative Clock 图 10-56 微分方程求解的 Simulink框图 400 第 10章 数学问题的非传统解法 0 2 4 6 8 10 −0.2 0 0.2 0.4 0.6 0.8 图 10-57 两种数值解方法比较 【例 10-41】试用近似方法求解下面的分数阶非线性微分方程。 3D0.9y(t) 3 + 0.2D0.8y(t) + 0.9D0.2y(t) + ∣∣2D0.7y(t)∣∣1.5 + 4 3 y(t) = 5 sin(10t) 【求解】根据方程本身，可以容易地写出 y(t)函数的显式表达式为 y(t) = 3 4 [ 5 sin(10t)− 3D 0.9y(t) 3 + 0.2D0.8y(t) + 0.9D0.2y(t) − ∣∣2D0.7y(t)∣∣1.5] 根据得出的 y(t)可以绘制出如图 10-58a所示的仿真模型，从得出的仿真模型可见，信号的各个分 c10mfod2.mdl y(t) 1 Out1 Sine Wave Product 3 Gain4 0.75 Gain3 2 Gain2 0.9 Gain1 0.2 Gain Fractional Der s^0.7 Fractional Der s^0.2 Fractional Der s^0.8 Fractional Der s^0.9 abs(u)^1.5 Fcn 3 Constant a) Simulink仿真模型 0 0.5 1 1.5 2 2.5 3 3.5 4 −2.5 −2 −1.5 −1 −0.5 0 0.5 b) 仿真结果 图 10-58 非线性分数阶微分方程的 Simulink描述及仿真结果 数阶微分信号可以由前面设计的模块获得，因此仿真的精度取决于滤波器对微分的拟合效果，选 择不同的拟合频段和滤波器阶次对求解精度将有一定的影响，图 10-58b对不同的滤波器频段、阶 次组合进