1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Alone Programmer Mar 24 at 19:06
  • $\begingroup$ I am not trying to solve a PDE but to differentiate a function twice. $\endgroup$ – econmajorr Mar 24 at 19:12
  • $\begingroup$ So the accuracy in terms of $\mathcal{O}(h^{n})$ could be find easily by Taylor expansion. $\endgroup$ – Alone Programmer Mar 24 at 19:44
  • $\begingroup$ @Alone Programmer can you say a little more on that? $\endgroup$ – econmajorr Mar 24 at 20:10
  • $\begingroup$ 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. $\endgroup$ – Alone Programmer Mar 24 at 20:12
6
$\begingroup$

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 fifth-order accurate $\mathcal{O}(h^{5})$.

$\endgroup$
  • $\begingroup$ 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$? $\endgroup$ – econmajorr Mar 26 at 8:13
  • $\begingroup$ @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. $\endgroup$ – Alone Programmer Mar 27 at 13:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.