# Infinite Function Value on Dirichlet Boundary

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 )$

• Can you invert the problem such that the Direchlet BC decays to 0 at the boundaries? If so then choosing the Laguerre polynomials and using a spectral method could be of use. – Aurelius Oct 31 '14 at 18:25
• What's the underlying PDE? – Bill Barth Oct 31 '14 at 19:01
• Your question isn't particularly clear -- do your exact boundary values go to zero, or is it your numerical solution that (erroneously) goes to infinity even though it shouldn't? – Wolfgang Bangerth Nov 1 '14 at 11:48
• @Aurelius: Can you explain what you mean by 'invert the problem'? – pdevar Nov 1 '14 at 14:34
• I still don't quite understand. Are you posing Dirichlet boundary values of the form $h^i(x)|_{\Gamma}=h_\text{Dirichlet}(x)$ and similar for $f^i$ where $h_\text{Dirichlet}(x)$ is infinite at one point? – Wolfgang Bangerth Nov 2 '14 at 0:43

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.

• I'm actually using multigrid methods and have yet to try finite element methods. As far as I can tell, they are very similar, only differing in the grid imposed. About the well-posedness of the problem, you may be right. I am going to spend some more time studying the boundary and the limits of the underlying processes. – pdevar Nov 3 '14 at 8:18

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.

• Thanks Wolfgang, I will have to take some time to think about this solution a bit more, as it is a bit over my head. I had another idea to truncate the state space and consider a new problem with mixed Von Neumann and Dirichlet conditions, something like \begin{align*} -\Delta u &= f &in \hspace{2pt} \Omega \\ u &= G &on \hspace{2pt} \partial\Omega \end{align*} becomes \begin{align*} -\Delta u &= f &in \hspace{2pt} \hat{\Omega} \\ -\Delta u &= f &on \hspace{2pt} \partial\hat{\Omega}\\ u &= g &on \hspace{2pt} \partial \Omega_{homog} \end{align*} where $\hat{\Omega} \subset \Omega$. – pdevar Nov 5 '14 at 15:39
• That is terrible notation, but I'm in a bit of a hurry right now... I know that the functions equal a finite and small value, $g$, on one of the boundaries and so should be able to use this information to calculate the solution over the trunctated state space by defining the Von Neumann condition as the function itself on the new boundary... right? – pdevar Nov 5 '14 at 15:46
• I must admit that I don't follow what exactly you are trying to do, but basically the approach needs to be: replace the original problem by a different one that doesn't suffer from the pathology you have, and try to show that the error you introduce with this substitution is of the same order as the discretization error. If you can show this for your approach, then you're good. – Wolfgang Bangerth Nov 5 '14 at 19:08