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

分数阶微分方程数值实验实验题目：考虑分数阶扩散微分方程这里的，其中初值为，边值，其真解为，计算其数值解。实验算法：1.将空间区间等距剖分成段，个节点为 ；将时间区间等距剖分成段，个节点为。2将方程组中的用有限算子离散，即其中，其中 是分数阶。再对利用中心差分进行离散，则得到的离散格式将方程中的利用进行离散，其中为时间步长方程的离散格式为即 (1.2)，等价于下面的矩阵形式(1.3) 其中，这里的，要求方程的数值解，即求系统。程序代码：function gg_alph = g( M,alph )gg_alph = zeros( M+1,1 );gg_alph(1,1) = 1;for i = 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;for k=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);for i = 1:M-1D( i,1 ) = d( x( i+1 ) );enda = tau * D / ( 2 * halph );gg = g( M,alph );for i = 1:M-1for k = 1:N-1if k = i-1A( i,k ) = a( i,1 ) * gg( i-k+2,1 );elseif k = iA( i,k ) = a( i,1 ) * gg( 2,1 );elseif k = i+1A( i,k ) = a( i,1 ) * gg( 1,1 );else A( i,k ) = 0;endendendfor k = 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);figureX,Y=meshgrid(x,t);mesh(X,Y,error);该图是方程对于的数值解