2

I want to solve the following simple differential equation numerically on a grid:

$$\frac{df(x)}{dx}=b(x)$$

I tried to discretize the LHS using central differences, but the problem with that is, that the resulting linear system of equations $Ax=b$ has a singular matrix $A$ with a zero diagonal.

I also tried forward/backward difference, which gives me a (still singular) matrix $A$ with diagonal elements. I tried the Gauss-Seidel method, but got a wrong solution, because the single-sidedness of the discretization causes information to propagate only in a single direction (best seen when b is a discretized delta distribution).

My question now is: Is there a way to discretize the equation above into a linear system $Ax=b$ which can be solved and gives the correct solution?

Thanks in advance for any help or comments.

Winther
  • 24,478
dkoerner
  • 23
  • 1
    Welcome to Math.SE. This is a numerical integration problem "in disguise". The singularity of your matrix $A$ is expected because you need an "initial condition" to determine $f(x)$. That is the same as the $+C$ issue of indefinite integration in calculus, i.e. if $f(x)$ solves the equation, so will $f(x)+C$. Adding an equation to your system for the initial condition will make the solution unique, everything else being correct. – hardmath Mar 15 '17 at 11:58
  • Thanks for your comment, @hardmath. Adding an equation will make $A$ non-square, which means that there is no guarantee that I get an exact solution, right?

    I dont really understand for what I need to specify inital values and how this turns into an additional equation in my system. Specifying an initial condition for $f(x)$ would mean that I need to specify an additional equation for each element of my grid as far as I understand.

     – dkoerner Mar 15 '17 at 12:19
  • It is not so strange since it is not uniquely solvable (which translates into underdetermined system). You could add for example any constant to $f$ and get another solution. You could add regularization terms or boundary conditions to try and make it uniquely solvable. – mathreadler Mar 15 '17 at 12:23
  • 1
    Getting "an exact solution" is not the issue (you are doing numerical integration, approximation is baked in), but rather the issue is getting a unique solution. Do a smallish example by hand, adding the extra equation and solving the linear system by row reduction. As you've already observed $A$ is singular, you might find a redundant row of $Ax=b$ to discard if you care about keeping the system "square". – hardmath Mar 15 '17 at 12:24
  • 2
    "I need to specify an additional equation for each element of my grid" No, you only need to specify the value of $f(x)$ at one point of the grid, and the solution propagates from there. In effect you are removing one unknown from the system $Ax=b$. – hardmath Mar 15 '17 at 12:27
  • Yes one piece of information, it could be a mean value or a boundary value or maybe something third. – mathreadler Mar 15 '17 at 12:29
  • An alternative to solving it as the ODE is to treat it as an integration problem. We have $f(x) = f(a) + \int_a^{x}b(t){\rm d}t$ where $f(a)$ is an initial condition. Now if you have an array of $x$-values $a = x_0 < x_1 < \ldots <x_n = b$ then $f(x_{i+1}) = f(x_i) + B_i$ where $B_i = \int_{x_i}^{x_{i+1}}b(t){\rm d}t$. This can for example be approximated as (trapezoidal rule): $B_i \simeq (x_{i+1}-x_i)\frac{b(x_i) + b(x_{i+1})}{2}$. You can phrase this as a linear algebra problem with a bidiagonal $A$ if you want, but it's easier to just evaluate $f(x_i)$ directly using a loop. – Winther Mar 15 '17 at 13:28
  • @Winther: Thanks for your comment. I think mathematically this would be identical to specifying a boundary condition on the $a$ side of the domain for the linear system approach (but it would certainly be more efficient).

    The problem for me is, that I have a delta distribution within the domain and therefore can not provide dirichlet-BC. Maybe your approach could be extended. It seems that we then would get something akin to the fast marching method.

     – dkoerner Mar 15 '17 at 13:41
  • If you have a delta-function in $b(x)$ e.g. $b(x) = b_0(x) + C\delta(x-x_)$ then in this approach we simply add $C$ to $B_i$ if $x_ \in (x_i,x_{i+1})$ since $\int_{x_i}^{x_{i+1}}C\delta(t-x_){\rm d}t = C$ if $x_$ is in this interval and $0$ otherwise. This just needs some if-statements inside the loop computing the $f(x_i)$'s. – Winther Mar 15 '17 at 13:53
  • Adding a delta-distribution to $b(x)$ amounts to imposing a jump in the value of $f(x)$ there. If you wish to "smear out" the jump over one or more intervals of the grid, you can of course do it that way. – hardmath Mar 15 '17 at 13:54

1 Answers1

1

Here's one way to use your original idea of central differences to get a numerical solution by solving a linear system.

Imagine that our domain for $f(x)$ is to be $[0,1]$, and that $f(0)=0$ is specified as the initial condition. For a grid we will choose $n$ equally spaced points, $x_k = k/n$ for $k=0,\ldots,n$.

Now we assemble the following linear system based on the central difference approximation to the derivative:

$$ \frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} = b\left( \frac{x_k + x_{k-1}}{2} \right) $$

for $k = 1,\ldots,n$. Since $x_k - x_{k-1} = 1/n$, the above simplifies to:

$$ f(x_k) - f(x_{k-1}) = (1/n)\cdot b\left( \frac{x_k + x_{k-1}}{2} \right) $$

Now the system of unknowns $f(x_k), k=0,1,\ldots,n$ has a unique solution. Since $f(x_0) = f(0) = 0$ is known, the remaining unknowns are expressed as summations:

$$ f(x_k) = \sum_{i=1}^k (1/n)\cdot b\left( \frac{x_i + x_{i-1}}{2} \right) $$

In other words the solution is just keeping a running sum of the indicated terms. In a calculus course this would be called using the Midpoint Rule to approximate the integral.

Possibly you will object to evaluating the given function $b(x)$ at the midpoints of the intervals $[x_{k-1},x_k]$ on the grounds that these are not grid points. In that case we could replace the values of $b(x)$ at the midpoints with the average of values at the endpoints of intervals:

$$ b\left( \frac{x_k + x_{k-1}}{2} \right) \approx \frac{b(x_k) + b(x_{k-1})}{2} $$

With this substitution in the system of linear equations we obtain the numerical quadrature scheme known as the Trapezoid Rule. Since only right-hand side of the linear system is changed by this replacement, our sketch of an explicit solution is still valid with the necessary substitution inside the running summation.

hardmath
  • 37,015
  • 1
    thanks for your answer and all the valuable comments above. Your approach could also be understood as storing the $f(x_k)$ between the grid points $x_k$, similar to staggered grids in CFD. Thanks again. This was very helpful! – dkoerner Mar 15 '17 at 15:16
  • Since you mention CFD, the pressure variable in an incompressible fluid is a quantity which we only care to determine up to a constant (since adding a uniform constant leaves the pressure gradient unchanged). – hardmath Mar 15 '17 at 15:20