Finite difference method for non-uniform grid

It's been two days that I've been stuck on this problem:

Given a regular function $u$, and a non-uniform grid, where every node has a non-constant distance from another, I want to find $u'(x_i)$ and get some information about the error.

Here's a picture of the situation:

The standard way is to expand with Taylor in $x_i$. By doing this, I got this system, given from these four equations:

(i) $u(x_{i+1})=u(x_i) + u'(x_i)h_i + u''(x_i) (h_i)^2/2 + u'''(x_i)h_i^3/6 + u^{(4)}(z)h_i^4/24$

(ii) $u(x_i)=u(x_i)$

(iii) $u(x_{i-1})=u(x_i) - u'(x_i)h_{i-1} + u''(x_i) h_{i-1}^2/2 - u'''(x_i)h_{i-1}^3/6 + u^{(4)}(z)h_{i-1}^4/24$

(iv) $u(x_{i+2})=u(x_i) + u'(x_i)(h_i+h_{i+1}) + u''(x_i) (h_i+h_{i+1})^2/2 + u'''(x_i)(h_i+h_{i+1})^3/6 + u^{(4)}(z)(h_i+h_{i+1})^4/24)$

Now, as done here, i do a linear combination of the equations, in order to get $u'(x_i)$. The unknown coefficients are $\alpha,\beta,\gamma,\delta$, while $h_{i-1},h_i,h_{i+1}$ are given.

I got this (linear) system:

(i) $\alpha+ \beta + \gamma +\delta=0$

(ii) $\alpha h_i + \delta(h_i + h_{i+1}) - \gamma h_{i-1}=1$ ($u'(x_i)$ coefficients)

(iii) $\alpha h_i^2/2 + \beta h_{i-1}^2 /2 + \delta (h_i + h_{i+1})^2/2 =0$

(iv) $\alpha h_i^3/6 - \gamma h_{i-1}^3/6+ \delta(h_i + h_{i+1})^3/6 =0$

It has a unique solution, but it's really ugly, and I computed it with Maxima or other numerical tools. I don't know how to proceed...any hint?

• The solution you get is your formula for the coefficients of the linear combination. In general it can be ugly, especially if the $h_i$ are just symbolic.
– Ian
Oct 13, 2017 at 14:41
• So once I found this coefficients, I have to put them in the equation $\alpha u(x_{i+1} + \beta u(x_i) + \gamma u(x_{i-1}) + \delta (u(x_{i+2})= u'(x_i) + O(max(h_{i-1},h_i, h_{i+1}))$...right?
– VoB
Oct 13, 2017 at 14:50
• If you did it right then that equation will hold, yes. (Actually you should have some power of the maximum instead; for four points you should get much better than first order.)
– Ian
Oct 13, 2017 at 14:50
• The argument of the BIG-O notation is right? Thanks so much
– VoB
Oct 13, 2017 at 14:51
• Depends completely on how the signal is sampled. Oct 13, 2017 at 20:02

As has been observed, the answer the process is basic linear algebra but the answers are likely to be complicated. If you're just interested in the answer, here's a little Mathematica code that expresses $u'(x)$ in terms of three nearby values $u(x+h_1)$, $u(x+h_2)$, and $u(x+h_3)$:

eqs = {
u[x + h] == Normal[Series[u[x + h1], {h1, 0, 3}]],
u[x + h2] == Normal[Series[u[x + h2], {h2, 0, 3}]],
u[x + h3] == Normal[Series[u[x + h3], {h3, 0, 3}]]
};
dq = Simplify[u'[x] /. First@Solve[eqs, {u'[x], u''[x], u'''[x]}]];


If we pass the result to Mathematica's TeXForm and replace each hi with h_i, we get the following:

$$\frac{h_1^3 \left(h_2^2 \left(u\left(h_3+x\right)-u(x)\right)+h_3^2 \left(u(x)-u\left(h_2+x\right)\right)\right)+h_1^2 \left(h_2^3 \left(u(x)-u\left(h_3+x\right)\right)+h_3^3 \left(u\left(h_2+x\right)-u(x)\right)\right)-h_2^2 \left(h_2-h_3\right) h_3^2 \left(u(x)-u\left(h_1+x\right)\right)}{h_1 \left(h_1-h_2\right) h_2 \left(h_1-h_3\right) \left(h_2-h_3\right) h_3}$$

I guess the result should be $O(h)$, where $h$ is the maximum of the $h_i$s. If we replace $u(x)$ with the arbitrary cubic $ax^3+bx^2+cx+d$ and expand, we find that the result is exact.

Hope that helps!

• So, let me know if I make a mistake: you solved (with Mathematica) the first system I wrote and the following expression is $u'(x_i)$, right? Then, I can't understand how u say that the result is $O(max(h_{i}s))$
– VoB
Oct 14, 2017 at 10:00
• I guess my system is actually slightly different from yours, though it should be easy to translate back and forth. They are related via \begin{align} h_1 &= -h_{i-1} \\ h_2 &= h_i \\ h_3 &= h_{i+1} \end{align} I honestly didn't put a lot of thought into the error term and, in particular, I did not explicitly solve for any $O$ term. I'm thinking, though, that any $O$ terms should be at least $O(h^4)$ and, of course, they can be combined. In situations with special symmetry, you can sometimes do better but I don't see that happening here. Oct 14, 2017 at 13:59
• Ok, it's clear ! ;)
– VoB
Oct 15, 2017 at 9:02