|
主函数比较
% 根据文献自己编写
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
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-19 14:22
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社