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?

  • $\begingroup$ 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. $\endgroup$
    – Kyle Kanos
    Jan 19, 2017 at 22:33
  • $\begingroup$ You're right. How do I transfer it? Or should I delete this one and post in the other? $\endgroup$
    – alexvas
    Jan 19, 2017 at 22:55
  • $\begingroup$ You have to Flag is for moderator attention & request migration there $\endgroup$
    – Kyle Kanos
    Jan 19, 2017 at 22:58
  • $\begingroup$ My take on the question is in lectures 21.6 and 21.65 here: math.colostate.edu/~bangerth/videos.html $\endgroup$ Jan 20, 2017 at 19:52

1 Answer 1


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


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.

enter image description here

  • $\begingroup$ Thanks but this doesn't address the third condition: the discontinuity in the first derivative at some $x_0\in [0,L]$ $\endgroup$
    – alexvas
    Jan 20, 2017 at 20:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.