# Numerical solution to 2D divergence equation

Is there any way to numerically solve the following two-dimensional equation: $$\nabla_{xy} \cdot \vec{f}(x,y) = a(x,y)$$ on a rectangular grid, knowing that $\vec{f}(x,y)$ is zero on the boundaries?

The problem you have is not well-posed. In essence, you are looking for two functions $f_1(x,y),f_2(x,y)$ that should solve only one equation. That's not enough information.

To give you an idea why this can't work, assume that you had found one set of solutions $\vec f(x,y)=(f_1(x,y),f_2(x,y))^T$. Then it is not difficult to verify that that $\tilde {\vec f}(x,y) =(f_1(x,y)+g(y),f_2(x,y)+h(x))^T$ is also a solution for any functions $g(y),h(x)$ you may want to choose. With one stroke of a pen, I have therefore already found an infinity of solutions -- too many to do anything useful with.

In other words, it is a moot point whether there are numerical methods that can solve your problem: There are fundamental mathematical problems you first need to address, and then we can talk about numerical approximation.

(You may want to object that the boundary conditions require me to choose $g(y)=h(x)=0$. But that only shifts the problem: I can still construct a $\tilde{\vec f}(x,y)=\vec f(x,y) + \vec \varphi(x,y)$ for any $\vec\varphi$ that satisfies $\text{div}\vec\varphi=0$ and $\vec\varphi=0$ on the boundary. There are infinitely many such functions, for example all solutions of the Stokes equation with zero boundary values and arbitrary right hand sides.)

There are many ways to do this, either finite difference methods, finite volume methods, or finite element methods can be applied.

The above equation often known as the poison equation is often the first presented equation in a new topic in a numerical methods course.

• But if I remember correctly, Poisson equation is of the type: $$\nabla \cdot \nabla h(x,y) = a(x,y)$$ in which we are using (for example) the finite volume method to solve for a scalar function, and not a vector function. – Maziar Noei Nov 14 '16 at 14:47
• Sorry I misread – Intrepid Traveller Nov 14 '16 at 14:48
• My above comment still applies. However this divergence equation might need stabilization terms added to the standard methods listed above. – Intrepid Traveller Nov 14 '16 at 14:49
• Could you be more specific or provide a reference? I was thinking it cannot be solved because the number of unknowns exceeds the number of equations, but I am not sure about this. I am interested in as much information as we can get about $\vec{f}(x,y)$. – Maziar Noei Nov 14 '16 at 14:52

I have constructed the matrix associated with the finite difference discretisation of $u_{xy}$ in an answer to a stack overflow question Matrix to generate finite difference (the matrix is in the MATLAB code as A).

Alternatively, mentioned on wikipedia and referencing Evans (2010) (p. 400), you can transform $u_{xy}$ equation into the one dimensional wave equation. I've attempted this below. Boundary conditions could get trickyer.

If you want to get rid of your mixed derivatives you can use a transformation of coordinates

\begin{aligned} \hat{x} &= x+y\,, & x = \frac{\hat{x}+\hat{y}}2\,, \\ \hat{y} &= x-y\,, & y = \frac{\hat{x}-\hat{y}}2\,. \end{aligned} You can show that $$\begin{gathered} u_{\hat{x}\hat{x}} = u_{xx} + 2 u_{xy} + u_{yy}\,,\\ u_{\hat{y}\hat{y}} = u_{xx} - 2 u_{xy} + u_{yy}\,, \end{gathered}$$ then \begin{aligned} u_{xy} {}&{}= \frac{ u_{\hat{x}\hat{x}} - u_{\hat{y}\hat{y}} }4 \\ {}&{}= a\left(\frac{\hat{x}+\hat{y}}2,\frac{\hat{x}-\hat{y}}2\right)\,. \end{aligned}

• Thanks Steve, but as I mentioned, the equation is a divergence equation. That is, $$\vec{f} = f_x \hat{e}_x + f_y \hat{e}_y$$ $$\nabla_{xy}\cdot \vec{f} = \dfrac{\partial f_x(x,y)}{\partial x} + \dfrac{\partial f_y(x,y)}{\partial y} = a(x,y)$$ with said boundary conditions. – Maziar Noei Nov 14 '16 at 18:13

As Wolfgang explained, the issue is not to find a solution, but to pick one.

If just some solution is fine for you, go ahead and take $f_2\equiv 0$ (for example) and solve $\partial_x f_1(x,y) = a(x,y)$ for every $y$ (which means to integrate up, e.g. $f_1(x,y) = \int_{x_0}^x a(t,y)dt$). But if you are not happy with this solution for some reason, you need to think about what other conditions you have on a solution.

Just two examples:

1. You may want to look for solutions with minimal $L^2$-norm. This would lead to a quadratic optimization problem with linear equality constraints, namely $$\min \int f_1(x,y)^2 + f_2(x,y)^2 dxdy\quad\text{s.t.}\quad \operatorname{div} f = a.$$ The optimality system for this question is another system of PDEs…

2. You may be interested in a solution that is curl-free. Then make the ansatz $f = \nabla u$ (which implies $\operatorname{curl}f = \operatorname{curl}\nabla u = 0$) and solve $$\operatorname{div}\nabla u = \Delta u = a$$ which is indeed a Poisson equation.

• The direct integration suggested at the top does not work in general because it leads to a function that does not satisfy the boundary conditions. – Wolfgang Bangerth Nov 16 '16 at 0:19
• Option 1 likely also does not work because I'm pretty sure one can construct sequences of functions (all of which satisfy the equation) whose $L_2$ norm goes to zero, but whose limit does not satisfy the equation. A more useful minimization function would minimize the $H^1$ seminorm, for example. – Wolfgang Bangerth Nov 16 '16 at 0:21
• Option 2 is complicated by the fact that you need to satisfy boundary conditions on $\nabla u$ all along the boundary, but the solution of the Poisson equation only allows you to impose boundary conditions on the normal component of $\nabla u$. – Wolfgang Bangerth Nov 16 '16 at 0:22
• When I wrote about the integration I hadn't noticed the boundary conditions. Also the complications with the boundary conditions in option 3 is indeed real. I am not sure about the problem with 2 though. – Dirk Nov 16 '16 at 1:13
• These were all good ideas, though! – Wolfgang Bangerth Nov 16 '16 at 3:48