eular_fixed_a_co2(xiu).m

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

%对于指定的a,画出真解和近似解,红色的为真解
%用Eular法解方程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
global a k
a=0.6;   
n_step=100;
h=1/n_step;
y(1)=0;
%数值积分计算系数Omega2(k)
for k=0:n_step
        Omega2(k+1)=quadl(@Omega1,0,2*pi)/2/pi;%Omega1由Omega1.m定义
end
Omega2=real(Omega2); 
%%%%%%%%%%%%%%

for k=1:n_step
    t=(k-1)*h;
    %修正项系数计算
    w(n+1,1)=2*gamma(1)/gamma(1-a)*n^(-a)-gamma(2)/gamma(2-a)*n^(1-a)
    w(n+1,2)=gamma(2)/gamma(2-a)*n^(1-a)-gamma(1)/gamma(1-a)*n^(-a);
    for j=0:n
        w(n+1,1)=w(n+1,1)-Omega2(n-j+1)*(2-j);
        w(n+1,2)=w(n+1,2)-Omega2(n-j+1)*(j-1);
    end
    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-w(k,))*h+y(k);
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),'\leftarrow y_h',...
     'HorizontalAlignment','left')  
text(0.7,0.7^2,'t^2 \rightarrow',...
     'HorizontalAlignment','right')