# Discretization of Laplacian with boundary conditions

I am trying to solve the equation below numerically for various functions $I(x)$:

$$\frac{\partial^2}{\partial x^2} G(x) + I(x)G(x) = 0$$

Subject to the boundary contitions:

$$\frac{\partial}{\partial x} G(x)\Big|_{x=0} = C_1$$ $$\frac{\partial}{\partial x} G(x)\Big|_{x=L} = C_2$$ $$\frac{\partial}{\partial x} G(x)\Big|_{x=x_0+\epsilon} - \frac{\partial}{\partial x} G(x)\Big|_{x=x_0-\epsilon} = 1$$

Where $x_0 \in [0,L]$, $\epsilon\to 0$, and $C_1$ and $C_2$ are given.

Normally, if there were no boundary conditions, I would discretize $[0,L]$ into $N$ individual units, then I would approximate $\frac{\partial^2}{\partial x^2}$ using the finite difference method yielding a $N\times N$ matrix. I would then add the discretized version of $I(x)$ as a matrix by multiplying by the identity $N\times N$ matrix, yielding:

$$\left( \begin{bmatrix} -1 & 1 & 0 & \dots & 0 & 0 \\ 1 & -2 & 1 & \dots & 0 & 0 \\ \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & -2 & 1 \\ 0 & 0 & 0 & \dots & 1 & -1 \end{bmatrix} + \mathbb{1} \cdot \begin{bmatrix} I(x_1) \\ I(x_2) \\ \vdots \\ I(x_N) \\ \end{bmatrix} \right) \begin{bmatrix} G_0 \\ G_1 \\ \vdots \\ G_{N-1} \\ G_N \\ \end{bmatrix} = 0$$

Then finding $G$ is a matter of inverting the matrix on the left hand side.

My question: how do I incorporate the boundary conditions into he expression above? Should the derivative matrix (first matrix) change, or should the right hand side become a vector other than $0$? Is there a good reference for learning about how to numerically solve PDEs with weird boundary conditions?

• Personally, I would apply it separately rather than in the matrix. However, this isn't much of a physics question. It might be more appropriate for Mathematics, Computational Science or Stack Overflow. – Kyle Kanos Jan 19 '17 at 22:33
• You're right. How do I transfer it? Or should I delete this one and post in the other? – alexvas Jan 19 '17 at 22:55
• You have to Flag is for moderator attention & request migration there – Kyle Kanos Jan 19 '17 at 22:58
• My take on the question is in lectures 21.6 and 21.65 here: math.colostate.edu/~bangerth/videos.html – Wolfgang Bangerth Jan 20 '17 at 19:52

You are correct that you should augment the matrix to satisfy the boundary conditions.

I do not have the book infront of me, but I recall that in the first chapter of the book Fundamental Algorithms in Computational Fluid Dynamics by Thomas H. Pullman and David W. Zingg they discuss in detail the matrix formulation of finite difference methods. In particular they discuss how different boundary conditions change the form of the matrix.

In my experience I have always used periodic boundary conditions. In this case, when you write down the difference equations at the left boundary of domain, you will get a difference approximation of your differential operator of the form

$u_{N-1}-2u_{i}+u_{i+1}$

a similar form can be found for the right boundary of the domain. What this does it it changes the form of your matrix slightly so that it becomes somewhat like the matrix below up to some constants. • Thanks but this doesn't address the third condition: the discontinuity in the first derivative at some $x_0\in [0,L]$ – alexvas Jan 20 '17 at 20:10