How are 1st-order differential equations handled?

General
,
1

I’m on Dolfinx version 0.10.0, and I wanted to perform \int_0^x f(x) dx by converting the integral into a 1st-order ODE

\frac{du}{dx} = f(x) \ \ \ \ \forall x\in 0\le x\le L

with f(x) = 1-2x, u(L) = 0; the exact solution is u = x(1-x). When I obtain the weak form by convoluting both sides with a test function v,

\int_0^L v\frac{du}{dx}dx = \int_0^L fv dx

and attempt to solve this problem, the solution goes to infinity everywhere. Here is my code and an image of the results:

from dolfinx.fem.petsc import LinearProblem
import dolfinx as df
import ufl
import pyvista as pv
from mpi4py import MPI
import numpy as np

# --- Mesh creation ---
comm = MPI.COMM_WORLD
nx = 10 # number of points
start_point = 0 # starts at x = 0
L = 1 # ends at x = 1
mesh_1D = df.mesh.create_interval(comm, nx=nx, points = [start_point,L]) #1D line segment
x = ufl.SpatialCoordinate(mesh_1D)[0] # x and y coordinates. Symbolic expressions 


# --- Define Functionspace ---
family = "Lagrange" # The type of functionspace used
degree = 2 # The degree of the polynomial used to interpolate data to
V1D = df.fem.functionspace(mesh_1D, (family, degree))


# --- Define Galerkin of 1st-order ODE ---
# ODE: du/dx = 1-2x; u(1) = 0
# Galerkin: int_0^L [v du/dx] dx = int_0^L [fv] dx
u = ufl.TrialFunction(V1D) 
v = ufl.TestFunction(V1D)
f = 1-2*x

# bilinear term: a(u,v) = int_0^L [v du/dx] dx
dx = ufl.Measure("dx", mesh_1D) #defines differential along mesh
a = ufl.inner(v,u.dx(0)) * dx
# Linear term: L(u) = int_0^L [fv] dx
L_ = ufl.inner(f,v) * dx
# boundary condition: u(0) = 0
bc = df.fem.dirichletbc(value=0.0, dofs = df.fem.locate_dofs_geometrical(V1D, lambda x: np.isclose(x[0],L)), V = V1D) # Assigns dirichlet values to degree of freedom at x = L


# --- Solve ---
ksp_type = "preonly" # Preconditioning Matrix Option
solver_type = "lu" #LU factorization (a direct solver)

problem = LinearProblem(a = a, L = L_, petsc_options_prefix="first_order_ODE", bcs = [bc], petsc_options={"ksp_type": ksp_type, "pc_type": solver_type})
u_sol = problem.solve()

unstable_solution
unstable_solution1600×295 9.3 KB

The solution blows up to infinity at every point, so I’m guessing this variational-form is unstable. Are Galerkins/variational-forms of 1st-order differential equations generally unstable? How should they be handled?

2

You can use Least Squares Finite element methods, i.e.
Reformulate your problem as
\min_u \frac{1}{2}\int_\Omega (u' - f)^2~\mathrm{d}x
which ends up solving
\int_\Omega (u' - f) v' ~\mathrm{d}x

import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import mesh, fem, default_scalar_type
from dolfinx.fem.petsc import LinearProblem

# 1. Create mesh and function space
L = 1.0
nx = 50
domain = mesh.create_interval(MPI.COMM_WORLD, nx, [0.0, L])

# Using degree=2 since the exact solution x(1-x) is quadratic.
# This will capture the solution exactly up to machine precision.
V = fem.functionspace(domain, ("CG", 2))

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(domain)

# 2. Define the source term
u_exact = x[0] * (1.0 - x[0])
f = u_exact.dx(0)

# 3. The Least Squares Weak Form (notice the dx(0) on the test function)
# dx(0) is the derivative with respect to the spatial coordinate x
a = u.dx(0) * v.dx(0) * ufl.dx
L_form = f * v.dx(0) * ufl.dx


# 4. Define Dirichlet Boundary Condition at x = L
def boundary_L(x):
    return np.isclose(x[0], L)


# Locate boundary facets and DOFs
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, boundary_L)
bc_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
compiled_expr = fem.Expression(u_exact, V.element.interpolation_points)
u_bc = fem.Function(V)
u_bc.interpolate(compiled_expr)
bc = fem.dirichletbc(u_bc, bc_dofs)

# 5. Solve the linear problem
problem = LinearProblem(
    a,
    L_form,
    bcs=[bc],
    petsc_options={"ksp_type": "preonly", "pc_type": "lu"},
    petsc_options_prefix="ls_test_",
)
uh = problem.solve()

error_loc = fem.assemble_scalar(fem.form((uh - u_exact) ** 2 * ufl.dx))
error_glob = np.sqrt(domain.comm.allreduce(error_loc, op=MPI.SUM))
print(f"{nx} Error in L2 norm: {error_glob:.3e}")
3

Ah that worked! Thank you! The code comments are also appreciated! So, rather than convoluting with a test function like with a Poisson Equation, you’re minimizing the square of the L2-norm.

Do you know if this method is handled in standard introductory texts on Finite Element Methods, or if this is a more advanced technique?

4

As far as I can tell it is a classical method: