# What is the maximum attainable accuracy with a given set of $\alpha,\beta$?

I am using LeVeque's book: https://faculty.washington.edu/rjl/fdmbook/

Suppose I want to compute $$u''$$ using FDM with $$\alpha=\beta=2$$ (centered) so the FDM is

$$u''=\sum_{m = - \alpha}^\beta a_mu(x_i+mh)$$

Before acutally determining the coeficients ($$a_m$$'s) can we say anthying about how what the maximum attainable accuracy is?

Leveques does present methods the order of accuracy when the coefficients are already determined, but exactly how to determine the maximum order beforehand.

• The general attainable accuracy depends on lots of things and FDM stencil is just one of them. You need to express what is your application and possibly the PDE that you are trying to solve here. Mar 24 '20 at 19:06
• I am not trying to solve a PDE but to differentiate a function twice. Mar 24 '20 at 19:12
• So the accuracy in terms of $\mathcal{O}(h^{n})$ could be find easily by Taylor expansion. Mar 24 '20 at 19:44
• @Alone Programmer can you say a little more on that? Mar 24 '20 at 20:10
• Try to expand $u(x_i+mh)$ terms around $x_i$ due to the fact that $h$ is small. Then try to find $u^{“}$ and eliminate unwanted terms. Mar 24 '20 at 20:12

If you are looking for general accuracy of your stencil, it could be extracted by using Taylor expansion. Basically, if I want to write down your stencil explicitly, it contains terms of $$u(x_{i}-2h)$$, $$u(x_{i}-h)$$, $$u(x_{i})$$, $$u(x_{i}+h)$$, and $$u(x_{i}+2h)$$. Let's look at Taylor expansion of each of these terms:

$$u(x_{i}-2h) - u(x_{i}) = - u^{'}(x_{i})(2h) + \frac{1}{2} u^{''}(x_{i}) (4 h^{2}) - \frac{1}{6} u^{'''}(x_{i}) (8h^{3}) + \frac{1}{24} u^{''''}(x_{i}) (16h^{4}) + \mathcal{O}(h^{5})$$

$$u(x_{i}-h) - u(x_{i}) = - u^{'}(x_{i}) h + \frac{1}{2} u^{''}(x_{i}) h^{2} - \frac{1}{6} u^{'''}(x_{i})h^{3} + \frac{1}{24} u^{''''}(x_{i})h^{4} + \mathcal{O}(h^{5})$$

$$u(x_{i}+h) - u(x_{i}) = u^{'}(x_{i})h + \frac{1}{2} u^{''}(x_{i}) h^{2} + \frac{1}{6} u^{'''}(x_{i}) h^{3} + \frac{1}{24} u^{''''}(x_{i}) h^{4} + \mathcal{O}(h^{5})$$

$$u(x_{i}+2h) - u(x_{i}) = u^{'}(x_{i})(2h) + \frac{1}{2} u^{''}(x_{i})(4 h^{2}) + \frac{1}{6} u^{'''}(x_{i}) (8h^{3}) + \frac{1}{24} u^{''''}(x_{i}) (16h^{4}) + \mathcal{O}(h^{5})$$

You can write these equations in matrix form:

$$\begin{bmatrix} -2h & 2h^{2} & -\frac{4}{3} h^{3} & \frac{2}{3} h^{4} \\ -h & \frac{1}{2} h^{2} & -\frac{1}{6} h^{3} & \frac{1}{24} h^{4} \\ h & \frac{1}{2} h^{2} & \frac{1}{6} h^{3} & \frac{1}{24} h^{4} \\ 2h & 2h^{2} & \frac{4}{3} h^{3} & \frac{2}{3} h^{4} \end{bmatrix} \begin{bmatrix} u^{'}(x_{i}) \\ u^{''}(x_{i}) \\ u^{'''}(x_{i}) \\ u^{''''}(x_{i}) \end{bmatrix} = \begin{bmatrix} u(x_{i} - 2h) - u(x_{i}) \\ u(x_{i} - h) - u(x_{i}) \\ u(x_{i} + h) - u(x_{i}) \\ u(x_{i} + 2h) - u(x_{i}) \end{bmatrix}$$

When you solve this linear equation where you know $$h$$ and the values of $$u(x_{i}-2h)$$, $$u(x_{i}-h)$$, $$u(x_{i})$$, $$u(x_{i}+h)$$, and $$u(x_{i}+2h)$$, you would get unknowns of $$u^{'}(x_{i})$$, $$u^{''}(x_{i})$$, $$u^{'''}(x_{i})$$, and $$u^{''''}(x_{i})$$. You are interested particularly in $$u^{''}(x_{i})$$. The coefficients $$a_{m}$$ also will be determined here automatically.

In fact the inverse of above matrix could be extracted by sympy as:

$$\left[\begin{matrix}\frac{0.111111111111111}{0.444444444444444 h + 0.888888888888889} & - \frac{0.444444444444444 h + 0.444444444444444}{h \left(0.444444444444444 h + 0.888888888888889\right)} & \frac{0.444444444444444 \left(h + 1\right)}{h \left(0.444444444444444 h + 0.888888888888889\right)} & \frac{0.333333333333333}{- 1.33333333333333 h - 2.66666666666667}\\\frac{0.333333333333333}{h \left(- 1.33333333333333 h - 2.66666666666667\right)} & \frac{1.0 \left(0.740740740740741 h + 1.03703703703704\right)}{h^{2} \left(0.444444444444444 h + 0.888888888888889\right)} & \frac{1.0 \left(0.444444444444444 h + 1.33333333333333\right)}{h^{2} \left(0.444444444444444 h + 0.888888888888889\right)} & \frac{1.0 \left(0.111111111111111 h - 0.444444444444444\right)}{h^{2} \left(1.33333333333333 h + 2.66666666666667\right)}\\- \frac{0.666666666666667}{h^{2} \left(0.444444444444444 h + 0.888888888888889\right)} & \frac{1.33333333333333}{h^{2} \left(0.444444444444444 h + 0.888888888888889\right)} & - \frac{1.33333333333333}{h^{2} \left(0.444444444444444 h + 0.888888888888889\right)} & \frac{2.0}{h^{2} \left(1.33333333333333 h + 2.66666666666667\right)}\\\frac{4.0}{h^{3} \left(1.33333333333333 h + 2.66666666666667\right)} & \frac{10.6666666666667 h + 5.33333333333334}{h^{4} \left(- 1.33333333333333 h - 2.66666666666667\right)} & \frac{16.0}{h^{4} \left(- 1.33333333333333 h - 2.66666666666667\right)} & \frac{1.33333333333333 h - 5.33333333333333}{h^{4} \left(- 1.33333333333333 h - 2.66666666666667\right)}\end{matrix}\right]$$

So:

$$a_{-2} = \frac{0.333333333333333}{h \left(- 1.33333333333333 h - 2.66666666666667\right)}$$

$$a_{-1} = \frac{\left(0.740740740740741 h + 1.03703703703704\right)}{h^{2} \left(0.444444444444444 h + 0.888888888888889\right)}$$

$$a_{1} = \frac{\left(0.444444444444444 h + 1.33333333333333\right)}{h^{2} \left(0.444444444444444 h + 0.888888888888889\right)}$$

$$a_{2} = \frac{\left(0.111111111111111 h - 0.444444444444444\right)}{h^{2} \left(1.33333333333333 h + 2.66666666666667\right)}$$

$$a_{0} = - \frac{2.5}{h^{2}}$$

This is the code to generate these coefficients:

import sympy as sp
sp.init_printing(use_latex='mathjax')

h = sp.symbols('h')
A = sp.Matrix( [[-2*h,2*h**2,-(4/3)*h**2,(2/3)*h**4],[-h,(1/2)*h**2,-(1/6)*h**3,(1/24)*h**4],[h,(1/2)*h**2,(1/6)*h**3,(1/24)*h**4],[2*h,2*h**2,(4/3)*h**3,(2/3)*h**4]]) # Creates a matrix.
A_inverse = A.inv()

print(sp.printing.latex(sp.simplify(A_inverse)))

print(sp.printing.latex(sp.simplify(A_inverse[1,0])))
print(sp.printing.latex(sp.simplify(A_inverse[1,1])))
print(sp.printing.latex(sp.simplify(A_inverse[1,2])))
print(sp.printing.latex(sp.simplify(A_inverse[1,3])))

print(sp.printing.latex(sp.simplify(-A_inverse[1,0]-A_inverse[1,1]-A_inverse[1,2]-A_inverse[1,3])))


So the explicit answer of your question is that: this stencil is fourth-order accurate and the error is in the order of $$\mathcal{O}(h^{5})$$.

• thank you for this very thorough answer! Simple question: Can I use scipy package to actualy give me the optimal coefficients for any set of $\alpha, \beta$? Mar 26 '20 at 8:13
• @econmajorr what do you mean using scipy? for any given $\alpha$ and $\beta$ the procedure is the same as here but you need to solve $\alpha + \beta + 1$ linear equations. So, if you are asking to solve this linear equations numerically by using scipy.linalg, the answer is: It might be possible, but you need to try it to know for sure. Mar 27 '20 at 13:44
• @AloneProgrammer I understand the exact arguement for saying the order is 5? I know that our Taylors have error term of $O(h^5)$ but what is the reason that we can conclude that final stencils are fifth order accurate? Mar 23 at 21:58
• @AloneProgrammer I have implemented the method. It seems like it is order 4. Are certain that this is order 5? Mar 24 at 13:31
• @k.dkhk Well, a rough explanation is that the accuracy of your solution is in the order of $h^{4}$ and the error that you are getting due to truncating terms of higher orders is in the order of $h^{5}$. In another word, $h^{4}$ is a really rough estimation of your accuracy and $h^{5}$ is a rough estimation of your error getting out of this stencil. Mar 24 at 14:26