# Numerically find Greens function

I am trying to numerically evaluate a Greens function for this equation:

$$\left[\frac{\partial^2}{\partial x^2} + f(x) \right] G(x) = \delta (x-x_0)$$

With Neumann boundary conditions. Here, the function $f(x)$ is known. The difficulty in solving this lies in the numerical approximation of the $\delta$ function. Rather than approximate $\delta \sim \frac{1}{x}$ I attempted to resolve the $\delta$ by integrating both sides in a small interval around $x_0$, but this leaves me with conditions on the derivative of $G$ at $x_0$ that I don't know how to implement in a finite elements method.

I would really appreciate a reference for solving such equations (preferably with a finite elements method), or if someone could lend a guiding hand for solving such an equation.

• What is your domain? Is it the whole real line? If that's the case, why do you need an $x_0$? – nicoguaro Feb 2 '17 at 21:11
• The domain is bounded: $x\in [0,L]$ with $0<x_0<L$. The conditions of $G$ on the boundary are Neumann: $G'(0)=C_1$ and $G'(L)=C_2$ where both $C_1$ and $C_2$ are known. – alexvas Feb 2 '17 at 21:17
• I suppose that depending on $f(x)$ you can write your operator as a Sturm-Liouville operator, and find the eigenvectors, then express $G(x)$ as an eigenvector expansion. – nicoguaro Feb 2 '17 at 21:29
• I'm not sure I understand. To turn this into an eigenfunction equation I need to have $G$ on both sides. I could move $f(x)G(x)$ over to the RHS, but what do I do with the $\delta$? I don't see how to rewrite this as a SL equation. – alexvas Feb 2 '17 at 21:56
• No, check this. And also this. – nicoguaro Feb 2 '17 at 21:59

Why approximate the delta function? In the finite element method, you need to evaluate the weak form of the equation, for which the right hand side will look like this: $$\int \varphi_i(x) \delta(x-x_0),$$ where $\varphi_i$ is the $i$th shape function. But you know what this integral is: $$\int \varphi_i(x) \delta(x-x_0) = \varphi_i(x_0).$$ Of course, most shape functions will be zero at $x_0$, but the shape functions defined on the cell that contains $x_0$ will not, and those are going to be the only nonzero entries in your right hand side vector.