Why does the graph converge at different points in time to the same point? When “xi” is responsible for the spatial coordinate. The graph should with the right border tends to 0.
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import scipy.special
import math
import os
plt.figure(figsize=(10, 5), dpi=80)
set_log_level(LogLevel.ERROR)
alpha = 1
T = 1
M = 100
dt = T/M
nx = 65
t = -1
k0 = 1
A = 1
sigma = 2
t_f = 0
L_f = np.sqrt(2*k0*(A**sigma)*(sigma+2)/sigma)
print('L_f =', L_f)
n = -1./(0.5*sigma)
xi = Expression('x[0]/(pow(k0, 1/2)*pow(A, sigma/2)*pow(-t, (1+n*sigma)/2))', degree=10, sigma=sigma, n=n, t=t, k0=k0, A=A)
u0 = Expression('A*pow(-t,n)*pow((1.-xi/(pow(2*k0*(pow(A,sigma))*(sigma+2)/sigma, 0.5))),(2./sigma))', degree=10, sigma=sigma, t=t, xi=xi, n=n, A=A, k0=k0) #точное решение
mesh = IntervalMesh(nx, 0, L_f)
V = FunctionSpace(mesh, 'CG', 2)
def b1(x, on_boundary):
if x[0] <= L_f:
return near(x[0], 0)
if x[0] > L_f:
return 0
def b2(x, on_boundary):
if x[0] <= L_f:
return near(x[0], L_f)
if x[0] > L_f:
return 0
bc1 = DirichletBC(V, u0, b1)
bc2 = DirichletBC(V, 0, b2)
v = TestFunction(V)
u = Function(V)
def q(uk):
return uk**sigma
ut = [interpolate(u0, V) for j in range(M)]
for k in range(1,M):
t = t + dt
xi_values = xi.compute_vertex_values(mesh)
xis = xi(1,66)
xi.t = t
u0.t = t
u0.xi = xi
F = u*v*dx + (dt**alpha)*dot(q(ut[k-1])*grad(u), grad(v))*dx - (ut[0] - sum( (-1)**j * scipy.special.binom(alpha, j) * (ut[k-j]-ut[0]) for j in range(1, k)) )*v*dx
J = derivative(F, u)
problem = NonlinearVariationalProblem(F, u, [bc1,bc2], J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-5
prm['newton_solver']['relative_tolerance'] = 1E-4
prm['newton_solver']['maximum_iterations'] = 50000
prm['newton_solver']['relaxation_parameter'] = 1.0
solver.solve()
ut[k].assign(u)
error_L2 = errornorm(u0, u, 'L2')
error_H1 = errornorm(u0, u, 'H1')
plot(u, label='t=%.2f'%(t))
print(' ')
print (f'alpha={alpha}, \u0394t={dt}, \u0394x={1/nx}, M={M}, error_H1={error_H1}, error_L2={error_L2}')
plt.title(f'alpha={alpha}, t_f={t_f}, \u0394x={1/nx}')
plt.grid(True)
plt.ylim(-0.5, 25)
plt.show()
print('complete')