Galerkin Method¶
The Galerkin method is a popular way of solving (partial) differential equations by discretizing them and solving the resulting linear system. Here, we will do so using a probabilistic linear solver. We consider the Dirichlet problem for the Poisson equation given by
on a pre-generated mesh.
[1]:
# Make inline plots vector graphics instead of raster graphics
%matplotlib inline
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('pdf', 'svg')
# Plotting
import matplotlib.pyplot as plt
plt.style.use('../probnum.mplstyle')
Discretization Mesh¶
We used FEniCS to compute a mesh on a circular domain. The mesh results from a choice of basis functions comprising a discrete subspace of the function space \(V\) of solutions \(u\).
[2]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
# Setup
data_path = "../data/PDE_discretization/"
# Load mesh data
mesh_resolution_coarse = 6
mesh_resolution_fine = 128
mesh_xy = np.load(file=data_path + "mesh_xy_res{}.npy".format(mesh_resolution_coarse))
mesh_cells = np.load(file=data_path + "mesh_cells_res{}.npy".format(mesh_resolution_coarse))
triang = tri.Triangulation(mesh_xy[:, 0], mesh_xy[:, 1], mesh_cells)
# Plot mesh
fig, axes = plt.subplots(nrows=1, ncols=1, squeeze=False)
axes[0, 0].triplot(triang, linewidth=2.5)
axes[0, 0].axis('scaled')
axes[0, 0].axis("off")
fig.tight_layout(pad=0, w_pad=0, h_pad=0)
plt.show()
Problem in Matrix Form¶
The choice of basis functions which induced the mesh lead to a discretization of the Dirichlet problem. This discretization can be represented as a linear system \(A \mathbf{u} = \mathbf{f}\).
[3]:
import scipy.sparse
# Load linear system from file
A = scipy.sparse.load_npz(file=data_path + "matrix_poisson_res{}.npz".format(mesh_resolution_coarse))
f = np.load(file=data_path + "rhs_poisson_res{}.npy".format(mesh_resolution_coarse))
print(np.linalg.cond(A.todense()))
print("Problem dimension: {}".format(A.shape[0]))
24.789619046465486
Problem dimension: 74
The resulting matrix is sparse. The diagonal in the upper left arises from the boundary conditions.
[4]:
import matplotlib
# Plot discretized differential operator
vmax = np.max(np.abs(A.todense()))
vmin = -vmax
fig, axes = plt.subplots(nrows=1, ncols=1, squeeze=False)
axes[0, 0].imshow(A.todense(), cmap="bwr", vmin=vmin, vmax=vmax)
axes[0, 0].set_xticks([])
axes[0, 0].set_yticks([])
plt.tight_layout()
Solution using a Probabilistic Linear Solver¶
We will now solve this problem using a probabilistic linear solver. This allows us to quantify numerical uncertainty and to incorporate prior knowledge about the solution.
[5]:
from probnum.linalg import problinsolve
# Solve linear system
uhat, Ahat, Ainvhat, info = problinsolve(A, f)
print(info)
# Save solution to file
np.save(data_path + "solution_res{}.npy".format(mesh_resolution_coarse), uhat.mean)
{'iter': 21, 'maxiter': 740, 'resid_l2norm': 4.095169573053468e-05, 'trace_sol_cov': 20.197255365792024, 'conv_crit': 'resid_rtol', 'rel_cond': None}
Estimated Solution¶
We plot the estimated solution fo the solver and some samples from the estimated distribution over the solution.
[6]:
# Load mesh data
mesh_xy = np.load(file=data_path + "mesh_xy_res{}.npy".format(mesh_resolution_coarse))
mesh_cells = np.load(file=data_path + "mesh_cells_res{}.npy".format(mesh_resolution_coarse))
# Interpolation levels and color scaling
interpol_levels = None
vmin_tricont = np.min(uhat.mean.ravel())
vmax_tricont = np.max(uhat.mean.ravel())
# Sample from Posterior
np.random.seed(42)
n_samples = 2
usamples = uhat.sample(n_samples)
soldict = {"$\\bm{\\mathsf{u}}_1$" : usamples[0], "$\\bm{\\mathsf{u}}_2$" : usamples[1]}
# Figure with gridspecs
fig = plt.figure(figsize=(5.8, 4), constrained_layout=False)
gs0 = fig.add_gridspec(1, 2, width_ratios=[2,1])
gs00 = gs0[0].subgridspec(1, 1, wspace=0)
gs01 = gs0[1].subgridspec(2, 1, wspace=0)
# Estimated solution
fig.add_subplot(gs00[0, 0])
plt.tricontourf(triang, uhat.mean.ravel(), levels=interpol_levels, vmin=vmin_tricont, vmax=vmax_tricont)
plt.triplot(triang, linewidth=2.5, color="white", alpha=.25)
plt.title("$\mathbb{E}[\\bm{\\mathsf{u}}]$", fontsize=24)
plt.axis('scaled')
plt.axis("off")
# Samples
for i, (title, uvec) in enumerate(soldict.items()):
fig.add_subplot(gs01[i, 0])
plt.tricontourf(triang, uvec, levels=interpol_levels, vmin=vmin_tricont, vmax=vmax_tricont)
plt.triplot(triang, linewidth=1.25, color="white", alpha=.25)
plt.title(title, fontsize=24, y=.95)
plt.axis('scaled')
plt.axis("off")
fig.tight_layout(pad=0, w_pad=0, h_pad=0)
Estimated Linear Operators¶
We can also inspect the estimate of the discretized differential operator itself and its inverse.
[7]:
# Plot linear operators
matdict = {"$A$" : A.todense(), "$\mathbb{E}(\mathsf{A})$" : Ahat.mean.todense(),
"$\mathbb{E}(\mathsf{H})$" : Ainvhat.mean.todense()}
# Set imshow limits uniform across matrices
vmax = np.max([np.max(np.abs(mat)) for mat in list(matdict.values())])
vmin = -vmax
fig, axes = plt.subplots(nrows=1, ncols=len(matdict), figsize=(8, 4), squeeze=False, sharex=True, sharey=True)
for i, (title, mat) in enumerate(matdict.items()):
# Plot matrix
axes[0, i].imshow(mat, cmap="bwr", vmin=vmin, vmax=vmax)
axes[0, i].set_title(title, fontsize=24)
axes[0, i].set_xticks([])
axes[0, i].set_yticks([])
# Set axis color
for child in axes[0, i].get_children():
if isinstance(child, matplotlib.spines.Spine):
child.set_color("gray")
plt.tight_layout()
Error and Uncertainty Calibration¶
To see how well our computed solution approximates the “true” solution, we solved the same problem on a much finer mesh and interpolated to the coarse mesh. This allows us to compute errors and check the calibration of the uncertainty returned by the probabilistic linear solver.
[8]:
# Load (interpolated) solution from file
u = np.load(file=data_path + "solution_res{}.npy".format(mesh_resolution_coarse))
u_interpol = np.load(file=data_path + "solution_interpol_res{}tores{}.npy".format(mesh_resolution_fine,
mesh_resolution_coarse))
We compute the signed relative error \(\lVert u \rVert^{-1}(u - \mathbb{E}[\mathsf{u}])\) and the scaled error \((\Sigma^+)^{-\frac{1}{2}}(u - \mathbb{E}[\mathsf{u}])\), where \(\Sigma\) is the covariance matrix of \(\mathsf{u}\).
[9]:
# Compute error
error_u = u_interpol - u
rel_error_u = error_u / np.linalg.norm(u)
# Scaled error
Sigma = uhat.cov.todense()
Sigma_pseudoinv = np.linalg.pinv(Sigma, hermitian=True)
scaled_error_u = np.real_if_close(scipy.linalg.sqrtm(Sigma_pseudoinv) @ error_u, tol=10**8)
[10]:
# Colormap scale
vmax = np.max(np.abs(rel_error_u))
vmin = -vmax
# Plot
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(4.5, 4), squeeze=False)
triang = tri.Triangulation(mesh_xy[:, 0], mesh_xy[:, 1], mesh_cells)
tricon = plt.tricontourf(triang, rel_error_u, levels=50, cmap="bwr", vmin=vmin, vmax=vmax)
cbar = plt.colorbar(tricon, pad=0)
cbar.ax.tick_params(labelsize=18)
plt.triplot(triang, linewidth=2.5, color="white", alpha=.25)
plt.title("$\\lVert \\bm{u} \\rVert^{-1} (\\bm{u} - \mathbb{E}[\\bm{\\mathsf{u}}])$", fontsize=24)
plt.axis('scaled')
plt.axis("off")
fig.tight_layout(pad=0, w_pad=0, h_pad=0)
If the uncertainty of the probabilistic linear solver is well calibrated, then the scaled error \((\Sigma^+)^{-\frac{1}{2}}(u - \mathbb{E}[\mathsf{u}])\) is a sample of a standard normal distribution.
[11]:
# Colormap scale
vmax = np.max(np.abs(scaled_error_u))
vmin = -vmax
# Plot scaled error
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(4.5, 4), squeeze=False)
triang = tri.Triangulation(mesh_xy[:, 0], mesh_xy[:, 1], mesh_cells)
tricon = plt.tricontourf(triang, scaled_error_u, levels=50, cmap="bwr", vmin=vmin, vmax=vmax)
cbar = plt.colorbar(tricon, pad=0)
cbar.ax.tick_params(labelsize=18)
plt.triplot(triang, linewidth=2.5, color="white", alpha=.25)
plt.title("$(\\Sigma^{+})^{-\\frac{1}{2}} ( \\bm{u} - \mathbb{E}[\\bm{\\mathsf{u}}] )$", fontsize=24)
plt.axis('scaled')
plt.axis("off")
fig.tight_layout(pad=0, w_pad=0, h_pad=0)
[ ]: