This notebook illustrates the use of monotone finite difference schemes to compute viscosity solutions of non-linear PDEs, in two space dimensions. We consider Pucci's operator $$ \Lambda u(x) := \alpha(x) \lambda_{\max}(\nabla^2 u(x)) + \lambda_{\min}(\nabla^2 u(x)) $$ in the PDE $$ {-} \Lambda u(x) + \beta(x) = 0, $$ with Dirichlet boundary conditions. The PDE parameters are a positive function $\alpha$, and an arbitrary function $\beta$. We denote by $\lambda_{\max}(M)$ and $\lambda_{\min}(M)$ the largest and smallest eigenvalue of a positive definite tensor $M$. More details on this problem below.
We design two monotone numerical schemes:
The two schemes involves adaptive stencils, built using techniques from lattice geometry. The techniques developed are fairly general, and can be applied to a wide range of non-linear PDEs. Numerical implementation is kept simple thanks to the use of automatic differentiation.
Assume without loss of generality that $\alpha \leq 1$. Then for any positive definite matrix $M$ one has $$ \alpha \lambda_{\max}(M) + \lambda_{\min}(M) = \min_{0 \leq \theta \leq \pi} \mathrm{Tr}(D_\alpha(\theta) M), $$ where we denoted, with $e(\theta) := (\cos \theta, \sin \theta)$ $$ D_\alpha(\theta) = \alpha\, e(\theta) e(\theta)^T + e(\theta)^\perp (e(\theta)^\perp)^T, $$ the symmetric matrix whose eigenvalues are $\alpha$ and $1$, the former associated with the eigenvector $e(\theta)$.
Remark on the range of the variable $\theta$. For any $\theta\in \mathbb R$, one has $e(\theta+\pi) = -e(\theta)$, and therefore $D_\alpha(\theta+\pi) = D_\alpha(\theta)$. By periodicity, we may therefore limit our attention to the interval $[0,\pi]$.
Remark on the case $\alpha\geq 1$. This second case is handled by replacing the minimum over $\theta\in [0,\pi]$ with a maximum. This does not induce any additional difficulty from the theoretical or numerical standpoints. However, for the sake of simplicity, we make the assumption that $\alpha\leq 1$ in the following.
Let $K$ be a positive integer, and let $\theta_1 \leq \cdots \leq \theta_K$ be a sampling of the interval $[0,\pi]$. Then we may consider the approximate operator $$ \Lambda_K u(x) := \min_{1 \leq k \leq K} \mathrm{Tr} (D_\alpha(\theta_k) \nabla^2 u(x)). $$ Introduce decompositions of the tensors, obtained e.g. by Selling's method, $$ D_\alpha(\theta_k) = \sum_{1 \leq i \leq n} \mu_{ki} e_{ki} e_{ki}^T, $$ where $\mu_{ki} \geq 0$ and $e_{ki}$ has integer coordinates. Then we obtain the monotone numerical scheme $$ \min_{1 \leq k \leq K} \sum_{1 \leq i \leq n} \mu_{ki} \frac{ u(x+h e_{ki}) - 2 u(x) +u(x-h e_{ki})} {h^2}. $$ A consistency defect remains, which can be estimated in terms of the width of the sampling $\theta_1,\cdots,\theta_K$ of the control space $[0,\pi]$.
An additional problem is that the numerical scheme cost increases as $K$ increases. This issue becomes more acute in the case of a multi-dimensional control space.
In order to introduce this discretization, we need to recall some elements from lattice geometry. Selling's decomposition of a tensor $D$ involves a geometrical object, referred to as a $D$-obtuse superbase and here denoted $$ \mathrm{osb}(D). $$ The obtuse superbase $s=\mathrm{osb}(D)$ dictates the support $(e_{si})_{i=1}^n$ of Selling's decomposition of $D$, hence the stencil of the numerical scheme. We can take advantage of this fact to rewrite the operator as $$ \Lambda u(x) = \min_{s \in S} \Lambda_s u(x) $$ where $$ \Lambda_s u(x) := \min_{\theta, \mathrm{osb}(D_\alpha(\theta)) = s} \mathrm{Tr} (D_\alpha(\theta) \nabla^2 u). $$ Each operator $\Lambda_s$ admits the consistent discretization $$ \Lambda_s u(x) \approx \min_{\theta, \mathrm{osb}(D_\alpha(\theta)) = s} \sum_{1 \leq i \leq n} \mu_{si}(\theta) \frac{u(x+h e_{si}) - 2 u(x) + u(x-e_{si})} {h^2}, $$ and a closed form can be obtained for the r.h.s. by examining a simple optimization problem.
Summary of volume Non-Divergence form PDEs, this series of notebooks.
Main summary of the Adaptive Grid Discretizations book of notebooks, including the other volumes.
Acknowledgement. Some of the experiments presented in these notebooks are part of ongoing research with Ludovic Métivier and Da Chen.
Copyright Jean-Marie Mirebeau, Centre Borelli, ENS Paris-Saclay, CNRS, University Paris-Saclay
import sys; sys.path.insert(0,"..") # Allow imports from parent directory
#from Miscellaneous import TocTools; print(TocTools.displayTOC('NonlinearMonotoneSecond2D','NonDiv'))
from agd import Selling
from agd import LinearParallel as lp
from agd import AutomaticDifferentiation as ad
from agd import Domain
from agd.Plotting import savefig; #savefig.dirName = "Figures/NonlinearMonotoneSecond2D"
import numpy as np
import matplotlib.pyplot as plt
Some utility functions
newton_root = ad.Optimization.newton_root
norm = ad.Optimization.norm
def BoundaryNeighborhood(interior,width=1):
bd=interior.copy()
bd[0,:]=False; bd[-1,:]=False; bd[:,0]=False; bd[:,-1]=False
directions = ( (0,0),(0,1),(0,-1),(1,0),(-1,0) )
neigh = np.stack(tuple(np.roll(bd,e,axis=(0,1)) for e in directions),axis=0)
neigh = np.logical_and(neigh.any(axis=0),np.logical_not(neigh).any(axis=0))
for i in range(width):
neigh = np.stack(tuple(np.roll(neigh,e,axis=(0,1)) for e in directions),axis=0).any(axis=0)
return neigh
We present a naive, non-monotone discretization of the addressed problem. This scheme can be used to check that a numerical solution (produced by other means) is correct, or to construct synthetic examples. However, using it to solve the PDE is usually bound to fail.
The naive scheme is based on a reconstruction of the Hessian matrix of the form $$ \begin{pmatrix} D^h_{00} u(x) & D^h_{01} u(x)\\ D^h_{01} u(x) & D^h_{11} u(x) \end{pmatrix}, $$ where $D_{00}$, $D_{01}$ and $D_{11}$ are finite-difference operators. Namely $$ D^h_{00} u(x) := \frac{u(x_0+h,x_1)-2 u(x_0,x_1) + u(x_0-h,x_1)}{h^2}, $$ likewise for $D^h_{11} u(x)$, and finally $$ D^h_{01} u(x) := \frac{u(x_0+h,x_1+h)-u(x_0-h,x_1+h)-u(x_0+h,x_1-h)+h(x_0-h,x_1-h)}{4 h^2}. $$
def SchemeNonMonotone(u,α,β,bc,sqrt_relax=1e-6):
# Compute the hessian matrix of u
uxx = bc.Diff2(u,(1,0))
uyy = bc.Diff2(u,(0,1))
uxy = 0.25*(bc.Diff2(u,(1,1)) - bc.Diff2(u,(1,-1)))
# Compute the eigenvalues
# The relaxation is here to tame the non-differentiability of the square root.
htr = (uxx+uyy)/2. # Half trace
Δ = ((uxx-uyy)/2.)**2 + uxy**2 # Discriminant of characteristic polynomial
sΔ = np.sqrt( np.maximum( Δ, sqrt_relax) )
λ_max = htr+sΔ
λ_min = htr-sΔ
# Numerical scheme
residue = β - α*λ_max - λ_min
# Boundary conditions
return np.where(bc.interior,residue,u-bc.grid_values)
Our next step is to define the parameters of our specific problem. Regarding the boundary conditions, we set $u=0$ on the square boundary, and $u=-1$ on some interior diamond. For well posedness, $d$ and $\alpha$ must be positive over the domain.
# Create the domain
aX0,dx = np.linspace(-1,1,100,retstep=True); aX1=aX0;
X = np.array(np.meshgrid(aX0,aX1,indexing='ij'))
# Set the boundary conditions
bc_grid_values=np.full(X.shape[1:],np.nan)
bc_grid_values[ad.Optimization.norm(X,ord=1,axis=0) < 0.4] = -1
bc = Domain.MockDirichlet(bc_grid_values,dx,padding=0.)
# Choose the PDE parameters
α = 0.25
β = 1
The naive and non-monotone discretization scheme is consistent, but lacks any other sort of theoretical guarantees. It is pure luck that the Newton method does converge in this simple instance, and that the result looks reasonable.
params = (α,β,bc); guess = np.zeros(bc.shape);
solution = ad.Optimization.newton_root(SchemeNonMonotone,guess,params)
Iteration: 1 Residue norm: 26.028138597195365 Iteration: 2 Residue norm: 0.9663557973411793 Iteration: 3 Residue norm: 0.058931864127965694 Iteration: 4 Residue norm: 0.0004164739393537342 Iteration: 5 Residue norm: 7.70993752174931e-08 Iteration: 6 Residue norm: 1.801669924361704e-12 Target residue reached. Terminating.
plt.title("Solution obtained with non-monotone scheme")
plt.contourf(*X,solution); plt.axis('equal');
We present a numerical scheme based on sampling the control space, which is quite simple and generic. Given fields of positive definite diffusion tensors $D_k(x)$, $1 \leq k \leq K$, without any specific assumption on their origin, we compute the decompositions $$ D_k(x) = \sum_{1 \leq i \leq n} \mu_{ki}(x) e_{ki}(x) e_{ki}(x)^T. $$ We then implement the scheme $$ \beta(x) - \min_{1\leq k \leq K} \sum_{1 \leq i \leq n} \mu_{ki}(x) \frac{u(x+h e_{ki}(x))-2 u(x) +u(x-h e_{ki}(x))}{h^2}. $$
def SchemeSampling(u,diffs,β,bc):
# Tensor decomposition
μ,e = Selling.Decomposition(diffs)
# Numerical scheme
μ = bc.as_field(μ)
residue = β - (μ*bc.Diff2(u,e)).sum(0).min(0)
# Boundary conditions
return np.where(bc.interior,residue,u-bc.grid_values)
The tensors involved in our PDE take the following form.
def Diff(α,θ):
e0 = np.array(( np.cos(θ),np.sin(θ)))
e1 = np.array((-np.sin(θ),np.cos(θ)))
if isinstance(α,np.ndarray):
e0,e1 = (as_field(e,α.shape) for e in (e0,e1))
return α*lp.outer_self(e0) + lp.outer_self(e1)
We also choose a discretization of the control space.
nθ = 20
θs = np.linspace(0,np.pi,nθ,endpoint=False)
We next solve the PDE and display the solution.
params = (Diff(α,θs), β,bc)
solution = newton_root(SchemeSampling,guess,params)
Iteration: 1 Residue norm: 69.11435896930232 Iteration: 2 Residue norm: 1.6659074845554858 Iteration: 3 Residue norm: 0.25230862927071485 Iteration: 4 Residue norm: 0.05531949820048643 Iteration: 5 Residue norm: 0.0034489933490459146 Iteration: 6 Residue norm: 1.4779288903810084e-12 Target residue reached. Terminating.
plt.axis('equal'); plt.title('Solution to lambda_max/2 + lambda_min = 1')
plt.contourf(*X,solution);
Since the solution is not explicit, we use the non-monotone numerical scheme to test the result. We eliminate a layer around the neighborhood of the boundary conditions, where the solution is not smooth.
residue_non_monotone = SchemeNonMonotone(solution,α,β,bc)
residue_non_monotone[BoundaryNeighborhood(bc.interior,width=5)] = 0.
print("Max cross-residue of the sampling based numerical solution:",norm(residue_non_monotone,ord=np.inf))
plt.axis('equal'); plt.title("Cross-residue of the sampling based numerical solution:")
plt.contourf(*X,residue_non_monotone);
Max cross-residue of the sampling based numerical solution: 0.05897075265381302
If one chooses $\alpha=1$, then the PDE becomes linear, namely $-\Delta u + \beta = 0$. In this very specific case, the sampling based and non-monotone scheme coincide, with the usual discretization of the laplacian operator. As a result the Newton method converges in one iteration, and the cross-residue vanishes.
params = (Diff(1.,θs), β,bc)
solution = newton_root(SchemeSampling,guess,params)
plt.axis('equal'); plt.title('Solution to λ_max + λ_min = 1')
plt.contourf(*X,solution);
Iteration: 1 Residue norm: 6.366462912410498e-12 Target residue reached. Terminating.
print("Cross residue in special case α=1 :", norm(SchemeNonMonotone(solution,1.,β,bc),ord=np.inf) )
Cross residue in special case α=1 : 6.366462912410498e-12
In contrast, if one chooses a very small value of $\alpha$, then the PDE becomes more and more non-linear, which raises numerical difficulties discussed below. If in addition $\beta = 0$, then we recover PDE characterization of the convex envelope: $$ -\lambda_{\min}(\nabla^2 u) = 0 $$
Numerical challenges. As $\alpha\to 0$, the condition number of the tensors $D_\alpha(\theta)$ increase. A finer sampling of the interval $[0,\pi]$ is required, which increases the numerical cost of the method. In addition, the width of the discretization stencil increases, and therefore the effective discretization scale is reduced.
Note on computing the convex envelope. The computation of convex envelopes is one of the most central problems in algorithmic geometry. For instance, Voronoi diagrams are deduced from a convex envelope computation in higher dimension. Extremely efficient software packages are available for this problem, and PDE methods are not the recommended way to go.
α_small = 0.01
nθ_small = 50
θs_small = np.linspace(0,np.pi,nθ_small,endpoint=False)
β_cvx_env = 0.
The above parameters turn the solution into the convex envelope of the boundary conditions. Recall that we imposed:
params = (Diff(α_small,θs_small), β_cvx_env,bc)
solution = newton_root(SchemeSampling,guess,params)
plt.axis('equal'); plt.title('Convex envelope of the boundary conditions.')
plt.contourf(*X,solution);
Iteration: 1 Residue norm: 551.0630896254976 Iteration: 2 Residue norm: 9.991945726638988 Iteration: 3 Residue norm: 12.63487784876113 Iteration: 4 Residue norm: 3.510662245585972 Iteration: 5 Residue norm: 0.732778858518305 Iteration: 6 Residue norm: 0.10918596570734207 Iteration: 8 Residue norm: 0.0029157371512060527 Iteration: 10 Residue norm: 1.8987381866395828e-05 Iteration: 12 Residue norm: 3.823126363027734e-08 Iteration: 14 Residue norm: 9.966120374337267e-09 Target residue reached. Terminating.
The solution is piecewise affine, and its gradient is piecewise constant. Note in particular that the solution is not twice differentiable, and the equation $-\lambda_{\min}(\nabla^2 u) = 0$ here only has meaning in the sense of viscosity solutions.
grad = np.array(np.gradient(solution,bc.gridscale))
plt.axis('equal'); plt.title("The solution to the convex envelope problem is piecewise affine")
plt.contourf(*X,norm(grad,ord=2,axis=0));
We next negate the boundary conditions, imposing:
This raises an apparent incompatibility: the convex envelope should be $u=0$ on the whole square, yet we impose a distinct value in the diamond. What to expect ?
bc_negated = Domain.MockDirichlet(-bc_grid_values,dx,padding=0.)
params = (Diff(α_small,θs_small), β_cvx_env,bc_negated)
solution = newton_root(SchemeSampling,guess,params)
plt.axis('equal'); plt.title('Convex envelope of the boundary conditions.')
plt.contourf(*X,solution);
Iteration: 1 Residue norm: 806.0820465420255 Iteration: 2 Residue norm: 487.7132470778603 Iteration: 3 Residue norm: 200.2083809007435 Iteration: 4 Residue norm: 164.31875614668522 Iteration: 5 Residue norm: 111.67932913072441 Iteration: 6 Residue norm: 75.49541847970949 Iteration: 8 Residue norm: 29.94570387126228 Iteration: 10 Residue norm: 1.461053500406706e-13 Target residue reached. Terminating.
The solution to the second order PDE $-\lambda_{\min}(\nabla^2 u)$ is still unique and well defined, in the sense of viscosity solutions, with these boundary conditions. It is piecewise constant, with value $u=0$ except on the interior diamond where we impose $u=-1$.
The numerical scheme implemented in this section involves a maximization over a set of diffusion tensors. The size of this set dictates the accuracy of the method, and it may therefore become rather large. In combination with the overhead of sparse AD, this can increase the numerical cost.
We can significantly limit the numerical cost using a technique based on the envelope theorem, which proceeds in two steps:
Important : other optimization opportunities. The "optimization" presented in this subsection only serves to illustrate the envelope theorem mechanism. It is not effective in terms of computation time, because the optimized part is not dominant. There are other optimization opportunities here, the most obvious one being to avoid recomputing the tensor decompositions at each call of the iterative solver. The choice of linear solver may also be of importance.
def SchemeSampling_OptInner(u,diffs,bc,oracle=None):
# Select the active tensors, if they are known
if not(oracle is None):
diffs = np.take_along_axis(diffs, np.broadcast_to(oracle,diffs.shape[:2]+(1,)+oracle.shape),axis=2)
print("Has AD information :", ad.is_ad(u), ". Number active tensors per point :", diffs.shape[2])
# Tensor decomposition
coefs,offsets = Selling.Decomposition(diffs)
# Return the minimal value, and the minimizing index
return ad.min_argmin( lp.dot_VV(coefs,bc.Diff2(u,offsets)), axis=0)
def SchemeSampling_Opt(u,diffs,β,bc):
# Evaluate the operator using the envelope theorem
result,_ = ad.apply(SchemeSampling_OptInner, u,bc.as_field(diffs),bc, envelope=True)
# Boundary conditions
return np.where(bc.interior, β-result, u-bc.grid_values)
params = (Diff(α,θs), β,bc)
solution = newton_root(SchemeSampling_Opt,guess,params)
Has AD information : False . Number active tensors per point : 20 Has AD information : True . Number active tensors per point : 1 Has AD information : False . Number active tensors per point : 20 Has AD information : True . Number active tensors per point : 1 Iteration: 1 Residue norm: 69.11435896930232 Has AD information : False . Number active tensors per point : 20 Has AD information : True . Number active tensors per point : 1 Iteration: 2 Residue norm: 1.6659074845554858 Has AD information : False . Number active tensors per point : 20 Has AD information : True . Number active tensors per point : 1 Iteration: 3 Residue norm: 0.25230862927071485 Has AD information : False . Number active tensors per point : 20 Has AD information : True . Number active tensors per point : 1 Iteration: 4 Residue norm: 0.05531949820048643 Has AD information : False . Number active tensors per point : 20 Has AD information : True . Number active tensors per point : 1 Iteration: 5 Residue norm: 0.0034489933490459146 Has AD information : False . Number active tensors per point : 20 Has AD information : True . Number active tensors per point : 1 Iteration: 6 Residue norm: 1.4779288903810084e-12 Target residue reached. Terminating.
Setting up a monotone and consistent discretization requires a bit more work, but is worthwhile in the end if performance and accuracy are a target. Let us recall that the diffusion tensors take the form, $$ D_\alpha(\theta) = \alpha\, e(\theta) e(\theta)^T + e(\theta)^\perp (e(\theta)^\perp)^T, $$ where $0< \alpha \leq 1$ is a fixed parameter in the following, and $\theta \in [0,\pi]$.
The first, and main, difficulty is to construct a sequence of angles $0 = \theta_0 \leq \cdots \leq \theta_N = \pi$ and of superbases $s_0,\cdots, s_{N-1}$ such that $$ s_k \text{ is } D_\theta(\theta) \text{-obtuse, for all } \theta \in [\theta_k, \theta_{k+1}]. $$ For that purpose, we remark that $$ D_\alpha(\theta) = D_0 + D_1 \cos(2 \theta) + D_2 \sin(2 \theta), $$ where (omitting the dependency on $\alpha$ for readability) $$ D_0 = \frac{\alpha+1} 2 \begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix}, \quad D_1 = \frac{\alpha-1} 2 \begin{pmatrix} 1 & 0\\ 0 &-1 \end{pmatrix}, \quad D_2 = \frac{\alpha-1} 2 \begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix}. $$ Then, for given $u,v \in R^2$, one has $$ <u,D_\alpha(\theta) v> = a_0 + a_1 \cos(2 \theta) + a_2 \sin(2 \theta) = r (\cos(2\theta-\phi) - c). $$ where $a_i = <u,D_i v>$. Then $r e(\phi) = (a_1,a_2)$, and $c=-a_0/r$. We assume that $r$ is positive. Eventually, the above scalar product is
def MakeD(α):
return np.moveaxis(0.5*np.array([
(α+1)*np.array([[1,0],[0,1]]),
(α-1)*np.array([[1,0],[0,-1]]),
(α-1)*np.array([[0,1],[1,0]])
]), 0,-1)
def NextAngleAndSuperbase(θ,sb,D):
pairs = np.stack([(1,2), (2,0), (0,1)],axis=1)
scals = lp.dot_VAV(np.expand_dims(sb[:,pairs[0]],axis=1),
np.expand_dims(D,axis=-1), np.expand_dims(sb[:,pairs[1]],axis=1))
ϕ = np.arctan2(scals[2],scals[1])
cst = -scals[0]/np.sqrt(scals[1]**2+scals[2]**2)
θ_max = np.pi*np.ones(3)
mask = cst<1
θ_max[mask] = (ϕ[mask]-np.arccos(cst[mask]))/2
θ_max[θ_max<=0] += np.pi
θ_max[θ_max<=θ] = np.pi
k = np.argmin(θ_max)
i,j = (k+1)%3,(k+2)%3
return (θ_max[k],np.stack([sb[:,i],-sb[:,j],sb[:,j]-sb[:,i]],axis=1))
def AnglesAndSuperbases(D,maxiter=200):
sb = Selling.CanonicalSuperbase(np.eye(2)).astype(int)
θs=[]
superbases=[]
θ=0
for i in range(maxiter):
θs.append(θ)
if(θ>=np.pi): break
superbases.append(sb)
θ,sb = NextAngleAndSuperbase(θ,sb,D)
return np.array(θs), np.stack(superbases,axis=2)
The above code is a bit intricate, but its purpose is simple : split the interval $[0,\pi]$ into sub-intervals on which the support of Selling's decomposition of the tensors is fixed and known.
α=0.1
θs,superbases = AnglesAndSuperbases(MakeD(α))
θs
array([0. , 0.12912139, 0.65627678, 0.91451955, 1.44167494, 1.57079633, 1.69991771, 2.2270731 , 2.48531588, 3.01247127, 3.14159265])
superbases
array([[[-1, -1, -1, 0, 0, 1, 0, 0, -1, -1], [ 1, -1, 1, 1, -1, 0, 1, -1, -1, 1], [ 0, 2, 0, -1, 1, -1, -1, 1, 2, 0]], [[-1, -1, -1, 1, 1, 0, -1, -1, 1, 1], [ 0, 0, 0, 1, -1, -1, -1, 1, 0, 0], [ 1, 1, 1, -2, 0, 1, 2, 0, -1, -1]]])
θs_sampled = np.linspace(0,np.pi,200)
decomp = Selling.GatherByOffset(θs_sampled,*Selling.Decomposition(Diff(α,θs_sampled)))
fig = plt.figure(figsize=(20,10)) # Paper:(10,5)
plt.title(f"Decomposition of a rotated matrix of eigenvalues {α} and {1}")
plt.xlabel("θ"); plt.ylabel("Coefficient")
for offset,(angle,coef) in decomp.items():
plt.plot(angle,coef)
plt.legend(decomp.keys());
for θ in θs: # Show a vertical line for each angle theta where the stencil changes
plt.axvline(x=θ)
savefig(fig,"DecompositionCoefficients.png")
A value of $\alpha$ closer to $1$ yields a smaller number of superbases.
α=1/4
θs_sampled = np.linspace(0,np.pi,200)
θs,superbases = AnglesAndSuperbases(MakeD(α))
decomp = Selling.GatherByOffset(θs_sampled,*Selling.Decomposition(Diff(α,θs_sampled)))
fig = plt.figure(figsize=(20,10)) # Paper:(10,5)
plt.title(f"Decomposition of a rotated matrix of eigenvalues {α} and {1}")
plt.xlabel("θ"); plt.ylabel("Coefficient")
for offset,(angle,coef) in decomp.items():
plt.plot(angle,coef)
plt.legend(decomp.keys());
for θ in θs: # Show a vertical line for each angle theta where the stencil changes
plt.axvline(x=θ);
Whereas a smaller values yields more superbases and increases the numerical scheme complexity.
α=1/100
θs_sampled = np.linspace(0,np.pi,400)
θs,superbases = AnglesAndSuperbases(MakeD(α))
decomp = Selling.GatherByOffset(θs_sampled,*Selling.Decomposition(Diff(α,θs_sampled)))
fig = plt.figure(figsize=(20,10)) # Paper:(10,5)
plt.title(f"Decomposition of a rotated matrix of eigenvalues {α} and {1}")
plt.xlabel("θ"); plt.ylabel("Coefficient")
for offset,(angle,coef) in decomp.items():
plt.plot(angle,coef)
plt.legend(decomp.keys());
for θ in θs: # Show a vertical line for each angle theta where the stencil changes
plt.axvline(x=θ);
The second step is discretize $$ \min_{\theta \in [\theta_0,\theta_1]} \mathrm{Tr}(D_\alpha(\theta) \nabla^2 u(x)), $$ when $D_\alpha(\theta)$ admits the same known obtuse superbase $b = (b_0,b_1,b_2)$ for each $\theta \in [\theta_0,\theta_1]$.
The discretization reads $$ \min_{\theta \in [\theta_0, \theta_1]} - \sum_{1 \leq i \leq 3} <b_{i+1}, D_\alpha(\theta) b_{i+2}> \frac{u(x+h b_i^\perp)-2 u(x) + u(x-h b_i^\perp)}{h^2}. $$ where indices of the superbase are understood modulo $3$.
For that purpose, we rely on the expression $D_\alpha(\theta) = D_0 + D_1 \cos(2 \theta) + D_2 \sin(2 \theta)$, and on the explicit solution $$ \min_{\phi \in [2\theta_0,2\theta_1]} d_0 + d_1 \cos \phi + d_2 \sin \phi = d_0 - \sqrt{d_1^2+d_2^2} $$ if $(\cos \phi,\sin \phi)$ is proportional to $-(d_1,d_2)$ for some $\phi \in [2\theta_0,2\theta_1]$. Otherwise the minimum is attained at $2\theta_0$ or $2 \theta_1$.
Theoretical issue for the Newton method The lack of differentiability of the term $\sqrt{d_1^2+d_2^2}$ is a theoretical issue for the Newton method. It does not raise any difficulty from a practical standpoint, although numpy does raise a warning on the matter.
def MinimizeTrace(u,α,bc,sqrt_relax=1e-16):
# Compute the tensor decompositions
D=MakeD(α)
θ,sb = AnglesAndSuperbases(D)
θ = np.array([θ[:-1],θ[1:]])
# Compute the second order differences in the direction orthogonal to the superbase
sb_rotated = np.array([-sb[1],sb[0]])
d2u = bc.Diff2(u,sb_rotated)
d2u[...,bc.not_interior]=0. # Placeholder values to silent NaNs
# Compute the coefficients of the tensor decompositions
sb1,sb2 = np.roll(sb,1,axis=1), np.roll(sb,2,axis=1)
sb1,sb2 = (e.reshape( (2,3,1)+sb.shape[2:]) for e in (sb1,sb2))
D = D.reshape((2,2,1,3,1)+D.shape[3:])
# Axes of D are space,space,index of superbase element, index of D, index of superbase, and possibly shape of u
scals = lp.dot_VAV(sb1,D,sb2)
# Compute the coefficients of the trigonometric polynomial
scals,θ = (bc.as_field(e) for e in (scals,θ))
coefs = -lp.dot_VV(scals, np.expand_dims(d2u,axis=1))
# Optimality condition for the trigonometric polynomial in the interior
value = coefs[0] - np.sqrt(np.maximum(coefs[1]**2+coefs[2]**2,sqrt_relax))
coefs_ = ad.remove_ad(coefs) # removed AD information
angle = np.arctan2(-coefs_[2],-coefs_[1])/2.
angle[angle<0]+=np.pi
# Boundary conditions for the trigonometric polynomial minimization
mask = np.logical_not(np.logical_and(θ[0]<=angle,angle<=θ[1]))
t,c = θ[:,mask],coefs[:,mask]
value[mask],amin_t = ad.min_argmin(c[0]+c[1]*np.cos(2*t)+c[2]*np.sin(2*t),axis=0)
# Minimize over superbases
value,amin_sb = ad.min_argmin(value,axis=0)
# Record the optimal angles for future use
angle[mask]=np.take_along_axis(t,np.expand_dims(amin_t,axis=0),axis=0).squeeze(axis=0) # Min over bc
angle = np.take_along_axis(angle,np.expand_dims(amin_sb,axis=0),axis=0) # Min over superbases
return value,angle
def SchemeConsistent(u,α,β,bc):
value,_ = MinimizeTrace(u,α,bc)
residue = β - value
return np.where(bc.interior,residue,u-bc.grid_values)
The scheme is efficiently solved by the Newton method. The warning (possibly) raised is related with the lack of differentiability of sqrt, as mentioned above.
%%time
params = (α,β,bc)
guess2 = 0.5*(X[0]**2 +2.*X[1]**2)
solution = newton_root(SchemeConsistent,guess2,params)
Iteration: 1 Residue norm: 261.3488583788808 Iteration: 2 Residue norm: 18.317895546589604 Iteration: 3 Residue norm: 12.641943601329235 Iteration: 4 Residue norm: 7.204295363133955 Iteration: 5 Residue norm: 7.265577373440989 Iteration: 6 Residue norm: 4.462833234235047 Iteration: 8 Residue norm: 2.58828303024454 Iteration: 10 Residue norm: 0.9007355797389778 Iteration: 12 Residue norm: 0.0009478570336929337 Iteration: 14 Residue norm: 5.0045809718568535e-09 Target residue reached. Terminating. CPU times: user 5.35 s, sys: 1.97 s, total: 7.32 s Wall time: 7.33 s
fig = plt.figure(figsize=(4,4)); plt.axis('equal')
plt.title("Solution from the monotone and consistent scheme.")
plt.contourf(*X,solution);
savefig(fig,"SolutionMonotoneConsistent.png")
For validation, we compute the residue of the naive scheme. It is small, as expected, except on the boundary of the obstacle.
fig = plt.figure(figsize=(6,4)); plt.axis('equal')
plt.title("Residue from the monotone and consistent scheme.")
res = SchemeNonMonotone(solution,*params)
res[BoundaryNeighborhood(bc.interior,width=5)] = 0.
plt.contourf(res); plt.colorbar();
savefig(fig,"ResidueMonotoneConsistent.png")
Similarly to the sampling based scheme, we propose an enhanced implementation taking advantage of the envelope theorem, which we recall lets one differentiate functions of the form $$ F(x) = \min_{\alpha \in A} F_\alpha(x). $$ However, there is one important distinction, relative to the nature of the optimization parameter $\alpha \in A$. Indeed, it is:
def MinimizeTrace_Opt(u,α,bc,oracle=None):
if oracle is None: return MinimizeTrace(u,α,bc)
# The oracle contains the optimal angles
diffs=Diff(α,oracle.squeeze(axis=0))
coefs,sb = Selling.Decomposition(diffs)
value = lp.dot_VV(coefs,bc.Diff2(u,sb))
return value,oracle
def SchemeConsistent_Opt(u,α,β,bc):
value,_ = ad.apply(MinimizeTrace_Opt,u,α,bc,envelope=True)
residue = β - value
return np.where(bc.interior,residue,u-bc.grid_values)
%%time
params = (α,β,bc)
guess2 = 0.5*(X[0]**2 +2.*X[1]**2)
solution = newton_root(SchemeConsistent_Opt,guess2,params)
Iteration: 1 Residue norm: 261.3488583788039 Iteration: 2 Residue norm: 18.317895546579795 Iteration: 3 Residue norm: 12.641943601156168 Iteration: 4 Residue norm: 7.204295363132027 Iteration: 5 Residue norm: 7.265577373441897 Iteration: 6 Residue norm: 4.462833234235095 Iteration: 8 Residue norm: 2.588283030244323 Iteration: 10 Residue norm: 0.9007355797388015 Iteration: 12 Residue norm: 0.0009478570337320136 Iteration: 14 Residue norm: 6.481196690444335e-10 Target residue reached. Terminating. CPU times: user 1.32 s, sys: 279 ms, total: 1.6 s Wall time: 1.6 s
We illustrate (first order) accurate (Dirichlet) boundary conditions on a general domain. This is in contrast with the numerical experiments presented in the above subsections, which rely on a rather crude implementation of the boundary conditions.
The chosen domain is ring shaped, with a non-smooth boundary including reentrant corners. The chosen boundary condition is $0$ on the inner boundary, and $1$ on the outer boundary.
outer = Domain.Union(Domain.Ball(),Domain.Box([[0,1],[-1,1]]) )
inner = Domain.AffineTransform(outer,0.4*lp.rotation(np.pi/3))
domain_ring = Domain.Complement(outer,inner)
def bc_value_ring(x):
"""0 on inner boundary, 1 on outer boundary."""
return outer.level(x)+inner.level(x) > 0
bc_ring = Domain.Dirichlet(domain_ring,bc_value_ring,X)
plt.contourf(*X,domain_ring.contains(X)); plt.axis('equal');
%%time
params = (0.5**2,0.,bc_ring)
guess2 = 0.5*(X[0]**2 +2.*X[1]**2)
solution_05 = newton_root(SchemeConsistent_Opt,guess2,params)
fig = plt.figure(figsize=[4,4]); plt.axis('equal')
plt.title(r"$\lambda_{min}(\nabla^2u)+\alpha\lambda_{\max}(\nabla^2u)=0$, $\alpha=1/4$.")
plt.contourf(*X,np.where(bc_ring.interior,solution_05,np.nan));
savefig(fig,"Consistent_Ring_05.png")
Iteration: 1 Residue norm: 87.71412160918004 Iteration: 2 Residue norm: 3.5886119180257348 Iteration: 3 Residue norm: 0.20858071879401102 Iteration: 4 Residue norm: 0.009504676465152953 Iteration: 5 Residue norm: 0.00041185762405138073 Iteration: 6 Residue norm: 7.523386892629644e-07 Iteration: 7 Residue norm: 3.1645797093915462e-12 Target residue reached. Terminating. CPU times: user 253 ms, sys: 26.8 ms, total: 280 ms Wall time: 280 ms
grad = np.array(np.gradient(solution_05,bc_ring.gridscale))
grad[:,np.logical_not(bc_ring.domain.contains_ball(X,1.5*bc.gridscale))]=0.
fig = plt.figure(figsize=[4,4]); plt.axis('equal')
plt.title(r"Norm of the solution gradient, $\alpha=1/4$")
plt.pcolormesh(*X,norm(grad,ord=2,axis=0)); #plt.colorbar();
savefig(fig,"Consistent_Ring_05_GradNorm.png")
fig = plt.figure(figsize=[4,4]); plt.axis('equal')
plt.title(r"Norm of the solution gradient, $\alpha=1/4$")
plt.contourf(*X,norm(grad,ord=2,axis=0)); #plt.colorbar();
savefig(fig,"Consistent_Ring_05_GradNorm_Levels.png")
%%time
params = (0.05**2,0.,bc_ring)
guess2 = 0.5*(X[0]**2 +2.*X[1]**2)
solution_005 = newton_root(SchemeConsistent_Opt,guess2,params)
fig = plt.figure(figsize=[4,4]); plt.axis('equal')
plt.title(r"$\lambda_{min}(\nabla^2u)+\alpha\lambda_{\max}(\nabla^2u)=0$, $\alpha=1/400$.")
#Paper : plt.title(r"$\lambda_{min}(\nabla^2u)+\mu\lambda_{\max}(\nabla^2u)=0$, $\mu=1/400$.")
plt.contourf(*X,np.where(bc_ring.interior,solution_005,np.nan))
savefig(fig,"Consistent_Ring_005.png")
Iteration: 1 Residue norm: 194.99025511872824 Iteration: 2 Residue norm: 38.274702920760774 Iteration: 3 Residue norm: 15.351992164905365 Iteration: 4 Residue norm: 6.635734520481244 Iteration: 5 Residue norm: 2.2868865303011217 Iteration: 6 Residue norm: 0.7152594416761943 Iteration: 8 Residue norm: 0.026452529915581367 Iteration: 10 Residue norm: 0.0007699501225768814 Iteration: 12 Residue norm: 3.2508683114450736e-05 Iteration: 14 Residue norm: 7.098774820069139e-07 Iteration: 15 Residue norm: 6.2055968763630104e-09 Target residue reached. Terminating. CPU times: user 13.3 s, sys: 2.68 s, total: 16 s Wall time: 16 s
grad = np.array(np.gradient(solution_005,bc.gridscale))
grad[:,np.logical_not(bc_ring.domain.contains_ball(X,1.5*bc.gridscale))]=np.nan
fig = plt.figure(figsize=[5,4]); plt.axis('equal')
plt.title(r"Norm of the solution gradient, $\alpha=1/400$")
#Paper : plt.title(r"Norm of the solution gradient, $\mu=1/400$")
plt.pcolormesh(*X,norm(grad,ord=2,axis=0)); plt.colorbar();
savefig(fig,"Consistent_Ring_005_GradNorm.png")
Journal version of figures
The experiments presented below aim at informally validating the numerical implementation, by
The consistent scheme, here denoted $F$, is arguably significantly more complex to implement than numerical scheme based on a sampling of the control space, here denoted $F_n$ where $n$ is the number of samples (angles).
For cross validation, one can observe that for any discrete map $u$, one has $$ F_n(u) = F(u) + O(1/n). $$ Note, crucially, that the test function $u$ is fixed, and so is the grid scale. The continuous limit only takes place in the control space.
np.random.seed(42)
u_random = np.random.uniform(-1,1,guess.shape)
bc_unit = Domain.MockDirichlet(guess.shape,1,padding=0.)
params=(α,β,bc_unit)
residue_consistent = SchemeConsistent(u_random,*params)
def error(nθ):
θs = np.linspace(0,np.pi,nθ,endpoint=False)
params=(Diff(α,θs),β,bc_unit)
residue_sampling = SchemeSampling(u_random,*params)
return norm(residue_consistent-residue_sampling,ord=np.inf)
nθ_validation = np.array([2**n for n in range(1,10)])
error_validation = np.array([error(n) for n in nθ_validation])
Convergence of the sampling based scheme toward the continuous one is observed, with the expected convergence rate.
plt.title("Convergence of the sampling scheme toward the consistent scheme")
n_inv = 1./nθ_validation
plt.loglog(n_inv,error_validation,
n_inv,n_inv);
We rely on automatic differentiation to compute the derivatives of an analytic function, and evaluate the PDE operator of interest. We then compute the numerical scheme residue on a synthetic problem with a known solution.
def Pucci_ad(u,α,x):
"""
Computes alpha*lambda_max(D^2 u) + lambda_min(D^2 u),
at the given set of points, by automatic differentiation.
"""
x_ad = ad.Dense2.identity(constant=x,shape_free=(2,))
hessian = u(x_ad).hessian()
Δ = ((hessian[0,0]-hessian[1,1])/2.)**2 + hessian[0,1]**2
sΔ = np.sqrt(Δ)
mean = (hessian[0,0]+hessian[1,1])/2.
λMin,λMax = mean-sΔ,mean+sΔ
return λMin+α*λMax
def Residue_ad(u,α,dom,X):
bc = Domain.Dirichlet(dom,u,X)
rhs = Pucci_ad(u,α,X)
residue = SchemeConsistent(u(X),α,rhs,bc)
residue[bc.not_interior]=0
return residue
def test_quadratic(x): return x[0]**2+x[1]**2
def test_polynomial(x): return 0.5*(x[0]**2+x[1]**2)**2
dom_convex = Domain.Union(Domain.Ball(),Domain.Box())
dom_ball = Domain.Ball()
dom_square = Domain.Box()
Because the scheme is second order consistent its residue is zero up to essentially machine precision on the quadratic function. (Note that the relaxation introduced for differentiability of a square root is mainly responsible for the small error).
[norm(Residue_ad(test_quadratic,0.5,dom,X),ord=np.inf) for dom in (dom_convex,dom_ball,dom_square)]
[1.0002291883637326e-08, 1.0003609052233742e-08, 1.0002291883637326e-08]
On the non-quadratic function, we get essentially second order convergence in the $L^1$ averaged norm. A slower convergence rate is achieved in the $L^\infty$ norm, because the finite differences are only first order accurate at the boundary.
aX_25 = np.linspace(-1,1,25)
X_25 = np.array(np.meshgrid(aX_25,aX_25,indexing='ij'))
aX_50 = np.linspace(-1,1,50)
X_50 = np.array(np.meshgrid(aX_50,aX_50,indexing='ij'))
X_100=X
print("Domain (line): convex, ball, square.")
print("Resolution (column) : 25,50,100.\n")
print("Mean residue of exact polynomial solution.")
print(np.array([[
norm(Residue_ad(test_polynomial,0.5,dom,X),ord=1,averaged=True)
for X in (X_25,X_50,X_100)]
for dom in (dom_convex,dom_ball,dom_square)]))
print("Residue of exact polynomial solution in the L^Infinity norm.")
print(np.array([[
norm(Residue_ad(test_polynomial,0.5,dom,X),ord=np.inf)
for X in (X_25,X_50,X_100)]
for dom in (dom_convex,dom_ball,dom_square)]))
Domain (line): convex, ball, square. Resolution (column) : 25,50,100. Mean residue of exact polynomial solution. [[0.01149697 0.00395542 0.00096948] [0.01172253 0.00432875 0.00105138] [0.00251758 0.00068901 0.00017824]] Residue of exact polynomial solution in the L^Infinity norm. [[0.12456117 0.11640453 0.06062434] [0.12456117 0.11640453 0.06062434] [0.01388889 0.00333195 0.00081624]]