bdf_fixed_a.m

BDF法解分数阶微分方程

%对于指定的a,画出真解和近似解,红色的为真解
%解方程y'=D^(1-a)(-y+g(t))  y(0)=0  其中g(t)=t^2+2t^(2-a)/gamma(3-a);
%真解是y=t^2;

clear
begin=clock;
global a
a=0.6    
n_step=100;
h=1/n_step;
y(1)=0;
%%%%%%%%%%%%%%%用一次Eular法计算第二步的值%%%%%%%%%%%%
k=1
 t=(k-1)*h;
 fy(k)=-y(k)+t^2+2*t^(2-a)/gamma(3-a); 
 out=0;
 for j=0:k-1
       out=out+b(j,k-1)*fy(j+1);
 end
 y(k+1)=h^(a-1)/gamma(1+a)*out*h+y(k);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 for k=2:n_step
     t=(k-1)*h;
     fy(k)=-y(k)+t^2+2*t^(2-a)/gamma(3-a); 
     out=0;
     for j=0:k-1
       out=out+b(j,k)*fy(j+1);
     end
     yTemp=y(k);
     fyTemp=-yTemp+(t+h)^2+2*(t+h)^(2-a)/gamma(3-a);
     out1=out+b(k,k)*fyTemp;
     y(k+1)=2/3*(h^(a-1)/gamma(1+a)*out1*h-1/2*y(k-1)+2*y(k));
     while(norm(yTemp-y(k+1))>1e-5)
        yTemp=y(k+1);
        fyTemp=-yTemp+(t+h)^2+2*(t+h)^(2-a)/gamma(3-a);
        out1=out+b(k,k)*fyTemp;
        y(k+1)=2/3*(h^(a-1)/gamma(1+a)*out1*h-1/2*y(k-1)+2*y(k));
     end
 end
x=0:h:1;
plot(x,y);
hold on;
y1=x.^2;
plot(x,y1,'r');
xlabel('t')
ylabel('y')
title('Eular')
text(0.9,y(91),'y_h \rightarrow',...
     'HorizontalAlignment','right')  
text(0.7,0.7^2,'\leftarrow t^2 ',...
     'HorizontalAlignment','left')
 ending=clock;
 timeCost=ending-begin;
 error=norm(y-y1);
 save mydata error timeCost