现代科学运算-matlab语言与应用

$\qquad \qquad$ 东北大学 http://www.icourses.cn/home/ (薛定宇)
《高等应用数学问题的MATLAB求解（第四版）》
代码可在matlab r2016b 运行。

10 数学问题的非传统解法

10.01 人工神经网络(上)

人工神经元的数学模型
人工神经元是一个信息处理单元
接受n路输入信号
对其处理得出输出信号
权值、阈值、激活函数(传输函数)

常用的传输函数
Sigmoid函数: $f(x) = \dfrac{2}{1+e^{-2x}} -1 = \dfrac{1 - e^{-2x}}{1 + e^{-2x}}$
对数Sigmoid函数: $f(x) = \dfrac{1}{1 + e^{-x}}$
也可以试用简单的饱和函数和阶跃函数等作为传输函数

x = -2: 0.01: 2;
y = tansig(x);
plot(x, y)

人工神经网络
人工神经网络：由若干各人工神经元按某种规则互连建立起来的网络
不同的神经网络类型
前馈型神经网络、径向基神经网络、Hopfield网络、自组织神经网络
神经网络的应用领域
非线性函数拟合、聚类分析、模式识别与分类

前馈型神经网络
结构：节点、隐层、权值
两层网络(含一个隐层)
$u_j = \sum \limits_{i=1}^{n} w_{ij}x_i + b_{1j}, u^{\prime}_j = F_1(u_j), j = 1, \cdots, p$
$y_j^{\prime} = \sum \limits_{i=1}^{p} v_{ji} u_i^{\prime} + b_2j,y_j= F_2(y_j^{\prime}), j = 1, \cdots, m$

神经网络使用三个步骤
建立网络模型
确定神经网络的结构(类型、隐层、节点个数)
训练网络
由样本数据对神经网络进行训练，得出权值
权值是使得误差最小化计算出的，无物理意义
泛化网络
其作用类似于插值
神经网络本身也是一种数学模型

创建神经网络模型
用于拟合的网络
一般选择前馈网络
选择隐层个数于节点个数
调用fitnet函数建立神经网络框架
net = fitnet([$n_1, n_2, \cdots$])
用于模式识别的网络
选择函数 patternnet

训练样本
随机选择训练用样本
$\begin{array}{l} x_{11} & x_{12} & \cdots & x_{1n} \\ x_{21} & x_{22} & \cdots & x_{2n} \\ \vdots & \vdots & \vdots & \vdots \\ x_{N1} & x_{N2} & \cdots & x_{Nn} \end{array} \begin{array}{l} \Rightarrow \\ \Rightarrow \\ \vdots \\ \Rightarrow \end{array} \begin{array}{l}\hat y_{11} & \cdots & \hat y_{1m} \\ \hat y_{21} & \cdots & \hat y_{2m} \\ \vdots & \vdots & \\ \hat y_{N1} & \cdots & \hat y_{Nm} \end{array}$
通过训练更新权值
$W_{ij}^{l+1} = W_{ij}^{l} + \beta e_j^{l}\alpha_i^l , V_{jt}^{l+1} = V_{jt}^{l} + \alpha d_t^l \gamma_j^l$

训练的目标函数
$\min \limits_{W,V} \sum \limits_{l=1}^{N} \sum \limits_{i=1}^{m}(y_{li} - \hat y_{li})^2$
优化目标函数，由样本得出权值
matlab: net = train(net, X, Y)

神经网络的泛化
将输入信号传给神经网络，计算出输出信号，称为神经网络的泛化
matlab: $Y_1 = sim(net, X_1), Y_1 = net(X_1)$
其它函数：
显示网络结构: view(net)
训练的效果：I = perform(net, $Y_1, X_1$)

例10-19 神经网络数据拟合

% 生成一组数据
x = 0: .5: 10;
y = 0.12*exp(-0.213*x) + 0.54* ...
    exp(-0.17*x).*sin(1.23*x);
% 试用神经网络对其进行拟合，样本点生成
% 选择1个隐层，隐层节点选择为5
x0 = [0: 0.1: 10];
y0 = 0.12*exp(-0.213*x0) + ...
    0.54*exp(-0.17*x0).*sin(1.23*x0);
net = fitnet(5);
net = train(net, x, y);
y1 = net(x0);
plot(x, y, 'o', x0, y0, x0, y1, ':');

神经网络的其它信息
显示结构: view(net)
训练后的net是结构体数据

w1 = net.IW{1}
w2 = net.LW{2, 1}

改变参数对神经网络拟合的影响
是不是节点个数越多越好？节点个数20

net = fitnet(20);
net = train(net, x, y);
y1 = net(x0);
plot(x, y, 'o', x0, y0, x0, y1, ':');

再次运行可能得出不同结果
生成Simulink模块 : gensim(net)

10.02 人工神经网络(下)

例10-20 二元函数的拟合
散点样本数据生成: $z = (x^2 - 2x)e^{-x^2-y^2-xy}$
如果需要，也可以生成网格点

