bdf_fixed_a_co2.m

BDF法解分数阶微分方程,分数阶导数定义系数通过fourier法计算

%对于指定的a,画出真解和近似解,红色的为真解
%用BDF法解方程y'=D^(1-a)(-y+g(t))  y(0)=0
%其中g(t)=t^2+2t^(2-a)/gamma(3-a);D^(1-a)的计算系数是fourier系数
%真解是y=t^2;
%误差可通过norm(y-y1)自己计算
clear
begin=clock;
global a k
a=0.6 ;   
n_step=100;
h=1/n_step;
y(1)=0;
%数值积分计算系数Omega2(k)
for k=1:n_step+1
        Omega2(k)=quadl(@Omega1,0,2*pi)/2/pi;%Omega1由Omega1.m定义
end
Omega2=real(Omega2); 
%%%%%%%%%%%%%%%用一次Eular法计算第二步的值%%%%%%%%%%%%
 k=1;
 t=(k-1)*h;
 fy(k)=-y(k)+t^2+2*t^(2-a)/gamma(3-a); 
 D_a=0;
 for j=0:k-1
       D_a=D_a+Omega2(j+1)*fy(k-j);
 end
 y(k+1)=h^(a-1)*D_a*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); 
     D_a=0;
     for j=1:k
         D_a=D_a+Omega2(j+1)*fy(k-j+1);
     end
     yTemp=y(k);
     fyTemp=-yTemp+(t+h)^2+2*(t+h)^(2-a)/gamma(3-a);
     D_a1=D_a+Omega2(1)*fyTemp;
     y(k+1)=2/3*(h^(a-1)*D_a1*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);
        D_a1=D_a+Omega2(1)*fyTemp;
        y(k+1)=2/3*(h^(a-1)*D_a1*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