4
$\begingroup$

I am solving the following PDE;

$$ \nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = \rho, $$

where $\rho(0.5,0.5) = 2$ (zero elsewhere), $0\leq x,y\leq1$ and the boundaries will satisfy the condition $u = 0$.

I am solving by finite differences using both a 5 point second order stencil and a 9 point fourth order stencil. It is my understanding that the solution to this PDE is the zero function and this is what I am using to calculate rates of convergence. I have used various error measures including the average difference, max difference and the value at the point $(0.5,0.5)$.

For each error measure I am getting convergence rates of approximately order 2 for both stencils, is there a reason why the fourth order stencil would not converge with order 4?

This is a repeat of a question I asked on MSE - https://math.stackexchange.com/questions/1073945/the-rate-of-convergence-for-finite-difference-methods-for-poissons-equation-wit/1074255#1074255

Although the answers there have helped a lot I would appreciate anyone who could elaborate a little further?

$\endgroup$
  • $\begingroup$ The answer you got on MSE looks good to me, why post a duplicate question? That said, have you absolutely made sure your implementation is correct? When you want to test a method's implementation, you should start by picking some known functions $u(x,y)$, finding the corresponding b.c. and r.h.s. $\rho$, and comparing with numerical approximations. $\endgroup$ – Kirill Dec 19 '14 at 18:38
  • $\begingroup$ No I am not completely certain...mainly because it does not converge how I expected! I am trying to test it on a known solution but I am having issues deciding how to deal with boundary points since the nine point stencil extends out two nodes...Duplicate question is because it was suggested to do so... $\endgroup$ – Josh Greenhalgh Dec 19 '14 at 18:38
  • $\begingroup$ See here about cross-posting; it's a slight breach of etiquette. For boundary points, if you know an exact solution, set the boundary points and the points just beyond the boundary to their exact known values; this gives you the linear equations you need to solve. An alternative is to use a smaller stencil at points next to the boundary, so that no extra ghost points are needed. $\endgroup$ – Kirill Dec 19 '14 at 18:40
8
$\begingroup$

Let's examine the one-dimensional three-point stencil case in detail, because I think it's important to be clear just how this behaviour arises, and what it means to set a point to a certain value in a finite-difference grid when the underlying function is discontinuous. The equation will be $$ u''(x) = \rho(x). $$ Instead of using the interval $[-1,1]$ with discontinuity at $\frac12$, I will use the interval $[-\frac12,\frac12]$ with the discontinuity placed at $0$. The grid size will be $h$, and I will only have to consider the interval $[-\frac h2,\frac h2]$ around the grid point $0$.

First, in a finite-difference approximation, we approximate $$ u''(0) = \frac{\hat u(-h) - 2\hat u(0) + \hat u(h)}{h^2} $$ and solve $$ \frac{\hat u(-h) - 2\hat u(0) + \hat u(h)}{h^2} = \hat \rho(0). $$ (Here the variables with hats are the numerical approximations on the grid to the variables without hats.)

But this is a very bad specification of the problem, because your function $\rho(0) = 2$ at $x=0$, and $\rho(x) = 0$ everywhere else, is discontinuous. In particular, if we shift the grid point $0$ to either side by some tiny amount $\epsilon$, the numerical solution changes entirely and becomes exactly zero. This means that this is a terribly misspecified problem.

We can make sense of it by converting it to an equivalent finite-volume formulation, where it will make much more sense.

In a finite-volume method, we solve $$ \int_{-h/2}^{h/2} u''(x)\,dx = \int_{-h/2}^{h/2}\rho(x)\,dx, $$ by picking $u(x)$ to be a suitable approximation to the unknown function. Let's pick, on the interval $[-\frac h2,\frac h2]$ the approximation $$ u(x) = \hat u(-h) \phi(x+h) + \hat u(0) \phi(x) + \hat u(h) \phi(x-h), $$ where $\phi(x)$ is the basis function $$ \phi(x) = \max\left(1-\frac{|x|}{h}, 0\right), $$ (it's a piecewise linear function that goes from $0$ at $-h$ to $1$ at $0$ to $0$ at $+h$, thus interpolating between grid points.). The approximation $u(x)$ is a weighted sum of three basis functions that look like this:

Basis functions

We can then compute $$ \phi''(x) = \frac1h\delta(h-|x|) - \frac2h \delta(x), $$ so that the approximation to the integral is $$ \int_{-h/2}^{h/2} u''(x) = \frac{\hat u(-h) -2\hat u(0)+\hat u(h)}{h}, $$ and the finite-volume approximation to our equation becomes $$ \frac{\hat u(-h) -2\hat u(0) + \hat u(h)}{h} = \int_{-h/2}^{h/2} \rho(x)\,dx = h \hat \rho(0). $$

