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

1、精选优质文档-倾情为你奉上精选优质文档-倾情为你奉上专心-专注-专业专心-专注-专业精选优质文档-倾情为你奉上专心-专注-专业分数阶微分方程数值实验实验题目：考虑分数阶扩散微分方程 这里的，其中初值为，边值，其真解为，计算其数值解。实验算法：1.将空间区间等距剖分成段，个节点为 ；将时间区间等距剖分成段，个节点为。2将方程组中的用有限算子离散，即其中，其中 是分数阶。再对利用中心差分进行离散，则得到的离散格式将方程中的利用进行离散，其中为时间步长方程的离散格式为即 （1.2），等价于下面的矩阵形式 （1.3） 其中，这里的， 要求方程的数值解，即求系统。程序代码：function gg_alp

2、h = 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 ) ); 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,

3、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 );

4、 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 * halph ); 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

5、 ); 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