分数阶常微分方程 微分变换法

微分方程 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

四、结果