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

分数阶微分方程数值实验 实验题目： 考虑分数阶扩散微分方程 这里的，，其中初值为，边值，其真解为，计算其数值解。 实验算法： 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:M gg_alph( i+1,1 ) = gamma( i-alph ) / ( gamma( -alph ) * gamma( i+1 ) ); end End 主程序 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+1 ue(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); %构造矩阵A A=zeros(M-1,M-1); for i = 1:M-1 D( i,1 ) = d( x( i+1 ) ); end a = tau * D / ( 2 * h^alph ); gg = g( M,alph ); for i = 1:M-1 for k = 1:N-1 if k = i-1 A( i,k ) = a( i,1 ) * gg( i-k+2,1 ); elseif k == i A( i,k ) = a( i,1 ) * gg( 2,1 ); elseif k == i+1 A( i,k ) = a( i,1 ) * gg( 1,1 ); else A( i,k ) = 0; end end end for k = 1:N b = ( 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;%数值解 end error = abs(u-ue); figure [X,Y]=meshgrid(x,t); mesh(X,Y,error); 该图是方程对于的数值解