There are two important things here. First, this is equivalent to the finite-difference formulation in that we end up solving the same equations. Second, the discontinuity in $\rho$ is given a very precise meaning: when we use the value $2$ for $\hat \rho(0)$, we are saying that this is the average value of $\rho$ on the interval $[-h/2,h/2]$: $$2 = \hat\rho(0) = \frac1h \int_{-h/2}^{h/2} \rho(x)\,dx.$$ This interpretation is not available in the finite-difference formulation. It is also not so sensitive to the location of the discontinuity: if the discontinuity were at some small distance $\epsilon$ away from $0$, the average value would be almost the same, but the value at $0$ might be completely different.

But if we say that $2$ is the average value of $\rho$ near the grid point, we can then go back and compute the exact solution of the equation with the right-hand side given by $$ \tilde \rho(x) = 2[-h/2 < x < h/2]. $$ (We pick a function of our own choice that gives the right average.) In this case, the exact solution will be $$ \tilde u(x) = 2 \int_{-h/2}^{h/2} G(x; u)\,du \approx 2h G(x; 0), $$ in terms of the Green function for the Poisson equation. In the two-dimensional case it will be $$ \approx 2h^2 G(x, y; 0, 0), $$ as in the other (correct) answer on MSE.

Finally, the outcome of all this is that when you say that you compare your numerical approximate solution with the exact solution $u(x)=0$, this is wrong. The exact solution should not be zero, it should be $$ \approx 2h^2 G(x, y; 0,0). $$ It therefore should make perfect sense that the fourth-order solution does not converge to zero with order $4$: it should converge to the correct solution, which is not zero but has magnitude of order $O(h^2)$.

If you do want to get zero as the numerical solution and compare with the mathematical solution, you should set $\hat\rho(0)$ to be the average value of $\rho$, which is $0$, not $2$.

The fact that the exact solution depends on the chosen grid size indicates that this is not a good way to check whether you implemented the method correctly. A very straightforward technique known as the method of manufactured solutions is better for this.

$\endgroup$
  • $\begingroup$ I really appreciate the work you have put into this answer, it has made things much clearer for me. The method I am using does converge as expected (order 4) for a manufactured solution with no discontinuities - so I am happy that my implementation is correct. In terms of this particular source term - the reason it seems to be converging to the zero function is because of the scaling of the greens function by $h^2$ which in the limit is zero? I still have a bit of work to do to fully understand the convergence I am seeing but I am much closer to being able to do this - thanks again!! $\endgroup$ – Josh Greenhalgh Dec 19 '14 at 21:13
  • $\begingroup$ @JoshGreenhalgh Upvote answers that you find helpful. Yes, if the exact solution is $O(h^2)$, it will converge to zero as $h\to0$, making it look like you see your numerical method converges to the "right" solution at order 2. $\endgroup$ – Kirill Dec 19 '14 at 21:17
  • $\begingroup$ I am not yet worthy to upvote...but it is done in spirit! $\endgroup$ – Josh Greenhalgh Dec 19 '14 at 21:42
  • $\begingroup$ @Kirill Beautiful ! +1 $\endgroup$ – Sektor Dec 20 '14 at 10:52
0
$\begingroup$

Think about what "second order" and "fourth order" approximations actually mean. A first order method can exactly reproduce a first order (linear) function, a second order can reproduce a quadratic, etc.

So say if the exact solution to a given problem is linear - if so, a 100th order method won't give you any benefit over a 1st order method, because both can exactly represent the solution. I suspect this might be a similar case with your zero function solution. Try something with a sinusoidal forcing term or boundary condition.

$\endgroup$
  • $\begingroup$ I don't see how this answers the question, which is about a very specific right-hand side; the answer shouldn't be "try a different r.h.s.". Can you explain the observed results directly? $\endgroup$ – Kirill Dec 19 '14 at 18:44
  • $\begingroup$ The mathexchange answer is effectively the same thing, just with some rigor and without my hand-waving. He showed the actual answer to the problem is itself $O(h^2)$. If the analytical answer to a problem is $O(h^2)$ then an $O(h^4)$ approximation will never be any better than an $O(h^2)$ approximation. The answer really is "try a different RHS" or "try different BCs", because the example given won't show a difference between 2nd and higher order schemes. Truncation errors will be identically 0. $\endgroup$ – Aurelius Dec 19 '14 at 18:55
  • $\begingroup$ The question is "is there a reason why the fourth order stencil would not converge with order 4?", which I interpret as asking about this precise r.h.s. I agree that if the question were "how should I test these methods", the answer would be "try a better r.h.s.". $\endgroup$ – Kirill Dec 19 '14 at 19:00
  • $\begingroup$ "is there a reason why the fourth order stencil would not converge with order 4?" - yes, and the reason is that the actual solution is order 2, so higher convergence is not possible. One can write out a mathematical proof, but it should be intuitively obvious. If you're trying to approximate a function which is a straight line, a tenth order approximation cannot be more accurate than a first order one, since the first order would be exact. $\endgroup$ – Aurelius Dec 19 '14 at 19:11
  • $\begingroup$ The fact that if the solution is order 2 a method of order 4 wound not show any improvement is clear to me...however why the 'solution' is order 2 is not! But thanks anyway...my method does converge with order 4 for a problem involving a trigonometric source term just not this one $\endgroup$ – Josh Greenhalgh Dec 19 '14 at 19:29

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.