[Math] Finite difference method for non-uniform grid

finite differencesnumerical methods

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: enter image description here

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?

Best Answer

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!

Related Question