主函数比较

% 根据文献自己编写

q = 0.5; %分数阶阶数

fdefun = @(t,y) [2/gamma(3-q)*t^(2-q)-1/gamma(2-q)*t^(1-q)-y+t^2-t]; %一元微分方程

y0 =0;%初值y0为列向量

h = 2^(-6);%步长

tspan = [0,2];

[t,y] = fdewfc(q,fdefun,tspan,y0,h);

figure(1)

plot(t,y(1,1:end)) ;

xlabel('t'); ylabel('y(t)');

hold on

plot(t,t.^2-t+0.5,'r-.')%解析解

%plot(t,-2./(t+1),'r')

%%%%%%与意大利 Roberto Garrappa, University of Bari, Italy 结果比较

t0=tspan(1); tfinal = tspan(2);

[t, y_fde12] = fde12(q,fdefun,t0,tfinal,y0,h);

plot(t,y_fde12+1,'k.')

Jfdefun = @(t,y)[-1];%Jacobi阵 此处退化为导数

[t, y_flmm2] = flmm2(q,fdefun,Jfdefun,t0,tfinal,y0,h);

plot(t,y_flmm2+1.5,'m--')

legend('本方法','解析解','Roberto Garrappa方法1','Roberto Garrappa方法2','location','northwest') ;

title('预估-校正ABM法与其他方法比较');

子函数

function [t,y] = fdewfc(q,fdefun,tspan,y0,h)

% 根据文献自己编写

%q = 0.5;

%fdefun = @(t,y) [2/gamma(3-q)*t^(2-q)-1/gamma(2-q)*t^(1-q)-y+t^2-t]; %一元微分方程

%初值y0为列向量

%步长h = 2^(-6) ;

t0 = tspan(1) ; tfinal = tspan(2);

t = t0:h:tfinal;% 时间向量

itn = length(t)-1;%ceil((tfinal-t0)/h);

rd = length(y0); %向量维数

y = zeros(rd,itn+1); y(:,1) = y0; %赋初值 按列向量排

gq1 = 1./gamma(q); %预估用常数

gq2 = (h.^q)./gamma(q+2);%校正用常数

%%%%%%%%%%主循环

for n = 0:itn-1

     % 生成B b_{j,n+1} = h^{\alpha}/\alpha((n+1-j)^{\alpha}-(n-j)^{\alpha})

   % 共n+1项，这里用 q 表示 \alpha

   jv=0:n; B = h^q/q*((n+1-jv).^q-(n-jv).^q);% n+1个元素；必须注意：每步迭代更新B

   % 生成A a_{0,n+1} = n^{\alpha+1}-(n-\alpha)(n+1)^{\alpha}, j=0,

   %       a_{j,n+1} =

   %       (n-j+2)^{\alpha+1}+(n-j)^{\alpha+1}-2(n-j+1)^{\alpha+1},1\leqslant

   %       j\leqslant n,

   %       a_{n+1,n+1}=1,j=n+1.

   % 共n+1项，这里用 q 表示 \alpha

   jv = 1:n; A = [n^(q+1)-(n-q)*(n+1)^q,(n-jv+2).^(q+1)+(n-jv).^(q+1)-2*(n-jv+1).^(q+1)]; % n+1个元素；每步迭代更新 A

    %%%% 预估 yp

    fdefunp = y0; %B = fdeB(q,h,n);%必须注意：每步迭代更新B

    for j = 1:n+1

       fdefunp = fdefunp + gq1.*B(:,j).*fdefun(t(j),y(:,j));

    end

   %%%%后一个子循环求和 %第1列和第3列合并

   fdefunsumn = 0; %fdefunsumn 应为与y0同型的零列向量，但此处写为0，无妨

   %A = fdeA(q,n); %必须注意：每步迭代更新A

   for j = 1:n+1

       fdefunsumn = fdefunsumn + A(:,j).*fdefun(t(j),y(:,j));%0~~n项

   end

   %%%% 校正

   y(:,n+2) = y0 + gq2.*(fdefun(t(n+2),fdefunp)+ fdefunsumn); %n

end