分數(shù)階微分方程數(shù)值實驗MATLAB編碼

分數(shù)階微分方程數(shù)值實驗實驗題目：考慮分數(shù)階擴散微分方程這里的，，其中初值為，邊值，其真解為，計算其數(shù)值解。實驗算法：1.將空間區(qū)間等距剖分成段，個節(jié)點為；將時間區(qū)間等距剖分成段，個節(jié)點為。2．將方程組中的用有限算子離散，即其中，其中是分數(shù)階。再對利用中心差分進行離散，那么得到的離散格式將方程中的利用進行離散，其中為時間步長方程的離散格式為即〔1.2〕，等價于下面的矩陣形式〔1.3〕其中，這里的，要求方程的數(shù)值解，即求系統(tǒng)。程序代碼：functiongg_alph=g(M,alph)gg_alph=zeros(M+1,1);gg_alph(1,1)=1;fori=1:Mgg_alph(i+1,1)=gamma(i-alph)/(gamma(-alph)*gamma(i+1));endEnd主程序T=1;M=100;%空間步數(shù)N=M;%時間步數(shù)h=1/M;%空間步長tau=T/N;%時間步長x=0:h:1;t=0:tau:T;alph=1.8;ue=zeros(M+1,N+1);u=ue;D=zeros(M-1,1);a=D;f=@(x,t)-(1+x).*exp(-t).*x.^3;%右端函數(shù)initial_condation=@(x)x.^3;left_boundary=@(t)0;right_boundary=@(t)exp(-t);exact=@(x,t)exp(-t).*x.^3;d=@(x)gamma(2.2)*x.^2.8/6;fork=1:N+1ue(1:end,k)=exact(x(1:end),t(k));%真解end%問題初邊值條件u(1:end,1)=initial_condation(x);u(1,1:end)=left_boundary(t);u(end,1:end)=right_boundary(t);%構(gòu)造矩陣AA=zeros(M-1,M-1);fori=1:M-1D(i,1)=d(x(i+1));enda=tau*D/(2*h^alph);gg=g(M,alph);fori=1:M-1fork=1:N-1ifk<=i-1A(i,k)=a(i,1)*gg(i-k+2,1);elseifk==iA(i,k)=a(i,1)*gg(2,1);elseifk==i+1A(i,k)=a(i,1)*gg(1,1);elseA(i,k)=0;endendendfork=1:Nb=(eye(M-1)+A)*u(2:end-1,k)+tau*f(x(2:end-1),t(k)+tau/2)'+...a.*(gg(3:end)*(u(1,k+1)+u(1,k)).*ones(M-1,1));b(end,1)=b(end,1)+a(end,1)*gg(1)*(u(end,k+1)+u(end,k));u(2:end-1,k+1)=(eye(M-1)-A)\b;%數(shù)