分数阶微分方程数值解法

matlab 微分方程数值解 微分变换法

分数阶常微分方程数值解法-微分变换法

一、 题目

考虑方程

$$D^\mu =y^2 +1, \qquad m-1<\mu \leq m,\qquad 0<x<1$$ 初始条件为： $$y^{(k)}(0)=0, \qquad k=0,\cdots,m-1$$

二、 分析

对方程进行微分变换，获得的结果如下：

$$Y(k + \mu \alpha)=\frac{\Gamma(1 + \frac{k}{\alpha})[\sum\limits_{k_1=0}^kY(k_1)Y(k-k_1) + \delta(k)]}{\Gamma(1+\mu+\frac{k}{\alpha})}$$ 对初始条件进行微分变换，得到： $$Y(k) = 0,\qquad k=0,1,\cdots,\mu\alpha-1$$ 现在可以利用由黎曼-刘维尔公式定义的分数阶微分扩展连续函数 $f(x)$ 的分数幂级数上的形式： $$f(x)=\sum\limits_{k=0}^\infty F(K)(x-x_0)^{\frac{k}{\alpha}}$$ 其中 $F(k)$ 是微分变换下相应的系数，所以只需要计算出相应的 $F(k)$ 便可以利用以上的公式求出函数在某点处的近似值，编写程序代码如下计算

三、代码

我们取$\quad N=30, \quad \alpha = 2,\quad \mu = 0.5\&1.5,\quad t= 0.0,0.1,\cdots,1.0$
Matlab 代码如下:

function [] = main(N,alpha,miu )
% 分数阶常微分方程数值方法之微分变换法 方程为 D(u)(y(x))=y^2 +1；
% 输入的参数为 N,alpha miu 三种其中 N 是 大于 0的整数，alpha 为大于0的整数，miu 为大于零的小数；
%% 首先对数据进行初始化；startp 代表微分方程的定义域的左端，endp 代表微分方程定义域的右端点；ydao 代表 存储 F(k)的数值；N 代表级数的个数；
[startp,endp] = deal(0,1);
ydao = zeros(N,1);
xarr = startp:0.1:endp;
miualpha= miu*alpha;
%% 通过循环计算函数的导数的值并存入向量中；
for t = miualpha:1:N+3;
    k = t-miualpha;
    upgamma1 = gamma(1+k/alpha);
    upright=0;
    if k==0;
        upright = upright +1;
    end
    for k1=0:k
        upright = upright + ydao(k1+1)*ydao(k-k1+1);
    end    
    ydao(t+1)=upgamma1*upright/gamma(1+miu+k/alpha);
end
%disp(ydao);
%% 通过循环计算函数在定义域的近似数值；
leng = length(xarr);
y=zeros(leng,1);
for m=1:1:leng;
    for k = 1:N+1;
        y(m) = y(m)+ydao(k)*(xarr(m)-startp)^((k-1)/alpha);
    end
end
disp(y);
end

四、结果

调用函数为：

     main(30,2,1.5)

此时取$\alpha =2, \mu =1.5$ 获得的结果为：

    0
    0.0238
    0.0673
    0.1239
    0.1914
    0.2689
    0.3562
    0.4540
    0.5630
    0.6851
    0.8225

画出的图像为：

分数阶微分方程数值解法-多步微分变换法

一、题目

给定方程为

$$D^{\mu}y(x)=y^2+1\qquad m-1<\mu \leq m,\qquad 0<x<5$$ 初始条件为： $$y^{(0)}(k)=0;\qquad k=0,\dots,m-1$$

二、分析

对方程进行微分变换，我们得到：

$$Y(k + \mu \alpha)=\frac{\Gamma(1 + \frac{k}{\alpha})[\sum\limits_{k_1=0}^kY(k_1)Y(k-k_1) + \delta(k)]}{\Gamma(1+\mu+\frac{k}{\alpha})}$$ 对初始条件条件进行微分变换，得到： $$Y^{(0)}(0)=0\quad k=0,\dots,\mu\alpha-1$$ 现在可以利用由黎曼-刘维尔公式定义的分数阶微分扩展连续函数 $f(x)$ 的分数幂级数上的形式： $$f(x)=\sum\limits_{k=0}^\infty F(K)(x-x_0)^{\frac{k}{\alpha}}$$ 其中 $F(k)$ 是微分变换下相应的系数，所以只需要计算出相应的 $F(k)$ 便可以利用以上的公式求出函数在某点处的近似值.
当获得近似值后，可以对定义域中由分成五个小段逐个求得在 $x=0,1,2,3,4$ 的数值，使用微分变换法便可求得近似的数值解。
对函数从 $0-5$ 分成5段，依次进行微分变换，求得方程的数值解；
代码如下:

三、代码

主函数的代码：

