Use of Taylor series expansion to find second derivative for sixth order method

Question

Use Taylor's expansion to derive sixth order method (i.e $\mathcal{O}(h^6)$) for approximating the second derivative ($f '' (x_0)$ ) for given sufficiently smooth function $f(x)$.

I have this things in my mind:

• I can use this seven points: $x_0 - 3h,x_0 - 2h, x_0 - h, x_0, x_0 + h, x_0 + 2h, x_0 + 3h$

$$ f(x_0 + h) = f(x_0) + (h) f'(x_0) + \frac{(h)^2}{2} f''(x_0) + \frac{(h)^3}{3!} f'''(x_0) + \frac{(h)^4}{4!} f''''(x_0) + \frac{(h)^5}{5!} f'''''(x_0) + \frac{(h)^6}{6!}f''''''(x_0) + \frac{(h)^7}{7!}f'''''''(\xi) $$

$$ f(x_0 + 2h) = f(x_0) + (2h) f'(x_0) + \frac{(2h)^2}{2} f''(x_0) + \frac{(2h)^3}{3!} f'''(x_0) + \frac{(2h)^4}{4!} f''''(x_0) + \frac{(2h)^5}{5!} f'''''(x_0) + \frac{(2h)^6}{6!}f''''''(x_0) + \frac{(2h)^7}{7!}f'''''''(\xi) $$

$$ f(x_0 + 3h) = f(x_0) + (3h) f'(x_0) + \frac{(3h)^2}{2} f''(x_0) + \frac{(3h)^3}{3!} f'''(x_0) + \frac{(3h)^4}{4!} f''''(x_0) + \frac{(3h)^5}{5!} f'''''(x_0) + \frac{(3h)^6}{6!}f''''''(x_0) + \frac{(3h)^7}{7!}f'''''''(\xi) $$

$$ f(x_0 - h) = f(x_0) - (h) f'(x_0) + \frac{(h)^2}{2} f''(x_0) - \frac{(h)^3}{3!} f'''(x_0) + \frac{(h)^4}{4!} f''''(x_0) - \frac{(h)^5}{5!} f'''''(x_0) + \frac{(h)^6}{6!}f''''''(x_0) - \frac{(h)^7}{7!}f'''''''(\xi) $$

$$ f(x_0 - 2h) = f(x_0) - (2h) f'(x_0) + \frac{(2h)^2}{2} f''(x_0) - \frac{(2h)^3}{3!} f'''(x_0) + \frac{(2h)^4}{4!} f''''(x_0) - \frac{(2h)^5}{5!} f'''''(x_0) + \frac{(2h)^6}{6!}f''''''(x_0) - \frac{(2h)^7}{7!}f'''''''(\xi) $$

$$ f(x_0 - 3h) = f(x_0) - (3h) f'(x_0) + \frac{(3h)^2}{2} f''(x_0) - \frac{(3h)^3}{3!} f'''(x_0) + \frac{(3h)^4}{4!} f''''(x_0) - \frac{(3h)^5}{5!} f'''''(x_0) + \frac{(3h)^6}{6!}f''''''(x_0) - \frac{(3h)^7}{7!}f'''''''(\xi) $$

• Is there any shortcut in matlab or julia to find this?

Answer 1

One shortcut is that $f''(x_0)$ is symmetric under reflection about $x_0$, so it's reasonable that your estimator should also be symmetric. Thus it should be of the form $$ F(f) = a_0 f(x_0) + a_1 (f(x_0+h) + f(x_0 - h)) + a_2 (f(x_0 + 2h) + f(x_0 - 2h)) + a_3 (f(x_0 + 3h) + f(x_0 - 3h))$$ It should give $0$ for constants, $(x-x_0)^4$ and $(x-x_0)^6$ and $2$ for $(x - x_0)^2$. That makes four linear equations to solve for four variables $a_0, \ldots, a_3$.

Answer 2

Considering any set of points, you could first write the general formula $$f(x+a_ih)=f(x)+a_i h f'(x)+\frac{1}{2} a_i^2 h^2 f''(x)+\frac{1}{6} a_i^3 h^3 f^{(3)}(x)+\frac{1}{24} a_i^4 h^4 f^{(4)}(x)+\frac{1}{120} a_i^5 h^5 f^{(5)}(x)+\frac{1}{720} a_i^6 h^6 f^{(6)}(x)+\frac{1}{5040}a_i^7 h^7 f^{(7)}(x)+O\left(h^8\right)$$ where, in the present case, $a_i=-3+i$ with $(0\leq i\leq 6)$.

Now write $$\sum_{i=0}^6 A_i f(x+a_ih)$$ Set all coefficients equal to $0$ except the coefficient of $f''(x)$ that you make equal to $1$. This would lead to the following set of equations $$A_0+A_1+A_2+A_3+A_4+A_5+A_6=0\tag 1$$ $$-3 A_0-2 A_1-A_2+A_4+2 A_5+3 A_6=0\tag 2$$ $$\frac{1}{2} (9 A_0+4 A_1+A_2+A_4+4 A_5+9 A_6)=1\tag 3$$ $$ \frac{1}{6} ( -27 A_0-8 A_1-A_2+A_4+8 A_5+27 A_6)=0\tag 4$$ $$\frac{1}{24} (81 A_0+16 A_1+A_2+A_4+16 A_5+81 A_6) =0\tag 5$$ $$\frac{1}{120} (-243 A_0-32 A_1-A_2+A_4+32 A_5+243 A_6)=0\tag 6$$ $$\frac{1}{720} (729 A_0+64 A_1+A_2+A_4+64 A_5+729 A_6)=0\tag 7$$ which is easy to solve using matrix calculations.

This leads to $$A_0=\frac 1{90}\quad A_1=-\frac 3{20}\quad A_2=\frac 3{2}\quad A_3=-\frac {49}{18}\quad A_4=\frac 3{2}\quad A_5=-\frac 3{20}\quad A_6=\frac 1{90}$$ This makes $$\sum_{i=0}^6 A_i f(x+a_ih)=h^2 f''(x)+O\left(h^8\right)$$

It is sure that, for the specific values you selectes for the $a_i$, taking into account the symmetry as Robert Israel answered makes things faster.

For sure, you could use the same equations setting to $0$ the rhs of $(3)$ and to $1$ the rhs of $(4)$ to get the third derivative and so on.

For example, the third derivative would be given using $$A_0=\frac 1{8}\quad A_1=-1\quad A_2=\frac {13}{8}\quad A_3=0\quad A_4=-\frac {13}{8}\quad A_5=1\quad A_6=-\frac 1{8}$$ This makes $$\sum_{i=0}^6 A_i f(x+a_ih)=f^{(3)}h^3+O\left(h^7\right)$$

Answer 3

I came up with this Mathematica code

m = 3;
a = Table[ToExpression["a" <> ToString[n]], {n, m, 0, -1}];
vals = Table[
   Normal[Series[f[x], {x, x0, 8}]] /. {x -> x0 + n h}, {n, -m, m}];
vals // TableForm
vals = Join[
   Table[vals[[i]] + vals[[2 m + 2 - i]], {i, m}], {vals[[m + 1]]}];
vals // TableForm
vals = Collect[ExpandAll[Dot[a, vals]], 
   Table[Derivative[n][f][x0], {n, 0, 2 m + 2}], Simplify];
list = Table[
   CoefficientList[vals, Derivative[n][f][x0]][[2]], {n, 0, 2 m + 2, 
    2}];
soln = Flatten[
   Solve[Table[list[[n]] == {0, 1, 0, 0}[[n]], {n, m + 1}], a]];
soln // TableForm
vals /. soln