Impose PDE itself as Boundary Condition?

Consider, for example, the elliptic PDE $u_{x}+u_y+u_{xx}+u_{yy}=0$ for $(x,y)\in[0,\infty)^2$.

Solution methods often require me to impose boundary conditions. Often, these arise naturally from applications (physics, biology, economics etc.). But sometimes, one may not know ex ante how the function should behave on a particular boundary.

Suppose my modeling setting imposes a behavior on $u$ as $x\to\infty$, $y=0$ and $y\to\infty$. However, my application in mind does not give me any priors on what $u$ does as $x=0$. Nonetheless, I have to impose something on this boundary. Can I simply impose that, as $x\to0$, the following inhomogeneous ODE holds $$u_x(0,y)+u_y(0,y)+u_{xx}(0,y)+u_{yy}(0,y)=0.$$ Because we are on the lower boundary of the $x$-domain, we can impose this conditions via finite (forward) differences, \begin{align*} \frac{u_{1,j}-u_{0,j}}{\Delta x}+\frac{u_{0,j+1}-u_{0,j-1}}{2\Delta y}+\frac{u_{2,j}-2u_{1,j}+u_{0,j}}{(\Delta x)^2}+\frac{u_{0,j+1}-2u_{0,j}+u_{0,j-1}}{(\Delta y)^2}=0, \end{align*} where $u_{i,j}=u(x_i,y_j)$.

So, instead of coming up with some Dirichlet or Neumann conditions, I essentially just impose the PDE itself on one (or several) of the boundaries. Is this something one has heard of? Do you know of problems with similar boundary conditions?

Solution 1:

The underlying problem is that the boundary condition you have provided does not assist in evaluating the solution.

Consider that, for your PDE, if $u=f(x,y)$ is a solution, then so is $u=f(x,y)+C$ for some constant, $C$. This constant cannot be determined using your boundary condition.

Indeed, consider that the following expression is a solution for any choice of constants $C_i$: $$ u=C_1+C_2(x-y)+C_3e^{-x}+C_4e^{-y}+C_5e^{-x-y} $$ None of these constants can be determined by the boundary condition as provided.

That being said, there are often implied conditions that can be used to obtain boundary conditions. For example, if you are working in a real-world context, you may need to apply an energy-conservation requirement. Alternatively, boundedness and related concepts often imply some form of boundary condition. For instance, consider $$ xy'+y=0 $$ In the limit as $x\to0$, we get $y=0$ unless $y'$ tends to infinity. Indeed, the general solution is $y=C/x$, and so if $y$ must remain bounded, then $y(0)=0$, and thus $y=0$.

However, this is actually a separate condition on the problem that generates a boundary condition. You cannot simply apply the DE to the boundary, as it provides no additional information, and boundary conditions are applied to provide information that is otherwise missing.

The mistake you make is the claim that you must apply something at the boundary. This is untrue. You must apply a condition that allows you to identify the correct solution (unless you want a family of solutions, that is).

The nature of the condition will depend on the goal of the differential equation. As mentioned above, sometimes you need to apply some energy-conservation condition. Other conditions are possible - minimisation conditions are common. For example, you might want the solution to your problem that has the least amount of variation (that is, minimise $\iint u_{x}^2+u_{y}^2\ dA$).

There is one other category - cases where your choice won't matter. This is much like how, when performing integration by parts, you don't need to include the constant of integration for the "$\int dv$" - it will cancel out anyway. In this situation, pick the simplest boundary condition for the problem - for instance, if $u(0,0)=1$ and $u(0,\infty)=0$, then you might pick $u(0,y)=e^{-y}$.

EDIT: Here, I'll show what your forward-difference application of the PDE at the boundary is actually assuming for the boundary. Assuming you use central difference for first derivatives when not at the boundary, you have... $$ \begin{align*} \frac{u_{1,j}-u_{0,j}}{\Delta x}+\frac{u_{0,j+1}-u_{0,j-1}}{2\Delta y}+\frac{u_{2,j}-2u_{1,j}+u_{0,j}}{(\Delta x)^2}+\frac{u_{0,j+1}-2u_{0,j}+u_{0,j-1}}{(\Delta y)^2}&=0\\ \frac{u_{2,j}-u_{0,j}}{2\Delta x}+\frac{u_{1,j+1}-u_{1,j-1}}{2\Delta y}+\frac{u_{2,j}-2u_{1,j}+u_{0,j}}{(\Delta x)^2}+\frac{u_{1,j+1}-2u_{1,j}+u_{1,j-1}}{(\Delta y)^2}&=0 \end{align*} $$ With a bit of algebra, we can cancel the $u_{2,j}$ terms to get $$\begin{align*} \left(\frac12+\frac{2+\Delta x}{(\Delta y)^2}\right)u_{0,j} - \left(\frac12+\frac2{(\Delta y)^2}\right)u_{1,j} & \\ + \left(\frac{\Delta x}{4\Delta y} - \frac{\Delta x}{2(\Delta y)^2} + \frac{1}{2\Delta y} - \frac{1}{(\Delta y)^2}\right)u_{0,j-1} - \left(\frac{1}{2\Delta y} - \frac{1}{(\Delta y)^2}\right)u_{1,j-1} & \\ - \left(\frac{\Delta x}{4\Delta y} + \frac{\Delta x}{2(\Delta y)^2} + \frac{1}{2\Delta y} + \frac{1}{(\Delta y)^2}\right)u_{0,j+1} + \left(\frac{1}{2\Delta y} + \frac{1}{(\Delta y)^2}\right)u_{1,j+1} & =0 \end{align*}$$ By Taylor expansion around the $(0,j)$ point, letting $u$ be the function at that point and stopping at $O(\Delta)$ (treating both $\Delta x$ and $\Delta y$ as proportional to $\Delta$), we get $$\begin{align*} \frac{\Delta x}{(\Delta y)^2} u - \frac{\Delta x}2 (u_x+u_y+2u_{xy}) - \frac{\Delta x^2}4 u_{xx} + O(\Delta) & =0 \end{align*}$$ where I have left in the terms that come out explicitly... but note that this is equivalent to $$ \frac{\Delta x}{(\Delta y)^2} u + O(\Delta) = 0 $$ And so, if $\Delta x$ and $\Delta y$ are both scaled down in proportion with each other, your boundary condition is effectively $u(0,y)=0$.

However, be careful, as this analysis applies for keeping $\Delta x$ and $\Delta y$ in proportion. I suspect that some issues may arise if you reduce $\Delta x$ without reducing $\Delta y$.

Solution 2:

Since you're using a finite difference method, you're turning the the PDE into a set of algebraic equations, I'm sure you know this.

It seems to me that you won't ever get a solution because you are just adding an equation (expanding the dimensions of the matrix $A$ in the equation $Ax = b$). So your issue is that you don't have the first or last entries of the $b$ vector in the above equation. So you're just finding the kernel of the system, but if the system is non-singular (which I think it is) then the kernel is just zero.

So just a few ideas as far as options go:

since you're working in a semi-infinite domain, maybe try other methods of solution like similarity or characteristics.

If you need to use a finite difference approach, then maybe use some scaling arguments. Does it really matter what value $u$ takes on at that boundary, or is it arbitrary and the relevant physical phenomena is really happening far away from that boundary. Maybe try a bunch of values for $u$ or $\frac{\partial u}{\partial n}$ at the boundary. It looks like you're in luck since your governing equation is linear, so there should be some pattern established when you vary the prescribed boundary condition.

Hope this helps!