Infinite Function Value on Dirichlet Boundary

Question

I have been working on a multigrid solution to a non-homogeneous Dirichlet boundary value problem. However, the function goes to infinity on the boundary. This causes numerical overflow errors to be raised in calculating the discretized derivatives. Does anyone have a suggestion on how to overcome this problem? I am not a mathemetician, but the Galerkin Method seems to be a way to transform the problem, but I am not sure.

Here are the pde's I'm solving (I've not included the formulae for the drift and diffusion coefficients as they are also long):

$\frac{1}{2} \omega_i^2 \sigma_{\omega_i}^2 \frac{\partial^2 f^i}{\partial \omega_i^2} + \frac{1}{2} \delta^2 \sigma_\delta^2 \frac{\partial^2 f^i}{\partial \delta^2} + \delta \omega_i \sigma_{\omega_i} \sigma_\delta \frac{\partial^2 f^i}{\partial \delta \partial \omega_i} + \delta \mu_\delta \frac{\partial f^i}{\partial \delta}$

$+ \omega_i \mu_{\omega_i} \frac{\partial f^i}{\omega_i} - \rho f^i + (\omega^i)^{-\gamma_i} \left ( tan \left [ \pi(\delta - \frac{1}{2}) \right ] \right )^{1-\gamma_i}= 0$

$\frac{1}{2} \omega_i^2 \sigma_{\omega_i}^2 \frac{\partial^2 h^i}{\partial \omega_i^2} + \frac{1}{2} \delta^2 \sigma_\delta^2 \frac{\partial^2 h^i}{\partial \delta^2} + \delta \omega_i \sigma_{\omega_i} \sigma_\delta \frac{\partial^2 h^i}{\partial \delta \partial \omega_i} + \delta \mu_\delta \frac{\partial h^i}{\partial \delta}$

$+ \omega_i \mu_{\omega_i} \frac{\partial h^i}{\omega_i} - \rho h^i + \omega^i(t)^{-\gamma_i} \left ( tan \left [ \pi(\delta(t) - \frac{1}{2}) \right ] \right )^{-\gamma_i} = 0$

where
$H_0(t) = exp \left ( -\int_0^t r(u) du \right) exp \left (- \int_0^t \frac{\sigma_D}{\xi(u)} du - \frac{1}{2} \int_0^t \left ( \frac{\sigma_D}{\xi(u)} \right)^2 du \right )$

Answer 1

Most of the standard methods for solving PDE numerically are not suited to handle infinite values, since they originate from physical problems. Infinite values of the solved variables may occur in physical problems solved numerically, for example, if one tries to solve the problem of a hairline crack in a fully elastic medium using the stress-displacement method (infinite stresses develop at the crack tip, which is physically wrong but does not affect the solution far from the tip).

To deal with this situation through the Finite Element method, one option would be to enrich the basis of the FE discretization, for example as done in the XFEM method, see this site, where in the present case the global enrichment function will be an exponent tending to infinity as one approaches your boundary.

I should note that it is not exactly clear what a Dirichlet boundary on which the solution is infinite actually means. It seems to be a badly-posed problem. Perhaps the rate at which the solution approaches infinity on that boundary would be a better candidate for this boundary condition? That would be ideal, as it would dictate the exponent of the enrichment function.

Answer 2

I don't know off the top of my head how to deal with the situation you describe, but let me outline the general approach one would take.

Let's pretend for a minute that you were solving the following, simpler problem: $$ -\Delta u = f \qquad\qquad \text{in $\Omega$} \\ \qquad u = 0 \qquad\qquad \text{on $\partial\Omega$} $$ and assume that $f$ has a singularity at a point $x_0$, i.e., it goes to infinity at this point. Now imagine that you also consider the following, perturbed problem $$ -\Delta \tilde u = \tilde f \qquad\qquad \text{in $\Omega$} \\ \qquad \tilde u = 0 \qquad\qquad \text{on $\partial\Omega$} $$ where $\tilde f$ is a perturbed version of $f$. Then $$ \| u -\tilde u \|_{H^1} \le \| f-\tilde f\|_{H^{-1}} $$ base on elementary a priori stability analysis. So if you discretize the second version using finite elements of finite differences (using for simplicity a lowest order method), you get an approximate solution $\tilde u_h$ that satisfies the estimate $$ \| u -\tilde u_h \|_{H^1} \le \| u -\tilde u \|_{H^1} + \| \tilde u - \tilde u_h \|_{H^1} \le \| f-\tilde f\|_{H^{-1}} + C h \| \tilde u \|_{H^2}. $$ Here, $h$ is the me size. In other words, if you solve a version of the equation perturbed in such a way that $$ \| f-\tilde f\|_{H^{-1}} \le Ch $$ then you get a method that has the correct convergence order. It's too late in the night right now to prove this with any certainty, but my best guess is that if you choose $\tilde f=f \ast G_h$ where $G_h$ is a Gaussian with width $h$, then you get at least this order. Note that this smoothed version of $f$ is always finite and so can serve as the basis for a discretization.

What you'll have to do is to go through a similar argument for boundary values instead of right hand side values. It shouldn't be too hard to do this. In the end, what you'll have to do is just convince yourself that if you smooth the offending function a bit, that you still converge to the correct solution and at the same convergence order as always.