Testing 1D Poisson Solver

Question

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?

Answer 1

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.

Answer 2

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 \begin{equation} 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,... \end{equation} 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.