Hello! Question: Why do charts at different times end at the same point? Below I will attach my schedule, and what schedule should be obtained
T = 1 #Maximum time (up to which we calculate u)
M = 5 #Number of intervals in the time grid
dt = T/M #Time step
nx = 65
t = 0.0
omega = 2
beta = 3 # beta = omega + 1 S-mode, beta > omega + 1 LS-mode, 1 < beta < omega + 1 HS-mode
m = 0.5*(beta-omega-1)/(beta-1)
print(f'm={m}')
L_f = (2*pi*(omega+1)**0.5)/omega
print(L_f)
tf = dt #exacerbation time
xsi = Expression('x[0]/(pow(1-t/tf,m))', degree=10, m=m, t=t, tf=tf)
u0 = Expression('(pow(1-t/tf,-1/(beta-1)))*pow(pow(cos((pi*xsi)/L_f), 2)*(2*(omega+1)/(omega*(omega+2))), 1/omega)', degree=10, t=t, L_f=L_f, omega=omega, beta=beta, m=m, xsi=xsi, tf=tf) #Initial condition
mesh = IntervalMesh(nx, -L_f/2, L_f/2) #Set up a spatial grid
V = FunctionSpace(mesh, 'CG', 2) #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)
u = TrialFunction(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 uk**omega
def f(uk):
return uk**beta
f_ = 1
mp.dps = 15
t_f = 0
#Time cycle
for j in range(1, M + 1):
t = dt * j
tf = mp.mpf(tf) + mp.mpf(dt)
xsi.t = t
xsi.tf = tf
u0.t = t
u0.tf = tf
#f.t = t
u_k = interpolate(u0, V) #Setting the value of u at the zero iteration
L = (ut[0] + (dt**alpha)*f_*f(u_k) - sum( (-1)**k * scipy.special.binom(alpha, k) * (ut[j-k]-ut[0]) for k in range(1, j)) )*v*dx
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)
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
error_L2 = errornorm(u0, u, 'L2')
error_H1 = errornorm(u0, u, 'H1')
plot(ut[j-1])
#Obtaining a local representation of the vector u
u_local = u.vector().get_local()
# Infinity test
has_nan = np.isnan(u_local)
#Infinity test
if has_nan.any():
print("Vector contains infinite values.")
print(f'tf={tf}')
break
plt.show()
My schedule
What should be