x = -3+6*rand(1, 200);
y = -2 + 4*rand(1, 200);
z = (x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
% 选择两个隐层，每层10个节点
net = fitnet([10, 10]);
net = train(net, [x; y], z);
[x2, y2] = meshgrid(-3: .1: 3, -2: .1: 2);
x1 = x2(:)';
y1 = y2(:)';
z1 = net([x1; y1]);
z2 = reshape(z1, size(x2));
surf(x2, y2, z2)

不同网络设置下的拟合

% 多选择节点
net = fitnet([20 20]);
[net, b] = train(net, [x;y], z);
z1 = net([x1; y1]);
z2 = reshape(z1, size(x2));
surf(x2, y2, z2);

% 4层网络的拟合
net = fitnet([7 7 7]);
[net, b] = train(net, [x; y], z);
z1 = net([x1; y1]);
z2 = reshape(z1, size(x2));
surf(x2, y2, z2);

例10-21 一个反例
第8章已知5个点的正弦函数拟合
已知: $x = [0, 0.4, 1, 2, \pi]$
函数计算: $y = \sin x$
直接建立网络并观察拟合效果

x = [0, 0.4, 1, 2, pi];
y = sin(x);
x1 = 0: 0.01: pi;
net = fitnet(10);
net = train(net, x, y);
y1 = net(x1);
plot(x1, y1, x, y, 'o');

径向基网络结构与应用
径向基函数(RBF)的数学描述
$\psi(x) = e^{-b \Vert x-c \Vert} = e^{-b(x-c)^T(x-c)}$
其中
c为聚类中心点
b > 0 为调节聚类效果的参数
神经网络工具箱中，radbas() 函数可以计算出标准径向基函数 $y_i = e^{-x_i^2}$ 的曲线参数

例10-21 径向基函数曲线
绘制不同参数(c, b)下的径向基函数曲线
取中心点 c = -2, 0, 2,并假设 b = 1

x = -4: 0.1: 4;
cc = [-2, 0, 2];
b = 1;
for c = cc
    y = exp(-b*(x-c).^2);
    plot(x, y);
    hold on;
end;

若隐层传输函数$F_1(x)$为径向基函数，输出层 $F_2(x)$为线性函数, 则此结构的网络称为径向基网络

径向基网络建模、训练与泛化
径向基网络的训练方式不是采用方向误差传播实现的，故不是BP网络
径向基网络的使用
由函数newrbe()建立并训练神经网络
net = newrbe(X, Y)
用sim()函数可以泛化神经网络
$Y_1 = sim(net, X_1); Y_1 = net(X_1)$

例10-24 一元函数的RBF拟合
用径向基神经网络重新拟合一元函数

x = 0: .5: 10;
y = 0.12*exp(-0.213*x) + ...
    0.54*exp(-0.17*x).*sin(1.23*x);
% 径向基网络的拟合(对本例优于fitnet)
x0 = [0: 0.1: 10];
y0 = 0.12*exp(-0.213*x0) + ...
    0.54*exp(-0.17*x0).*sin(1.23*x0);
net = newrbe(x, y);
y1 = sim(net, x0);
plot(x, y, 'o', x0, y0, x0, y1, ':');

例10-25 二元函数的RBF拟合
二元函数 $z = f(x, y) = (x^2-2x)e^{-x^2-y^2-xy}$

% 生成样本数据
x = -3+6*rand(1, 200);
y = -2+4*rand(1, 200);
z = (x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
% 曲面拟合
net = newrbe([x; y], z);
[x2, y2] = meshgrid(-3: .1: 3, -2: .1: 2);
x1 = x2(:)';
y1 = y2(:)';
z1 = sim(net, [x1; y1]);
z2 = reshape(z1, size(x2));
surf(x2, y2, z2)

例10-26 正弦的RBF网络拟合
第8章已知5个点的正弦函数拟合
已知: $x = [0, 0.4, 1, 2, \pi]$
函数计算: $y = \sin x$
网络fitnet失效, 改用RBF网络拟合

x = [0, 0.4, 1, 2, pi];
y = sin(x);
x1 = 0: 0.01: pi;
y0 = sin(x1);
net = newrbe(x, y);
y1 = net(x1);
plot(x1, y1, x1, y0, x, y, 'o');

神经网络界面
启动神经网络工具箱的图形用户界面的函数调用格式 nntool
该界面可以建立所需的神经网络模型，并可以由已知数据对该网络进行训练、泛化

10.03 进化算法与全局最优化(上)

遗传算法的基本概念
遗传算法是基于进化论，在计算机上模拟生命进化机制
它根据适者生存、优胜劣汰自然进化规则搜索和计算问题的解，并行的搜索方法
美国 Michigen 大学 John Holland 1975 年提出
更可能获得全局最优解

遗传算法的进化机制
基本概念
种群(初始选择的个体)
个体
染色体–编码(二进制编码、实数编码)
适应度函数–目标函数
生成下一代：复制、交叉、变异
代

遗传算法和传统优化算法比较
传统最优化算法从初始点开始搜索，遗传算法从种群(若干个初值)开始进行搜索
遗传算法并不依赖于目标函数导数信息或其它辅助信息来进行最优解搜索
遗传算法采用的是概率型规则而不是确定性规则，所以每次得出的结果不一定完全相同，有时甚至可能会有较大差异。

遗传算法的matlab函数
matlab工具箱的全局优化工具箱
早期版本称为遗传算法于直接搜索工具箱
大量其它寻优算法，包括粒子群优化算法、模拟退火算法与模式搜索算法等
一些方法号称能够求解有约束最优化问题
英国Sheffield大学的遗传算法工具箱
matlab central下有大量工具箱
GAOT: 最早的遗传算法最优化工具箱

GAOT工具箱中的函数调用格式：
最大值计算、目标函数的不同格式
[a, b, c] = gaopt(bound, fun)
matlab 全局优化工具箱
[x, f, flag, out] = ga(fun, n, opts)
打开遗传算法程序界面：gatoool()

例10-28 改进的Rastrigin多峰函数无约束优化
目标函数: $f(x_1, x_2) = 20 + (x_1/30 - 1)^2 + (x_2/20-1)^2 + (x_2-2x_3)^4 + 10(x_1-4)^4 - 10[cos(x_1/30-1)\pi + cos(x_2/20-1)\pi]$
目标函数的曲面表示
多峰函数，全局最优(30, 20)点
传统方法容易陷入局部最优解

ezsurf(['20+(x1/30-1)^2+(x2/20-1)^2', ...
    '-10*(cos(pi*(x1/30-1))', ...
    '+cos(pi*(x2/20-1)))'], [-100, 100]);

ezsurf(['20+(x1/30-1)^2+(x2/20-1)^2', ...
    '-10*(cos(pi*(x1/30-1))', ...
    '+cos(pi*(x2/20-1)))'], [-100, 100]);
view(0, 0);

GAOT工具箱使用
目标函数的描述
GAOT求解最大值问题
目标函数的格式与一般最优化不一样
$f(x_1, x_2) = 20 + (x_1/30 -1)^2 + (x_2/20-1)^2 + (x_2-2x_3)^4 + 10(x_1-x_4)^4 - 10[cos(x_1/30-1)\pi + cos(x_2/20-1)\pi]$

function [S, f] = c10mga7(S, options)
    s = S(1:2);
    f = -20-(x(1)/30-1)^2-(x(2)/20-1)^2+...
        10*(cos(pi*(x(1)/30-1))+cos(pi*(x(2)/20-1)));

GAOT求解

range=[-100 100; -100 100];
x = gaopt(range, 'c10mga7')

全局优化工具箱求解

f = @(x)20+(x(1)/30-1)^2+(x(2)/20-1)^2 ...
    -10*(cos(pi*(x(1)/30-1))+...
    cos(pi*(x(2)/20-1)));
x = ga(f, 2)

全局搜索函数(第6章)

[x, f0] = fminunc_global(f, -100, 100, 2, 50)

遗传算法在有约束最优化问题中的应用
GAOT工具箱不能直接用于有约束最优化问题的求解
对于有约束最优化问题，应该将其转变为无约束最优化问题
GAOT求解的是最大值问题
不等式约束则可以通过惩罚函数方法转换到目标函数中，如将目标函数设置成-20000
等式约束，则可以通过方程求解的方式将其中若干个自变量用其它自变量表示

例10-32 线性规划问题
线性规划问题是凸问题，没有必要用遗传算法
试用遗传算法求解线性规划问题
$\min \quad (x_1 + 2x_2 + 3x_3)$
$x \; s.t. \left \{ \begin{array}{l}-2x_1 + x_2 + x_3 \leq 9 \\ -x_1 + x_2 \geq -4 \\ 4x_1 -2x_2 -3x_3 = -6 \\ x_{1, 2} \leq 0, x_3 \geq 0 \end{array} \right.$
对于等式，用$x_1 和 x_2 表示 x_3$
$x_3 = \dfrac{6 + 4x_1 - 2x_2}{3}$

% 描述目标函数
function [sol, y] = c10mga4(sol, options)
    x = sol(1:2);
    x = x(:);
    x(3) = (6+4*x(1)-2*x(2))/3;
    y1 = [-2 1 1 ]*x;
    y2 = [-1 1 0]*x;
    if (y1 > 9 | y2 < -4 | x(3) < 0
        y = -100;
    else
        y = -[1 2 3]*x;
    end;

% 运用遗传算法求解
[a, b, c] = gaopt([-1000 0; -1000 0], ...
     'c10mga4', [], [], [], 'maxGenTerm', 1000);
x = a(1:2);
x(3) = (6+4*x(1)^2-2*x(2)/3,
f = a(3)

常见求解方法
用线性规划方法求解

f = [1 2 3];
A = [-2 1 1; 1 -1 0];
B = [9; 4];
Aeq = [4 2 -3];
Beq = -6;
x = linprog(f, A, B, Aeq, Beq, [-inf; -inf; 0], [0; 0; inf]);
x', fm = f*x

linprog()得出的结果更精确
线性规划是凸问题，没有必要用遗传算法
如果等式约束复杂，很难求解并转换

10.04 进化算法与全局最优化(下)

全局优化工具箱函数求解
有约束优化问题的直接求解
[x, a, key] = ga(f, n, A, B, $A_{eq}, B_{eq}, x_m, x_M, f_1$, IntCon)
[x, a, key] = ga(problem)
看似强大，但等式约束有时效果似乎不佳
求解的是最小值，而不像gaopt的最大值
得出的是可行解

例10-28 重解线性规划问题
用ga函数直接求解线性规划问题
$\min \quad (x_1 + 2x_2 + 3x_3)$
$x \; s.t. \left \{ \begin{array}{l}-2x_1 + x_2 + x_3 \leq 9 \\ -x_1 + x_2 \geq -4 \\ 4x_1 -2x_2 -3x_3 = -6 \\ x_{1, 2} \leq 0, x_3 \geq 0 \end{array} \right.$

f = @(x)[1 2 3]*x(:);
A = [-2 1 1; 1 -1 0];
B = [9; 4];
Aeq = [4 -2 -3];
Beq = -6;
xm = [-inf; -inf; 0];
xM = [0; 0; inf];
[x a k] = ga(f, 3, A, B, Aeq, Beq, xm, xM)

例10-29 整数规划问题
非线性整数规划
$\min \quad 2y_1^2/16 + y_2^2/100 - 4y_1 - y_2$
$y \; s.t. \; \left \{ \begin{array}{l}y_1^2/16 - 6y_1/4 + y_2/10 - 11 \leq 0 \\ -y_1y_2/40 + 3y_2/10 + e^{y_1/4-3} - 1 \leq 0 \\ y_2 \geq 30 \end{array} \right.$
直接求解：
[x, a, key] = ga(f, n, A, B, $A_{eq}, B_{eq}, x_m, x_M, f_1$, IntCon)

f = @(y)2*y(1)^2/16+y(2)^2/100-4*y(1)-y(2);
xm = [-200; 30];
xM = [200; 200];
x = ga(f, 2, [], [], [], [], xm, xM, @c6mdisp, [1, 2])

粒子群优化算法与求解
粒子群优化算法是一种进化算法，该算法是受生物界鸟群觅食的启发提出的搜索算法
假设某个区域内有一个全局最优点，和位于随机初始位置的粒子
每一个粒子有到目前为止自己的个体最优值$p_{i,b}$
整个粒子群有到目前为止群体的最优值$g_b$

粒子群算法核心公式
每个粒子需要更新自己的速度和位置
$v_i(k+1) = \phi(k)v_i(k) + \alpha_1\gamma_{1i}(k)[p_{i,b} - x_i(k)] + \alpha_2\gamma_{2i}(k)[g_b - x_i(k)]$
$x_i(k+1) = x_i(k) + v_i(k+1)$
其中,
$\gamma_{1a}, \gamma_{2a}为[0, 1]区间内均匀分布的随机数$
$\phi(k)为惯量函数$
$\alpha_1,\alpha_2为加速常数$

粒子群算法的求解函数
matlab2014b新的粒子群优化函数
[x, $f_m$, key] = particleswarm(problem),
[x, $f_m$, key] = particleswarm(f, n, $x_m, x_M$, opts)
无约束最优化问题
不能求解真正的有约束问题
算法的效率比其它已知粒子群优化算法函数高
早期版本可以搜索其它现成工具：例如，PSOt工具箱

例10-29 重新求解无约束最优化问题
改进的Rastrigin函数
$f(x_1, x_2) = 20 + (x_1/30-1)^2+(x_2/20-1)^2-10[\cos(x_1/30-1)\pi + \cos(x_2/20-1)\pi]$

% 粒子群算法直接求解
f = @(x)20+(x(1)/30-1)^2+(x(2)/20-1)^2- ...
    10*(cos(pi*(x(1)/30-1))+ ...
    cos(pi*(x(2)/20-1)));
[x g] = particleswarm(f, 2)

全局优化工具箱其它函数
无约束优化问题求解
模拟退火方法 simulannealbnd
x = simulannealbnd(f, $x_0, x_m, x_M$, opts)
x = simulannealbnd(problem)
有约束最优化问题
模式搜索方法 patternsearch
[x, $f_{opt}$, flag, c] = patternsearch(problem)
[x, $f_{opt}$, flag, c] = patternsearch(fun, $x_0, A, B, A_{eq}, B_{eq}, x_m, x_M, CFun, OPT, p_1, p_2, \cdots$)

无约束优化算法比较
求解改进的Rastrigin函数最优化问题
每种算法运行100次，评价效果
$\begin{array}{l}\hline solver \; function & toolbox & success \; rate & time & error \; norm \\ \hline fminunc\_global() & author's \; own & 97\% & 112 & 3.1 \times 10^{-12} \\ga() & Global & 95\% & 6.3 & 0.039 \\ gaopt() & GAOT \; Toolbox & 100\% & 21.2 & 3.28\times10^{-8} \\ particleswarm() & Global & 95\% & 2.94 & 4.54\times 10^{-8} \\ patternsearch() & Global & 34\% & 3.6 & 1.4\times 10^{-9} \\ simulannealbnd() & Global & 57\% & 41.9 & 0.14 \\ ievyPSO() & free & 100\% & 27 & 3.71\times 10^{-10} \\ \hline \end{array}$

有约束问题比较
比较的函数
全局最优化工具箱函数
遗传算法：ga()
模式搜索:patternsearch()
第六章的:fminunc_global()
比较的结果：运行100次
fmincon_global 成功 100%, 总耗时225秒
ga(), patternsearch()函数均低于3%

例10-30 第6章的例题
最优化问题:
$\min \quad x_5$
$x \; s.t. \; \left \{ \begin{array}{l}x_3 + 9.625x_1 x_4 + 16 x_2 x_4 + 16 x_4^2 + 12 - 4x_1 - x_2 - 78 x_4 = 0 \\ 16x_1 x_4 + 44 - 19 x_1 - 8 x_2 - x_3 - 24 x_4 = 0 \\ -0.25x_5 - x_1 \leq -2.25 \\ x_1 - 0.25x_5 \leq 2.25 \\ -0.5 x_5 - x_2 \leq -1.5 \\ x_2 - 0.5x_5 \leq 1.5 \\ -1.5x_5 - x_3 \leq -1.5 \\ x_3 - 1.5x_5 \leq 1.5 \end{array} \right.$

clear P;
P.objective = @(x)x(5);
P.nonlcon = @c6exnls;
P.Aineq = [-1 0 0 0 -0.25; 1 0 0 0 -0.25; 0 -1 0 0 -0.5;
    0 1 0 0 -0.5; 0 0 -1 0 -1.5; 0 0 1 0 -1.5];
P.bineq = [-2.25; 2.25; -1.5; 1.5; -1.5; 1.5];
P.options = gaoptimset;
P.solver = 'patternsearch';
P.X0 = 10*rand(5, 1);
patternsearch(P)

10.05 分数阶微积分定义与性质

分数阶微积分学概述
第三章的整数阶微积分问题 $d^n y /d x^n$
分数阶微积分理论已经有300年历史
早期主要侧重理论研究
近年来在很多领域已经开始应用
自动控制领域出现了分数阶控制理论等新分支
分数阶信号处理

为什么需要分数阶微积分
例：正弦函数的导数是什么？高阶导数是什么？
$\dfrac{d^n}{d t^n} \sin t = \sin(t + \dfrac{n \pi}{2})$
分数阶导数提供信息

n0 = 0: 0.1: 1.5;
t = 0: 0.2: 2*pi;
Z = [];
for n = n0
    Z = [Z; sin(t+n*pi/2)];
end;
surf(t, n0, Z)

分数阶微积分可能揭示出整数阶微积分视角下看不到的东西

分数阶微积分定义
四种常用的定义、互相关联
分数阶Cauchy积分公式
Gr$\ddot u$nwald-Letnikov分数阶微积分定义
Riemann-Liouville分数阶微积分定义
Caputo分数阶微积分定义

分数阶Cauchy积分公式
数学描述
$\mathscr{D}^{\gamma}f(t) = \dfrac{\Gamma(\gamma+1)}{2\pi j} \oint_{C} \dfrac{f(\tau)}{(\tau-t)^{\gamma+1}} d \tau$
其中, C为包围f(t)单值与解析开区域的光滑曲线
不易计算，不考虑

Gr$\ddot u$nwald-Letnikov分数阶微积分定义
数学描述
$_a\mathscr{D}_t^{\alpha} f(t) = \lim \limits_{h \to 0} \dfrac{1}{h^{\alpha}} \sum \limits_{j=0}^{[(t-a)/h]} (-1)^j \dbinom \alpha j f(t-jh)$
$其中，\dbinom \alpha j 为二项式系数 w_j^{(\alpha)} = (-1)^j \dbinom \alpha j$
$w_0^{(\alpha)} = 1, w_j^{(\alpha)} = \Big(1 - \dfrac{\alpha + 1}{j} \Big) w_{j-1}^{(\alpha)}, j = 1, 2, \cdots$

Riemann-Liouville分数阶微积分定义
积分数学描述
$_a \mathscr{D}_t^{-\alpha} f(t) = \dfrac{1}{\Gamma(\alpha)} \int_a^t (t-\tau)^{\alpha-1} f(\tau) d \tau$
$其中, 0 < \alpha < 1, 且 a 为初始时间$
分数阶微分定义的数学描述
$_a\mathscr{D}_t^{\beta} f(t) = \dfrac{d^n}{dt^n}\Big[ _a\mathscr{D}_t^{-(n-\beta)} f(t) \Big] = \dfrac{1}{\Gamma(n-\beta)} \dfrac{d^n}{dt^n} \Big[ \int_a^t \dfrac{f(\tau)}{(t-\tau)^{\beta - n+1}} \Big] d \tau$

Caputo分数阶微积分定义
微分的数学描述
$_0\mathscr{D}_t^{\alpha} y(t) = \dfrac{1}{\Gamma(1-\alpha)} \int_0^t \dfrac{y^{m+1}(\tau)}{(t-\tau)^{\alpha}} d \tau$
$其中，\alpha = m + \gamma, m 为整数， 0 \leq \gamma \leq 1$
$当 \gamma < 0 时, Caputo 分数阶积分的数学描述与Riemann-Liouville积分定义完全一致$
$_0 \mathscr{D}_t^{\gamma} y(t) = \dfrac{1}{\Gamma(-\gamma)} \int_0^t \dfrac{y(\tau)}{(t-\tau)^{1 + \gamma}} d \tau$

不同分数阶微积分定义的关系
Gr$\ddot u$nwald-Letnikov，Riemann-Liouville 等效
Caputo定义和Riemann-Liouville定义
$0 < \alpha < 1, \quad _{t_0}^{C} \mathscr{D}_t^{\alpha} y(t) = _{t_0}^{RL} \mathscr{D}_{t}^{\alpha}(y(t) - y(t_0))$
$\qquad \qquad = _{t_0}^{RL}\mathscr{D}_t^{\alpha} y(t) - \dfrac{y(t_0)(t-t_0)^{-\alpha}}{\Gamma(1-\alpha)}$
更一般地
$_{t_0}^{C} \mathscr{D}_t^{\alpha} y(t) = _{t_0}^{RL} \mathscr{D}_t^{\alpha} y(t) - \sum \limits_{k=0}^{m-1} \dfrac{y^{(k)}(t_0)}{\Gamma(k-\alpha+1)}(t-t_0)^{k-\alpha}$

分数阶微积分的性质
$解析函数f(t), _0\mathscr{D}_t^{\alpha} f(t) 对 t 和 \alpha 都是解析的$
$\alpha 为整数时，分数阶微分与整数阶微分的值完全一致，且 _0\mathscr{D}_t^{0} f(t) = f(t)$
分数阶微积分算子为线性的，对常数a、b
$_0\mathscr{D}_t^{\alpha}[a f(t) \pm b g(t)] = a \; _0\mathscr{D}_t^{\alpha} f(t) \pm b \; _0\mathscr{D}_t^{\alpha} g(t)$
分数阶微积分算子满足交换律
$_0\mathscr{D}_t^{\alpha} \Big[_0\mathscr{D}_t^{\beta} f(t) \Big] = _0\mathscr{D}_t^{\beta} \Big[ _0\mathscr{D}_t^{\alpha} f(t) \Big] = _0\mathscr{D}_t^{\alpha+\beta} f(t)$

函数分数阶导数的Laplace变换
函数的分数阶积分表达式的Laplace变换
$\mathscr{L} \Big[ \mathscr{D}_t^{-\gamma} f(t) \Big] = s^{-\gamma} \mathscr{L}[f(t)]$
Riemann-Liouville微分的Laplace变换为
$\mathscr{L}\Big[ \; _0^{RL} \mathscr{D}_t^{\alpha} f(t) \Big] = s^{\alpha} \mathscr{L}[f(t)] - \sum \limits_{k=1}^{n-1} s^k [ _0\mathscr{D}_t^{\alpha-k-1} f(t)]_{t = 0}$
Caputo微分的Laplace变换
$\mathscr{L} \Big[ \; _0^{C}\mathscr{D}_t^{\gamma} f(t) \Big] = s^{\gamma}F(s) - \sum \limits_{k=0}^{n-1} s^{\gamma-k-1} f^{(k)}(0)$

10.06 分数阶微积分数值计算

原函数或其采样点已知时
用Grunwald-Letnikov定义直接计算
Grunwald-Letnikov定义下高精度数值求解
Gaputo定义的高精度数值解
原函数未知时
通过高阶整数阶滤波器去模拟分数阶微积分行为
Oustaloup滤波器逼近 Grunwald-Letnikov微积分
Caputo微积分的滤波器逼近

用Grunwald-Letnikov定义求解分数阶微分
数学描述
$_a\mathscr{D}_t^{\alpha}f(t) = \lim \limits_{h \to 0} \dfrac{1}{h^{\alpha}} \sum \limits_{j=0}^{[(t-a)/h]} (-1)^{j} \dbinom \alpha j f(t-jh)$
$\qquad \qquad \approx \dfrac{1}{h^{\alpha}} \sum \limits_{j=0}^{[(t-a)/h]} w_j^{(\alpha)} f(t - jh)$
其中：
$w_0^{(\alpha)} = 1, w_j^{(\alpha)} = \Big(1 - \dfrac{\alpha + 1}{j} \Big) w_{j-1}^{(\alpha)}, j = 1, 2, \cdots$

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;

函数调用格式： $y_1$ = glfdiff(y, t, $\gamma$)

高精度分数阶微积分数值计算
高精度Grunwald-Letnikov微积分的数值计算
算法精度o($h^p$)
matlab函数 [$y_1, t$] = glfdiff9(y, t, $\gamma$, p)
Caputo微积分的高精度数值计算
$y_1$ = caputo9(y, t, $\alpha$, p)
代码下载–FOTF工具箱
http://cn.mathworks.com/matlabcentral/fileexchange/60874-fotf-toolbox

例10-47 常数的微积分
阶跃函数的导数和积分是什么？
整数阶的积分与微分
对不同的阶次

t = 0: 0.01:1.5;
gam = [-1, -0.5, 0.3, 0.5, 0.7];
y = ones(size(t));
dy = [];
for a = gam
    dy = [dy; glfdiff(y, t, a)];
end;
plot(t, dy)

例10-50 正弦信号的分数阶微分
已知函数为$f(t) = \sin(3t+1)$
0.3阶Riemann-Liouville导数计算
其0.3，1.3，2.3 阶Caputo 导数计算

t = 0: 0.01: pi;
y = sin(3*t + 1);
y1 = glfdiff9(y, t, 0.3, 4);
y2 = caputo9(y, t, 0.3, 4);
plot(t, y1, t, y2, '--', t, d, ':');
figure;
y1 = caputo9(y, t, 1.3, 4);
y2 = caputo9(y, t, 2.3, 4);
plotyy(t, y1, t, y2);

Caputo微积分定义的数值计算的精度评价
计算指数函数exp(-t)的0.6阶Caputo导数(h = 0.01)
解析解: $y_0(t) = _0^C \mathscr{D} _t^{\alpha} e^{-t} = \lambda^{q} t^{\gamma} E^{1, \gamma+1}(\lambda t) = -t^{0.4} E_{1, 1.4}(-t)$
Mittag-Leffler函数
y = ml_func$([\alpha, \beta], z) \qquad E_{\alpha, \beta}(z) = \sum \limits_{k=0}^{\infty} \dfrac{z^k}{(\alpha k + \beta)}$

t0 = 0.5: 0.5: 5;
q = 1;
gam = q - 0.6;
y0 = -t0.^0.4.*ml_func([1,1.4], -t0, 0, eps);
t = 0: 0.01: 5;
y = exp(-t);
ii = [51:50:501];
T = t0';
for p=1:6
    y1 = caputo9(y, t, 0.6, p);
    T = [T [y1(ii)-y0']];
end;

事先未知函数的分数阶导数
可以考虑采用高阶整数阶滤波器逼近分数算子
分数阶微分环节的频域响应(Bode图)
可以考虑引入整数滤波器去逼近
渐近线逼近，实际逼近效果会更好

Oustaloup滤波算法
如果函数f(t)的表达式不是事先已知
连续滤波器传递函数模型: Oustaloup算法
感兴趣区间($\omega_b, \omega_h$)
$G_f(s) = k\prod \limits_{k=1}^{N} \dfrac{s + \omega_k^{\prime}}{s + \omega_k}$
$\omega_k^{\prime} = \omega_b \omega_{u}^{2k-1-\gamma)/N}, \omega_k = \omega_b \omega_u^{(2k-1+\gamma)/N}, K = \omega_h^2,$
$\omega_u = \sqrt{\omega_h/\omega_b}$

function G = ousta_fod(gam, N, wb, wh)
    k = 1:N;
    wu = sqrt(wh/wb);
    wkp = wb*wu.^((2*k-1-gam)/N);
    wk = wb*wu.^((2*k-1+gam)/N);
    G = zpk(-wkp, -wk, wh^gam);
    G = tf(G);

函数调用格式： G = ousta_fod($\gamma, N, \omega_b, \omega_h$)

例10-51 由滤波器计算微分
函数$f(t) = e^{-t} \sin (3t + 1)$
感兴趣频率区间 $\omega_b = 0.01, \omega_h = 1000 rad/sec$
分数阶次 0.5， 5阶滤波器

G = ousta_fod(0.5, 5, 0.01, 1000),
bode(G)

绘制分数阶微分曲线

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);

非零初值Caputo微积分的Oustaloup滤波器逼近
标准Oustaloup滤波器可以用RL微积分逼近
直接方法构造Caputo微分器较麻烦，不宜采用
非零初值信号的Caputo微积分逼近可以根据两个定理
定理1：RL积分器
$_{t_0}^{C}\mathscr{D}_t^{\gamma} y(t) = _{t_0}^{RL}\mathscr{D}_t^{-(\lceil \gamma \rceil - \gamma)}\Big[ y^{(\lceil \gamma \rceil)} (t) \Big]$
$_{t_0}^{C}\mathscr{D}_t^{2.3}y(t) = _{t_0}^{RL}\mathscr{D}_t^{-0.7} \Big[y^{\prime \prime \prime}(t) \Big]$
定理2: RL微分器
$_{t_0}^{RL}\mathscr{D}_{t}^{\lceil \gamma \rceil - \gamma} \Big[ \; _{t_0}^{C} \mathscr{D}_{t}^{\gamma} y(t) \Big] = y^{(\lceil \gamma \rceil)}(t)$
$_{t_0}^{RL}\mathscr{D}_{t}^{0.7} \Big[ \; _{t_0}^{C}\mathscr{D}_t^{2.3}y(t) \Big] = y^{\prime \prime \prime}(t)$

10.07 分数阶微分方程的数值解

零初值分数阶线性微分方程的标准型
分数阶线性微分方程标准型:
$a_n \mathscr{D}_t^{\beta_n} y(t) + a_{n-1}\mathscr{D}_t^{\beta_{n-1}} y(t) + \cdots + a_1 \mathscr{D}_t^{\beta_1} y(t) + a_0 \mathscr{D}_t^{\beta_0} y(t) = \hat u(t)$
如果右侧不是u(t),可以通过下式计算
$\hat u(t) = b_1\mathscr{D}_t^{\gamma_1} u(t) + b_2 \mathscr{D}_t^{\gamma_2} u(t) + \cdots + b_m \mathscr{D}_t^{\gamma_m} u(t)$
假设 $\beta_n > \beta_{n-1} > \cdots > \beta_1 > \beta_0 > 0$
如果有负阶次，应该事先变换 $z(t) = \mathscr{D}_t^{\beta_{\alpha}} y(t)$

线性微分方程的闭式解
由Grunwald-Letnikov定义
$_a\mathscr{D}_t^{\beta_i} y(t) \approx \dfrac{1}{h^{\beta_i}} \sum \limits_{j=0}^{[(t-a)/h]} \omega_j^{(\beta_i)} y_{t-jh} = \dfrac{1}{h^{\beta_i}} \Big[ y_t + \sum \limits_{j=1}^{[(t-a)/h]} \omega_j^{(\beta_i)} y_{t-jh} \Big]$
$\omega_0^{(\beta_i)} = 1, \omega_j^{(\beta_i)} = (1 = \dfrac{\beta_i + 1}{j}) \omega_{j-1}^{(\beta_i)}, j = 1, 2, \cdots$
微分方程的闭式数值解为
$y_t = \dfrac{1}{\sum \limits_{i=0}^{n} \dfrac{a_i}{h^{\beta_i}}} \Big[ \hat u_t - \sum \limits_{i=0}^{n} \dfrac{a_i}{h^{\beta_i}} \sum \limits_{j=1}^{[(t-a)/h]} \omega_j^{(\beta_i)} y_{t-jh} \Big]$
直接求解: y = fode_sol($a, n_a, b, n_b, u, t$)

function y = fode_sol(a, na, b, nb, u, t)
    h = t(2) - t(1);
    D = sum(a./[h.^na]);
    nT = length(t);
    vec = [na nb];
    W = [];
    u = u(:).';
    D1 = b(:)./h.^nb(:);
    nA = length(a);
    y1 = zeros(nT,1);
    W = ones(nT, length(vec));
    for j=2:nT, W(j,:)=W(j-1,:).*(1-(vec+1)/(j-1)); end
    for i=2:nT, 
        A=[y1(i-1:-1:1)]'*W(2:i,1:nA); y1(i)=(u(i)-sum(A.*a./[h.^na]))/D;
    end
    for i=2:nT, y(i)=(W(1:i,nA+1:end)*D1)'*[y1(i:-1:1)]; end

例10-52 分数阶微分方程的数值解
求解零初值分数阶线性微分方程
$\mathscr{D}_t^{3.5} y(t) + 8 \mathscr{D}_t^{3.1} y(t) + 26 \mathscr{D}_t^{2.3} y(t) + 73 \mathscr{D}_t^{1.2} y(t) + 90 \mathscr{D}_t^{0.5} y(t) = 90 \sin t^2$
选择计算步长为 h = 0.002

a = [1, 8, 26, 73, 90];
n = [3.5, 3.1, 2.3, 1.3, 0.5];
t = 0: 0.002:10;
u = 90*sin(t.^2);
y = fode_sol(a, n, 1, 0, u, t);
subplot(211);
plot(t, y);
subplot(212);
plot(t, u);

分数阶微分方程的高精度求解
Riemann-Liouville 方程的高精度求解
y = fode_sol9($a, n_a, b, n_b, u, t, p$)
Caputo方程的高精度求解
y = fode_caputo9($a, n_a, b, n_b, y_0, u, t, p$)
非线性Caputo方程的高精度求解
预估值：[y, t] = nlfep(fun, $\alpha, y_0, t_n, h, p, \epsilon$)
矫正解: y = nlfec(fun, $\alpha, y_0, y_p, t, p, \epsilon$)

function y=fode_sol9(a,na,b,nb,u,t,p)
h=t(2)-t(1);  n=length(t); vec=[na nb]; u=u(:);
g=double(genfunc(p)); t=t(:); W=[];
for i=1:length(vec), W=[W; get_vecw(vec(i),n,g)]; end
D1=b(:)./h.^nb(:); nA=length(a); y1=zeros(n,1);
W=W.'; D=sum((a.*W(1,1:nA))./[h.^na]);
for i=2:n
   A=[y1(i-1:-1:1)]'*W(2:i,1:nA); y1(i)=(u(i)-sum(A.*a./[h.^na]))/D;
end
for i=2:n, y(i)=(W(1:i,nA+1:end)*D1)'*[y1(i:-1:1)]; end

function w=get_vecw(gam,n,g),
p=length(g)-1; b=1+gam; g0=g(1);
w=zeros(1,n); w(1)=g(1)^gam;
for m=2:p, M=m-1; dA=b/M;
   w(m)=-[g(2:m).*((1-dA):-dA:(1-b))]*w(M:-1:1).'/g0;
end
for k=p+1:n, M=k-1; dA=b/M;
   w(k)=-[g(2:(p+1)).*((1-dA):-dA:(1-p*dA))]*w(M:-1:(k-p)).'/g0;
end

function g=genfunc(p)
a=1:p+1; A=rot90(vander(a)); g=(1-a)*inv(sym(A'));

function dy=glfdiff9(y,t,gam,p)
if strcmp(class(y),'function_handle'), y=y(t); end, y=y(:);
h=t(2)-t(1); t=t(:); n=length(t); u=0; du=0; r=(0:p)*h; 
R=sym(fliplr(vander(r))); c=double(inv(R)*y(1:p+1)); 
for i=1:p+1, u=u+c(i)*t.^(i-1);
   du=du+c(i)*t.^(i-1-gam)*gamma(i)/gamma(i-gam);
end
v=y-u; g=double(genfunc(p)); w=get_vecw(gam,n,g);
for i=1:n, dv(i)=w(1:i)*v(i:-1:1)/h^gam; end
dy=dv+du'; if abs(y(1))<1e-10, dy(1)=0; end

function dy=caputo9(y,t,gam,p)
if gam<0, dy=glfdiff9(y,t,gam,p); return; end
h=t(2)-t(1); t=t(:); n=length(t); y=y(:); q=ceil(gam); 
r=max(p,q); R=sym(fliplr(vander((0:(r-1))'*h)));
c=double(inv(R)*y(1:r)); u=0; du=0;
for i=1:r, u=u+c(i)*t.^(i-1); end
if q<r
   for i=(q+1):p
      du=du+c(i)*t.^(i-1-gam)*gamma(i)/gamma(i-gam);
end, end
v=y(:)-u(:); dv=glfdiff9(v,t,gam,p); dy=dv(:)+du(:);

function [c,y]=caputo_ics(a,na,b,nb,y0,u,t)
na1=ceil(na); q=max(na1); K=length(t); 
p=K+q-1; y0=y0(:); u=u(:); t=t(:); d1=y0./gamma(1:q)';
I1=1:q; I2=(q+1):p; X=zeros(K,p);
for i=1:p, for k=1:length(a)
   if i>na1(k), 
      X(:,i)=X(:,i)+a(k)*t.^(i-1-na(k))*gamma(i)/gamma(i-na(k));
end, end, end
u1=0; for i=1:length(b), u1=u1+b(i)*caputo9(u,t,nb(i),K-1); end
X(1,:)=[]; u2=u1(2:end)-X(:,I1)*d1;
d2=inv(X(:,I2))*u2; c=[d1;d2];
y=0; for i=1:p, y=y+c(i)*t.^(i-1); end

function y = fode_caputo9(a, na, b, nb, y0, u, t, p)
    T=0; dT=0; t=t(:); u=u(:);
    if p>length(y0)
        yy0=caputo_ics(a,na,b,nb,y0,u(1:p),t(1:p)); 
        y0=yy0(1:p).*gamma(1:p)';
    elseif p==length(y0)
        yy0=caputo_ics(a,na,b,nb,y0,u(1:p+1),t(1:p+1)); 
        y0=yy0(1:p+1).*gamma(1:p+1)';
    end
    for i=1:length(y0), T=T+y0(i)/gamma(i)*t.^(i-1); end
    for i=1:length(na), dT=dT+a(i)*caputo9(T,t,na(i),p); end
    u=u-dT; y=fode_sol9(a,na,b,nb,u,t,p)+T';

例10-53 Caputo分数阶线性微分方程
已知线性微分方程
$y^{\prime \prime \prime}(t) + \dfrac{1}{16} _0^{C} \mathscr{D}_t^{2.5} y(t) + \dfrac{4}{5}y^{\prime \prime}(t) + \dfrac{3}{2} y^{\prime}(t) + \dfrac{1}{25} _0^C \mathscr{D}_t^{0.5} y(t) + \dfrac{6}{5} y(t) = \dfrac{172}{125} \cos \dfrac{4t}{5}$
初值: $y(0) = 1, y^{\prime}(0) = 4/5, y^{\prime \prime}(0) = -16/25, 0 \leq t \leq 30$
解析解: $y(t) = \sqrt{2} \sin(4t/5 + \pi/4)$

a = [1, 1/16, 4/5, 3/2, 1/25, 6/5];
na = [3, 2.5, 2, 1, 0.5, 0];
b = 1;
nb = 0;
t = [0: 0.05: 30];
u = 172/125*cos(4*t/5);
y0 = [1, 4/5, -16/25];
y1 = fode_caputo9(a, na, b, nb, y0, u, t, 4);
y = sqrt(2)*sin(4*t/5+pi/4);
max(abs(y-y1)),
plot(t, y, t, y1)

例10-54 非线性微分方程的高精度求解
Caputo定义下的微分方程
$_0^{C} \mathscr{D}_t^{1.455}y(t) = -t^{0.1} \dfrac{E_{1, 1.545}(-t)}{E_{1, 1.445}(-t)} e^{t} y(t)_0^C\mathscr{D}_t^{0.555} y(t) + e^{-2t} - [y^{\prime}(t)]^2$
初值 $y(0) = 1, y^{\prime}(0) = -1,解析解: y(t) = e^{-t}$

function [T,dT,w,d2]=aux_func(t,y0,alpha,p)
an=ceil(alpha); y0=y0(:); q=length(y0); d2=length(alpha); 
m=length(t); g=double(genfunc(p)); 
for i=1:d2, w(:,i)=get_vecw(alpha(i),m,g)'; end
b=y0./gamma(1:q)'; T=0; dT=zeros(m,d2);
for i=1:q, T=T+b(i)*t.^(i-1); end
for i=1:d2
   if an(i)==0, dT(:,i)=T;
   elseif an(i)<q
      for j=(an(i)+1):q
         dT(:,i)=dT(:,i)+(t.^(j-1-alpha(i)))*...
               b(j)*gamma(j)/gamma(j-alpha(i));
end, end, end

function [y,t]=nlfec(fun,alpha,y0,yp,t,p,err)
yp=yp(:); t=t(:); h=t(2)-t(1); m=length(t); ha=h.^alpha;
[T,dT,w,d2]=aux_func(t,y0,alpha,p);
[z,y]=repeat_funs(fun,t,yp,T,d2,alpha,dT,ha,w,m,p);
while norm(z)>err, yp=y; z=zeros(m,1);
   [z,y]=repeat_funs(fun,t,yp,T,d2,alpha,dT,ha,w,m,p);
end
% common code for repetitive computation
function [z,y]=repeat_funs(fun,t,yp,T,d2,alpha,dT,ha,w,m,p)
for i=1:d2, dyp(:,i)=glfdiff9(yp-T,t,alpha(i),p)'+dT(:,i); end
f=fun(t,yp,dyp(:,2:d2))-dyp(:,1); y=yp; z=zeros(m,1);
for i=2:m, ii=(i-1):-1:1;
   z(i)=(f(i)*(ha(1))-w(2:i,1)'*z(ii))/w(1,1); y(i)=z(i)+yp(i);
end

function [y,t]=nlfep(fun,alpha,y0,tn,h,p,err)
m=ceil(tn/h)+1; t=(0:(m-1))'*h; ha=h.^alpha; z=0; 
[T,dT,w,d2]=aux_func(t,y0,alpha,p); 
y=z+T(1); dy=zeros(1,d2-1);
for k=1:m-1
   zp=z(end); yp=zp+T(k+1); y=[y; yp]; z=[z; zp];
   [zc,yc]=repeat_funs(fun,t,y,d2,w,k,z,ha,dT,T);
   while abs(zc-zp)>err
      yp=yc; zp=zc; y(end)=yp; z(end)=zp;
      [zc,yc]=repeat_funs(fun,t,y,d2,w,k,z,ha,dT,T);
end, end
% common code for repetitive computation
function [zc,yc]=repeat_funs(fun,t,y,d2,w,k,z,ha,dT,T)
for j=1:(d2-1)
   dy(j)=w(1:k+1,j+1)'*z((k+1):-1:1)/ha(j+1)+dT(k,j+1);
end, f=fun(t(k+1),y(k+1),dy);
zc=((f-dT(k+1,1))*ha(1)-w(2:k+1,1)'*z(k:-1:1))/w(1,1); yc=zc+T(k+1);

f = @(t, y, Dy)-t.^0.1.*ml_func([1, 1.545], -t).*exp(t)./ ...
    ml_func([1, 1.445], -t).*y.*Dy(:, 1)+exp(-2*t)-Dy(:, 2).^2;
alpha = [1.455, 0.555, 1];
y0 = [1, -1];
tn = 1;
h = 0.01;
err = 1e-8;
[yp1, t] = nlfep(f, alpha, y0, tn, h, 1, err);
tic,
[y2, t] = nlfec(f, alpha, y0, yp1, t, 2, err);
toc,
max(abs(y2-exp(-t)))