function [ ] = fenbu_weifen(N,)
% 本函数使用多步微分变换法求解分数阶微分方程求解 方程为 D(u)(y)=3*y+1;
%   所需的参数为N，级数的级数，alpha大于零的整数 miu 是分数；
miu = 1;
alpha=2;
startydao=zeros(miu*alpha+1,1);
truy =[];
xarr=[];
for t=1:5;
    [startp,endp]=deal((t-1)*0.2,t*0.2);
    x=startp:0.1:endp;
    [ydaoend,y] = DeriTrans(N,startp,endp,startydao );
    startydao = ydaoend;
    disp(y);
    truy=[truy ;y];
    xarr=[xarr x];
    disp(ydaoend);
end
size(xarr)
size(truy)
plot(xarr,truy);
end

所用的 DeriTransh 子函数代码如下：

function [ydaoend,y] = DeriTrans(N,startp,endp,startydao )
% 分数阶常微分方程数值方法之微分变换法 方程为 D(u)(y(x))=y^2 +1；
%   输入的参数为 N,alpha miu 三种其中 N 是 大于 0的整数，alpha 为大于0的整数，miu 为大于零的小数；
%   函数的返回值为ydao 即函数在终点处的各阶导数的值；y为函数的数值解；
%% 首先对数据进行初始化；startp 代表微分方程的定义域的左端，endp 代表微分方程定义域的右端点；ydao 代表 存储 F(k)的数值；N 代表级数的个数；
alpha=2;
miu=1;
ydao = zeros(N,1);
ydao(1:miu*alpha+1)=startydao;
xarr = startp:0.1:endp;
miualpha= miu*alpha;
%% 通过循环计算函数的导数的值并存入向量中；
for t = miualpha:1:N+3;
    k = t-miualpha;
    upgamma1 = gamma(1+k/alpha);
    upright=0;
    if k==0;
        upright = upright +1;
    else
        for k1=0:k
            upright = upright + ydao(k1+1)*ydao(k-k1+1);
        end
    end
    ydao(t+1)=upgamma1*upright/gamma(1+miu+k/alpha);
end
%disp(ydao);
%% 通过循环计算函数在定义域的近似数值；
leng = length(xarr);
y=zeros(leng,1);
for m=1:1:leng;
    for k = 1:N+1;
        y(m) = y(m)+ydao(k)*(xarr(m)-startp)^((k-1)/alpha);
    end
end
%disp(ydao);
%dos('explorer http://mathnews.')
%% 计算导数在右端点处的数值，k-1代表求导的阶数；
ydaoend=zeros(miualpha,1);
val1=0;
val2=0;
for n=1:N    
    temp=(endp-startp)^((n-1)/alpha)*ydao(n);
    val1 = val1 + temp;
    temp = (endp-startp)^((n-1)/alpha-1)*ydao(n);
    val2 = val2 + temp*(n-1)/alpha;    
end
ydaoend(1) = val1;
ydaoend(3) = val2;
% syms x;
% expr=0;
% for n=1:N;
%     expr=expr+ydao(n)*(x-startp)^((n-1)/alpha);
% end
% 
% for k=0:miu  
%     res=diff(expr,x,k);
%     ydaoend(k*alpha+1)=subs(res,x,endp);
%    
% end
end

四、结果

Matlab 输出的结果如下:
调用函数：

fenbu_weifen(30,2)

输出结果为：

0
0.1003
0.2027

0.2027
     0
1.0411

0.2027
0.3051
0.4140

0.4140
     0
1.1303

0.4140
0.5186
0.6347

0.6347
     0
1.2315

0.6347
0.7417
0.8659

0.8659
     0
1.3469

0.8659
0.9755
1.1088

1.1088
     0
1.4796

画出的图像为：

在代码中 ，我们取 $\mu$ 的数值为2；$\alpha$ 的数值为2 这样获得的图像如上，与真实的值相差比较小。

分数阶微分方程数值解法-方程组解法

一、题目

考虑分数阶微分代数系统：

$$\begin{cases} D^{\alpha}_*x(t)-ty'(t)+t^2z'(t)+x(t)-(1+t)y(t)+(t^2+2t)z(t)=0\\ D^{\alpha}_*y(t)-tz'(t)-y(t)+(t-1)z(t)=0\\ z(t) = sint,0<\alpha\leq1 \end{cases}$$
满足初始条件： $$x(0)=y(0)=1,z(0)=0$$

二、分析

对方程进行微分变换，获得方程如下：

$$X(0)=1,Y(0)=1,Z(0)=0\\ \frac{\Gamma(\alpha (k+1)+1)}{\Gamma(\alpha k+1)}X(k+1)=(k+1)Y(k)-(k+1)Z(k-1)-X(k)+Y(K-1)-Z(k-2)\\ \frac{\Gamma(\alpha (k+1)+1)}{\Gamma(\alpha k+1)}Y(k+1)=(k+1)Z(k)+Y(k)-Z(k-1)\\ Z(2k)=0,Z(2k+1)=(-1)^k/(2k+1)!$$ 其中 $$Y(-1)=Z(-1)=Z(-2)=0$$

三、代码

四、结果