Anomalous diffusion equation

Hello everyone, I’m solving the anomalous diffusion equation. I consider a linear one-dimensional case. The problem is that at the initial time the solution coincides with the analytical one, but not at other times. Who knows where the error is? The code uses the discrete derivative
Grunwald-Letnikov order alpha as a difference analogue of the fractional product
water Caputo, taking into account the memory effect
Снимок экрана от 2022-12-18 14-22-09

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import scipy.special

alpha = 1 #Derivative exponent

T = 1 #Maximum time (up to which we calculate u)
M = 10 #Number of intervals in the time grid
dt = T/M #Time step

u0 = Expression('-1./2*sin(3*pi*x[0]) + 3./2*sin(pi*x[0])', degree=3) #Initial condition

mesh = UnitIntervalMesh(20) #Set up a spatial grid
V = FunctionSpace(mesh, 'CG', 4) #We define the space of functions

#Set the boundary conditions
def boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(V, u0, boundary)

#We introduce the functions u and v in order to determine a(u,v) and L(v)
v = TestFunction(V)

ut = [interpolate(Constant(0.0), V) for j in range(M+1)] #The value of the function u at different times
ut[0] = interpolate(u0, V) #initial #Initial condition

#The functions q(u) and f(u) from the problem statement
def q(uk):
	return 1

def f(uk):
	return 0


tol = 1.0E-5 #Accuracy at which we stop iterations
maxiter = 50 #The maximum allowable number of iterations (the condition on the number of iterations can be removed in principle)

t = dt #We start the calculations from the first moment in time

while t <= T: #Time cycle
	print("t=",t) #Can be removed
	j = int(t//dt) #Determine the node number of the time grid
	u_k = interpolate(Constant(0.0), V) #Setting the value of u at the zero iteration
	L = (ut[0] + (dt**alpha)*f(u_k) - sum( (-1)**k * scipy.special.binom(alpha, k) * (ut[j-k]-ut[0]) for k in range(1, j)) )*v*dx
	eps = 1.0 #Error ||u_k-u_k-1||
	itern = 0 #Iteration number
	while eps > tol and itern < maxiter: #Iterations
		itern += 1 #Iteration number
		u = TrialFunction(V)
		a = u*v*dx + (dt**alpha)*q(u_k)*inner(nabla_grad(u), nabla_grad(v))*dx
		u = Function(V)
		solve(a == L, u, bc)
		#Calculating differences between alternative iterations and observations
		diff = u.vector() - u_k.vector()
		eps = np.linalg.norm(diff, ord=np.Inf)
		u_k.assign(u) #Update the value of u
	ut[j].assign(u) #We write the value of u to the array of values ​​of the function u at different points in time
	t += dt #move on to the next point in time
    

print("complete")

#Building a graph
plt.title('q = 1, f = 0, alpha = 1')
plt.grid(True)
#plot(ut[0])
plot(ut[1])
#plot(ut[M-3])
#plot(ut[M-2])
#plot(ut[M-1])
#plt.legend(['ut[M-1]'])
#plt.show()

def u_analytic(x):
   return np.sin(3*pi*x)*(-1/2)*np.exp(-(9*pi**2)*0.1)+ np.sin(pi*x)*3/2*np.exp(-(pi**2)*0.1) #t=0.1
x = np.linspace(0, 1, 80)
plt.plot(x, u_analytic(x), color='red')

#plt.legend(['ut[M-1]'])
plt.show()