eular_vary_a.m

Eular法解分数阶微分方程,分数阶导数定义系数通过fourier法计算,阶数a变化

%当a变化时(取a=0.1,0.2,0.3,...,1),画出误差图
%解方程y'=D^(1-a)(-y+g(t))  y(0)=0  其中g(t)=t^2+2t^(2-a)/gamma(3-a);
%真解是y=t^2;
clear
global a
n_step=100;
step=1;
while n_step<=400  
     for a1=1:10 %a变化
%%%%%%%%%%%%%%%%对不同的a解方程%%%%%%%%%%%%%%%%%%%%%%%%
        a=a1*0.1;
        h=1/n_step;               
        y(1)=0;                  
        for k=1:n_step             
             t=k*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);
        end
        x=0:h:1;
        y1=x.^2;
        error(a1,step)=norm(y1-y);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     end
     n_step=n_step*2;%剖分变成两倍
     step=step+1;
end
%%%画图
a1=1:10;
k=1;
%subplot(1,step,k)
plot(a1,error(:,k));
xlabel('a')
ylabel('error')
title('Eular')
text(a1(2*k),error(2*k,k),'\leftarrow h=0.01',...
     'HorizontalAlignment','left')
 hold on
k=2;
%subplot(1,step,k)
plot(a1,error(:,k),'r');
text(a1(2*k),error(2*k,k),'\leftarrow h=0.005',...
     'HorizontalAlignment','left')
 hold on
k=3;
% subplot(1,step,k)
plot(a1,error(:,k),'g');
text(a1(2*k),error(2*k,k),'\leftarrow h=0.0025',...
     'HorizontalAlignment','left')
%%%%%%%%%误差比值
 error(:,3)./error(:,2)
 error(:,2)./error(:,1)