The numerical result is not cosistant with the exact solution.
I’m a very new user, can someone help me to check the following codes.
The problem is
-u’’ + u = x, 0 < x < 1,
u(0) = 0, u(1) = 0,
the exact solution is u(x) = x - sinh(x)/sinh(1).
import numpy as np
import ufl
from dolfinx import fem, io, mesh, plot
from ufl import ds, dx, grad, inner
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
# +
msh = mesh.create_interval(comm=MPI.COMM_WORLD, points=(0.0, 1.0), nx=50)
V = fem.FunctionSpace(msh, ("Lagrange", 1))
facets = mesh.locate_entities_boundary(msh, dim=0,
marker=lambda x: np.logical_or(np.isclose(x[0], 0.0),
np.isclose(x[0], 1.0)))
dofs = fem.locate_dofs_topological(V=V, entity_dim=0, entities=facets)
bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)
# +
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = x[0]
a = inner(grad(u), grad(v)) * dx
L = inner(f, v) * dx
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
# -
# +
with io.XDMFFile(msh.comm, "out_poisson/poisson.xdmf", "w") as file:
file.write_mesh(msh)
file.write_function(uh)
# try:
# import pyvista
# cells, types, x = plot.create_vtk_mesh(V)
# grid = pyvista.UnstructuredGrid(cells, types, x)
# grid.point_data["u"] = uh.x.array.real
# grid.set_active_scalars("u")
# plotter = pyvista.Plotter()
# plotter.add_mesh(grid, show_edges=True)
# warped = grid.warp_by_scalar()
# plotter.add_mesh(warped)
# plotter.show()
# except ModuleNotFoundError:
# print("'pyvista' is required to visualise the solution")
# print("Install 'pyvista' with pip: 'python3 -m pip install pyvista'")
# -
import matplotlib.pyplot as plt
cells, types, x = plot.create_vtk_mesh(V)
plt.figure(1,dpi=300)
exact = x[:,0] - np.sinh ( x[:,0] ) / np.sinh ( 1.0 )
plt.plot ( x[:,0], uh.x.array.real)
plt.plot ( x[:,0], exact)
filename = 'bvp_02_solution.png'
plt.savefig ( filename )
print ( ' Graphics saved as "%s"' % ( filename ) )
plt.show()
plt.close ( )
