分数阶微分方程数值实验MATLAB编码

分数阶微分方程数值实验实验题目：考虑分数阶扩散微分方程),(),()(),(txqxtxuxdttxu(1.1)这里的6)2.2()(1xxd，3)1(),(xextxqt，其中初值为30,xxu，边值tetutu),1(0),0(，其真解为3),(xetxut，计算其数值解。实验算法：1.将空间区间[0,1]等距剖分成N段，1N个节点为12101Nxxx；将时间区间]1,0[等距剖分成M段，1M个节点为1...0121Mttt。2．将方程组(1.1)中的xtxu,用有限GrunwardLetnikov算子离散，即210,210)1(),(jikkjkikjiGLFkikiughukhtxuD其中)1()1()1()1()1(,kkkgkkki1,2,...,1N，1,,2,1Mj其中是分数阶。再对21jkiu利用中心差分2121jkijkijkiuuu进行离散，则得到xtxu,的离散格式)(2110,210,jkijkiikkjikkuughughki将方程(1.1)中的ttxu,利用jijiuu-1进行离散，其中为时间步长方程(1.1)的离散格式为2110,12-jijkijkiikkijijiquughduu即21-0,10,122jijkiikkijijkiikkijiqughduughdu（1.2）i1,2,...,1N，1,,2,1Mj等价于下面的矩阵形式211)(jjjFUAIUAI（1.3）其中,01,11,1,11,02,12,22,01,11000gagagagagagagagagaANNNNNNN，这里的hdaii2，21212121121jNjjjqqqF要求方程(1.1)的数值解，即求系统3.1。程序代码：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;%空间步数N=M;%时间步数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;%右端函数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);%构造矩阵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;%数值解enderror=abs(u-ue);figure[X,Y]=meshgrid(x,t);mesh(X,Y,error');该图是方程1.1对于8.1的数值解