Prescribe solution of a PDE at specific points

Question

I am using MATLAB's PDE toolbox to solve the differential equation

$-\nabla\cdot\left(c(x)\nabla u(x)\right) + a(x)u(x) = f(x)$

The particular problem in question is an electrostatic problem, but has the solution (voltage) specified at several points inside the domain. Since the PDE tool is setup to apply boundary conditions along edges, I thought I might apply these point conditions by the following process:

For a boundary condition at point $x_0$, where we want to specify $u(x_0) = b_0$,

$c(x_0) = 0$

$a(x_0) = 1$

$f(x_0) = b_0$

Such that at the point $x = x_0$, the above differential equation simply becomes $u(x) = b_0$.

However, this does not appear to result in the correct solution. Is this a mathematically valid process? If not, is there a preferred method of setting the solution values at certain arbitrary points in the PDE tool? I am not using the GUI.

Edit

In my understanding, this is a fairly straightforward process using the FEM method. One merely specifies that the corresponding row of the stiffness matrix [K] is entirely 0, with the exception of a 1 at the index of the corresponding node. Then the value of that point is placed in the corresponding row of the solution [b] matrix. Unfortunately, these matrices are not accessible in the PDE tool, hence the question of how to proceed.

Resolution

Thank you all for your help. In the end, it appears that prescribed value points are not currently compatible with MATLAB's adaptmesh algorithm. I ended up writing my own mesh adaption algorithm and using the assempde commands to specify the matrix data and solve.

Answer 1

As said in the comments, it is not mathematically appropriate to constrain a single point value in the given BVP.

The simplest way to achieve something similar to this is to create small circular punctures into your domain and set the prescribed value at the boundary of those.

With small enough circles the solution should be very close to the solution given by the Lagrange multiplier method proposed in the comments. Furthermore, it is probably easier for you to apply adaptmesh with this approach.

If you wish to use the Lagrange multiplier method, you could implement the refinement iteration by yourself, compute error indicators and apply the function refinemesh to the mesh iteratively.

Answer 2

Here is the Lagrange multiplier approach alluded to by Christian Clason.

Structurally, I hope you agree that your problem can be put into the form, \begin{align} \text{argmin}_{u}&\quad \frac{1}{2}||Au - f||^2 \\ \text{such that}&\quad Bu=g, \end{align}

where

• $A$ is the PDE operator,
• $f$ is the PDE right hand side,
• $B$ is the "observation operator" that takes $u$ and observes it, returning the vector of values of the function at the points of interest, and
• $g$ is the list of values that the function $u$ is supposed to take at those points.

For elliptic problems like this, $A$ is symmetric and positive, so it has a square root. Thus one could instead replace the objective function with, \begin{align} \frac{1}{2}||Au-f||^2 &\rightarrow \frac{1}{2}||A^{1/2}u-A^{-1/2}f||^2 \\ &= \frac{1}{2} u^*Au - u^*f + \frac{1}{2}f^*A^{-1}f \\ \end{align} But since $f^*A^{-1}f$ is a constant value independent of $u$, it doesn't effect the optimal point. Thus we can instead consider the equivalent optimization problem,

\begin{align} \text{argmin}_{u}&\quad \frac{1}{2}u^* A u - u^*f \\ \text{such that}&\quad Bu=g, \end{align} The Lagrangian for this problem is, $$L = \frac{1}{2}u^* A u - u^*f + \lambda^*(Bu-g).$$ At the optimal point, the gradient of the Lagrangian must be zero, from which we get the following linear system to be solved, $$0 = \begin{bmatrix}\nabla_u L \\ \nabla_\lambda L\end{bmatrix} = \begin{bmatrix}Au-f + B^*\lambda \\ Bu - g\end{bmatrix},$$ which is the same as the following "KKT system", $$\begin{bmatrix}A & B^* \\ B\end{bmatrix}\begin{bmatrix}u \\ \lambda\end{bmatrix} = \begin{bmatrix}f \\ g\end{bmatrix}.$$

This KKT system can be solved directly, or it can be solved iteratively with Krylov methods that work for indefinite systems (e.g. MINRES). If you choose to go the iterative route, there is some well-established theory by Murphy, Golub, and Wathen that shows the following preconditioner, $$P = \begin{bmatrix}A \\ & B A^{-1}B^*\end{bmatrix}$$

clusters the eigenvalues of the KKT matrix onto precisely 3 points (I.e., $\lambda(P^{-1} \text{KKT})$ is a set with 3 elements), and so using this preconditioner, Krylov methods like MINRES will converge in at most 3 iterations. You can then replace $A$ with a preconditioner $\tilde{A}$ (e.g., multigrid). Furthermore, the rank of $B A^{-1}B^*$ is simply the number of points of interest, so you can replace it with the identity and use the preconditioner $$\tilde{P} = \begin{bmatrix} \tilde{A} \\ & I \end{bmatrix}$$ at the cost of taking a small number of extra number of Krylov iterations.

Now, in terms of discretizations of the problem with finite element methods, you will have to replace operators with their discrete versions, and put mass matrices in appropriate places, yielding the discrete KKT system, $$\begin{bmatrix}K & B^T \\ B\end{bmatrix}\begin{bmatrix}u_h \\ \lambda\end{bmatrix} = \begin{bmatrix}M f_h \\ g\end{bmatrix},$$ and $$\tilde{P} = \begin{bmatrix}\tilde{K} \\ & I\end{bmatrix}.$$ where $K,M$ are the stiffness and mass matrices, $u_h,f_h$ are the discrete versions of $u,f$ (e.g., vector of nodal values), $\tilde{K}$ is the discrete preconditioner, and $B$ is the observation operator that works on the discrete space. E.g., the action of each row of $B$ is to take some weighted linear combination of the nodal values for points in the mesh surrounding the observation node.