I'm very new to finite difference method and I am just introduced to methods of solving differential equation using finite difference method via sparse matrix method.

I find that the main idea is to approximate the differentials by central differences and writing the differential equation at different points in domain and making sparse matrices out of them.

For example, take the following 1-D problem:

$\frac{d^2 f(x)}{dx^2} = f(x)$,

with boundary conditions, $f(0)=1$, $f(2) = e^2$

At different points in domain with h(minimum step size) = 0.25,

$f(0) = 1$

$\frac{f(0) - 2f(0.25) + f(0.5)}{h^2} - f(0.25) = 0 $

$\frac{f(0.25) - 2f(0.5) + f(0.75)}{h^2} - f(0.5) = 0 $

$\frac{f(0.5) - 2f(0.75) + f(1)}{h^2} - f(0.75) = 0 $


$f(2) = e^2 $

Rewriting the above set of equations in sparse matrices: $a\cdot f = b$,

$ \begin{pmatrix} 1 & 0 & 0 & 0 & \cdots & 0\\ 1/h^2 & (-2/h^2 -1) & 1/h^2 & 0 & \cdots & 0\\ 0 & 1/h^2 & (-2/h^2 -1) & 1/h^2 & \cdots & \vdots\\ 0 & 0 & 1/h^2 & (-2/h^2 -1) & \cdots &\\ \vdots & \vdots & \vdots & \vdots &\ddots & 1/h^2\\ 0 &0 &\cdots & & & 1\\ \end{pmatrix}$ $ \begin{pmatrix} f(0)\\ f(0.25)\\ f(0.5)\\ \vdots\\ f(2)\\ \end{pmatrix}$ $=$ $\begin{pmatrix} 1\\ 0\\ 0\\ \vdots\\ e^2\\ \end{pmatrix}$

and then the solution becomes, $f = a^{-1} \cdot b $.

Even for nonlinear differential equations, the method takes the similar footing and these sparse matrices are to be built first before going through iteration methods for solving nonlinear systems.

Now, the question that I have now is, what about differential equations that involves sin and cos of the unknown function?

for example: $\frac{d^2 f(x)}{dx^2} + sin( f(x) ) = 0 $

or system that involves inverse?

for example: $\frac{d^2 f(x)}{dx^2} + \frac{1}{1+f(x)} f(x) = 0 $

how do we tackle problems like these using finite difference sparse matrix methods?

  • $\begingroup$ Note that you can inexpensively increase the accuracy by changing the method to Numerov's method. The order increases from 2 to 4. The idea is that you approximate the leading error term by the same second order finite difference formula, now applied to second derivatives, and replace the second derivatives by inserting the ODE. $\endgroup$ May 10, 2023 at 8:43
  • 2
    $\begingroup$ A better title would be something like "System solver based on finite differences for non-linear second order BVP". Note also the other usual approaches to BVP, single and multiple shooting, collocation methods. $\endgroup$ May 10, 2023 at 8:46

1 Answer 1


A common method is to spatially discretize your ODE/PDE, then use a Newton-Raphson method to solve the resulting non-linear system.

For example, let's start by replacing the second derivative with a 3 point finite difference as you've done, while we evaluate the remaining terms using a "one point stencil". Here I've written it in a slightly more compact form using a sub-index on $f$.

$$ f_0 = f(0)\\ f_1 = f(h)\\ f_2 = f(2h)\\ \vdots\\ f_n = f(2)\\ \frac{f_{i-1} - 2 f_{i} + f_{i+1}}{h^2} + \frac{1}{1+f_i} f_i = 0 $$ Note that the last equation is actually a whole series of equations for $i \in [1, n-1]$. We will need to provide boundary conditions for $f_0$ and $f_n$. We then have the coupled nonlinear system of equations $$ \vec{F}(\vec{f}) = \vec{0} = \begin{bmatrix} f_0 - f(0)\\ f_n - f(2)\\ \frac{f_{i-1} - 2 f_{i} + f_{i+1}}{h^2} + \frac{1}{1+f_i} f_i \end{bmatrix} $$ To solve a system of nonlinear equations, we can first Taylor expand $\vec{F}$ with respect to each component of $\vec{f}$, keeping just the first term: $$ \vec{F}(\vec{f}_0 + \Delta \vec{f}) = \vec{F}(\vec{f}_0) + \nabla_{\vec{f}} \vec{F}(\vec{f}_0) \Delta \vec{f} + \underbrace{\ldots}_{\mathrm{dropped~terms}} $$ where $\vec{f}_0$ here is our initial guess for the solution, and $$ \nabla_{\vec{f}} \vec{F} = \begin{bmatrix} \partial_{f_0} F_0 & \partial_{f_1} F_0 & \partial_{f_2} F_0 & \ldots & \partial_{f_n} F_0\\ \partial_{f_0} F_1 & \partial_{f_1} F_1 & \partial_{f_2} F_1 & \ldots & \partial_{f_n} F_1\\ \vdots\\ \partial_{f_0} F_n & \partial_{f_1} F_n & \partial_{f_2} F_n & \ldots & \partial_{f_n} F_n \end{bmatrix} $$ is the Jacobian matrix. We would like to find an update $\Delta \vec{f}$ which makes $\vec{F}(\vec{f} + \Delta \vec{f}) = \vec{0}$. So let's plug that in, and solve: $$ \nabla_{\vec{f}} \vec{F}_{\vec{f}_0} \Delta \vec{f} = -\vec{F}_0 $$ Notice that the Jacobian matrix will end up being a sparse matrix, but for nonlinear problems some terms are now a function of $\vec{f}$ itself, explicitly evaluated with our guess $\vec{f}_0$. If the system was linear, there would be no dependency on $\vec{f}$ and we can then exactly find the correct update $\Delta f$ to get $\vec{F}(\vec{f}) = 0$ in one sparse linear solve since our Taylor expansion would be exact.

For nonlinear problems the situation is slightly more complicated. If $\vec{f}_0 + \Delta \vec{f}$ is closer to a solution, then we can just repeat the process again until we're sufficiently close to the actual solution. Even if the first few updates actually takes us away from the true solution, sometimes it will still eventually end up converging to a solution. However, convergence to a solution is not guaranteed unless our initial guess is "sufficiently close", for a problem dependent definition of sufficiently close. Some problems almost any initial guess is sufficiently close. For others, only guesses which are exact roots are "sufficiently close" (in which case the Newton-Raphson method isn't useful anymore).

The Newton-Raphson method can be said to converge locally, but not globally. There are ways to extend its convergence using line searches or trust regions, though I'll leave it at this for now since the base Newton-Raphson method does work a lot of the times.


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.