分数阶微分方程数值实验
实验题目：
考虑分数阶扩散微分方程
),()
,()(),(t x q x t x u x d t t x u +∂∂=∂∂α
α (1.1) 这里的6
)2.2()(1
+Γ=αx x d ，3)1(),(x e x t x q t -+-=，其中初值为()30,x x u =，边值
⎩⎨⎧==-t
e t u t u ),1(0
),0(，其真解为3),(x e t x u t -=，计算其数值解。
实验算法：
1.将空间区间[0,1]等距剖分成N 段，1+N 个节点为
12101N x x x +=<<
<= ；将时间区间]1,0[等距剖分成M 段，1+M 个节点为
1...0121=<<<=+M t t t 。
2．将方程组(1.1)中的()ααx
t x u ∂∂,用有限G runward Letnikov -算子离散，即 2
10,210)1(),(+=-+=---∑∑=⎪⎪⎭
⎫ ⎝⎛-=j i k k j k
i
k j i GL F
k i k i u g h u k h t x u D ααα
αα
其中)
1()1()
1()1()1(,+Γ+-Γ+Γ-=⎪
⎪⎭
⎫ ⎝⎛-=k k k g k
k k αααα
i 1,2,...,1N =+，1,,2,1+=M j
其中α 是分数阶。 再对21
+-j k i u
利用中心差分212
1
+--+-+=j k i j k i j k i u u u 进行离散，则得到()α
αx
t x u ∂∂,的离散格式)(2110
,2
1
,+--=-+
=-+=∑∑-j k i j k i i
k k j i
k k u u g h u g
h
k
i ααα
α
将方程(1.1)中的()t
t x u ∂∂,利用τj i j i u u -1+进行离散，其中τ为时间步长 方程(1.1)的离散格式为
()2
1
10
,
12-++--=+++=∑j i
j k
i j k
i i
k k i j
i j i q
u
u
g h
d u u αατ
即ττταα
αα21
-0
,10
,1
22+=+-=+++=-∑∑j i j k i i
k k i j
i j k
i i
k k i j i
q u g h
d u u
g h
d u （1.2）
i 1,2,...,1N =+，1,,2,1+=M j
等价于下面的矩阵形式
()2
1
1
)(++++=-j j j F U A I U
A I τ （1.3）
其中⎪⎪
⎪⎪⎪
⎪
⎭
⎫
⎝
⎛=+-++++αα
α
αααα
αα,01,11,1,11,02,12,22,01,1100
0g a g a g a g a g a g a g a g a g a A N N N N N N N
，这里的α
τh d a i i 2=， ⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛=+++++212121
2
1
121j N j j j q q q F
要求方程(1.1)的数值解，即求系统()3.1。
程序代码：
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 );
