I found this method in this book for solving the Poisson equation with error converging with $\mathcal{O}(h^4)$. However, when I try to implement it for the 1D equation $-u''(x)=f(x)$, it only converges with around $\mathcal{O}(h^2)$. The math is as follows:

Using Taylor Series, we get $$\frac{1}{h^2}[-u_{i-1} + 2u_i - u_{i+1}] = -u''_i-\frac{h^2}{12}u^{(4)}_i + \mathcal{O}(h^4).$$ Substituting in $$-u^{(4)} = \frac{1}{h^2}[-u''_{i-1} + 2u''_i - u''_{i+1}] + \mathcal{O}(h^2)$$ we get $$\frac{1}{h^2}[-u_{i-1} + 2u_i - u_{i+1}] = f_i+\frac{1}{12}(f_{i-1} -2f_i +f_{1+1})_i + \mathcal{O}(h^4)$$ where $u_i =u(x_i)$ and $f_i = f(x_i)$.

This means that solving the system $$\frac{1}{h^2}Au = f-\frac{1}{12}Af$$ where $A$ is of the form $$ A = \left[ \begin{array}{c} 2&-1&0&0 \\-1&2&-1&0 \\ 0&-1&2&-1 \\ 0&0&-1&2\end{array}\right] $$ should result in an error that converges with $\mathcal{O}(h^4)$. But when I implement this (in Matlab), halving the step-size $h$ results in approximately quartering the error.

Is there something wrong with my math?

(If necessary, I can include the Matlab code here. But as of now I'm thinking the problem lies with the math.)


I haven't checked the exact coefficients in your formulas, but there are typically two difficulties with this scheme:

  • Is your solution smooth enough to do the Taylor expansion to this order? For example, if your function $f$ is only once continuously differentiable, you can't expect the solution to be four times continuously differentiable.

  • What do you do with boundary conditions? That is, what do you use for $f_{i-1}$ at the left boundary, and similarly $f_{i+1}$ at the right?

  • $\begingroup$ So, in the last few minutes that I've been playing with this, I've switched the function to $u(x) = \sin(2 \pi x)$, and now it works as intended (4th order error convergence), but before I was also using a $C^{\infty}$ function, $u(x) = xe^{x-1}-x$, and I can't figure out why that won't work. $\endgroup$
    – jake
    Apr 15 '13 at 3:01
  • $\begingroup$ Oh, and the BCs are $u(0)=u(1)=0$. I should have mentioned that. $\endgroup$
    – jake
    Apr 15 '13 at 3:03
  • $\begingroup$ @jake: The only issue I see is that your original $u$ seems incompatible with the right hand side of your equation system. $A$ is clearly designed with $u(0)=u(1)=0$ in mind. But you reuse this matrix for your right hand side. This is only OK if $f(0)=f(1)=0$ is true. But that doesn't seem right for $xe^{x-1}-x$ $\endgroup$
    – sellibitze
    Apr 20 '19 at 11:12

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.