# Testing 1D Poisson Solver

I'm trying to test a simple 1D Poisson solver to show that a finite difference method converges with $\mathcal{O}(h^2)$ and that using a deferred correction for the input function yields a convergence with $\mathcal{O}(h^4)$.

So, the equation is $- u'' = f$ with boundary conditions $u(0) = u(1) = 0$. The method I'm trying to use is using the discretized operator $$A = \left[ \begin{array}{c} 2&-1&0&0&0&0 \\-1&2&-1&0&0&0 \\ 0&-1&2&-1&0&0 \\ 0&0&-1&2&-1&0 \\ 0&0&0&-1&2&-1 \\ 0&0&0&0&-1&2\end{array}\right]$$

(the example matrix is for $h = \frac{1}{5}$.) Then solve for $Au=h^2f$. I've shown that theoretically this should converge with $\mathcal{O}(h^2)$, but when I test it on Matlab, I'm getting only $\mathcal{O}(h)$ convergence.

Then, I'm trying what my course instructor called "deferred correction", and altering $f$ before solving. I concluded that the correction should be $f \mapsto f + \frac{h^2}{12} Af$. I've shown that this should converge with $\mathcal{O}(h^4)$, but in Matlab I still get $\mathcal{O}(h)$.

Here's the Matlab script:

function [u err] = threeptsolve(ureal, du2, h)

% INPUT: 'ureal' is the function handle for the real solution.
%        'du2' is the function handle for the second derivative of 'ureal'
%        'h' is the step size
% OUTPUT: 'u' is the approximated solution
%         'err' is the error at each point

x = [0:h:1]';
n = length(x);
f = -du2(x);
realu = ureal(x);

A = 2 * eye(n);
A = A + diag(-1*ones(n-1,1), 1) + diag(-1*ones(n-1,1), -1);
A = (1/h^2) * A;

% uncomment if using deferred correction
% f = f + h^2/12 * A * f;

u = A\f;

err = (realu - u);

end


When I try this with some sample smooth functions (with appropriate boundary values), and then try again with $h/2$, I get a vector of (approximate) twos when I compute err1 ./ err2(1:2:end).

Is my math wrong, or is it my code?

• What is your ureal? If it isn't sufficiently smooth, you may lose high order accuracy. I doubt for a test case that this is the issue though. – Godric Seer Mar 12 '13 at 15:50
• @GodricSeer, I use polynomials that are zero at the boundaries. E.g. ureal = @(x)(-10*x.^4 + 5*x.^3 + 2*x.^2 + 3*x). – jake Mar 12 '13 at 15:58
• Try with this: you have make a minor mistake from the beginning in writing your equation to solve as $Au=h^2 f$ but in fact it is $Au=h f$ may be that is the problem. – atiliomor Nov 6 '16 at 16:25

Your formulation of $A$ assumes that $u_0$ and $u_{n+1}$ are zero, which is correct. However, you are then including your boundaries at $u_1$ and $u_n$. Your exact answer satisfies these later BC's but not the imposed BC's in $A$. These discontinuities floating around cause you to lose your higher order accuracy. To get second order convergence you only need to change one line:

x = [0:h:1]';


to

x = [h:h:1-h]';


And I get a column of 4's (for 2nd order conv.) when I execute err1 ./ err2(2:2:end-1) (Note the shift by 1 index since the solutions now line up at the even indices rather than the odd). I have not gotten 4th order yet from your "deferred correction", however this solve a part of your problem.

• Thank you! I've been looking at this all day and didn't catch that. Now, I'm finally on to fixing the deferred correction method. (It's more possible for that one that my math is off. I'll try re-calculating.) – jake Mar 12 '13 at 18:30

For defect correction you will need a 4th-order discretization. Let $A_2 u = f_2$ be your second-order discretization and let $A_4 u = f_4$ be a 4th-order discretization. Defect corrections then become $$A_2 u^{(0)} = f_2 \\ A_2 u^{(k)} = f_4 - A_4 u^{(k-1)} + A_2 u^{(k-1)} \quad k=1,2,...$$ As the defect corrections converge $u^{(k)} \approx u^{(k-1)}$ so the $A_2 u$ terms cancel and you have solved the 4th-order discretization without ever "inverting" $A_4$. There is ample theory on how many iterations are needed, see for example Hackbush, Multi-grid methods and applications, Springer 1985.