Finite difference method for non-uniform grid

Question

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?

Answer 1

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!