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.


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.


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.

  • $\begingroup$ Welcome to SciComp.SE! This is not the mathematically correct way to incorporate constraints; if $x_0$ is a point on the boundary, you need to define a sub-boundary that only consists of this point (and you should ask on Matlab Answers how to do this with the PDE toolbox; note that there were substantial changes in the most recent version). If $x_0$ is in the interior of the domain, this is not a boundary condition at all, and needs to be enforced via a Lagrange multiplier. $\endgroup$ Oct 23 '15 at 15:52
  • $\begingroup$ Thanks Christian. My use of the term 'boundary condition' is indeed not very appropriate here, I suppose I am referring to specified solutions at given given points. I've had limited success with Matlab Answers, so I thought I'd try here :) Thanks, I'll look into Lagrange multipliers! $\endgroup$ Oct 23 '15 at 16:24
  • $\begingroup$ Sorry, should have read the question more carefully. My suggestion then would be to edit the question (and the title) to reflect that, and also make the question more about the mathematical approach (and add a note at the end that you're particularly interested in implementing this using the PDE toolbox). This way, you get useful answers from people who know nothing about the toolbox, and possibly hints from those who do. (There are many more of the former than the latter here, so it might make sense to later ask more concrete questions about the implementation at Matlab Answers.) $\endgroup$ Oct 23 '15 at 16:44
  • $\begingroup$ Another piece of important information would be whether you can assume that $x_0$ is a vertex of the triangulation. Also, could you add some details about the physical background? Mathematically, it doesn't make sense to prescribe the solution at a point, since the PDE has a unique solution -- which either has the prescribed values or not. So something has to give, and that would be either the PDE or the right-hand side in this case, which might not be what you want. $\endgroup$ Oct 23 '15 at 16:51
  • $\begingroup$ Re your Edit: Yes, they are -- [K,M,F,Q,G,H,R] = assempde(___) outputs the stiffness and mass matrices K and M as well as a matrix H describing the Dirichlet conditions, which you can then modify at will. You can then solve the modified problem using u = assempde(K,M,F,Q,G,H,R). $\endgroup$ Oct 23 '15 at 17:06

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. PDE Toolbox, prescribing values inside the domain

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.

  • $\begingroup$ Thanks knl. I ended up writing my own mesh adaption algorithm, similar to what you suggested, and using the assempde commands to modify the matrices' solution points. $\endgroup$ Nov 6 '15 at 3:43

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}


  • $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.

  • $\begingroup$ Thanks Nick Alger for our detailed response! I ended up writing my own mesh adaption algorithm and using the assempde commands to modify the matrices' solution points -- no Lagrange multipliers necessary! $\endgroup$ Nov 6 '15 at 3:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.