1D inhomogeneous Poisson PDE with Dirichlet BCs, slow convergence

Question

For an assignment I have to implement a 1D Poisson PDE with inhomogeneous Dirichlet BC's $$\Delta_1 u = f, \quad u(a)=g(a), \: u(b) = g(b) $$ I have managed to make it work, but I am not seeing the expected convergence behaviour. I used the standard (?) discretisation on the unit interval [0,1], $$-u_{j+1} + 2u_j - u_{j-1} = h^2f_j$$

I have the following MATLAB code:

function [sol, xVals, L, fVals] = poissonSolveFD2(f, g, N)

h = 1/(N+2);
h2= h*h;

sol = zeros(N+2,1);
sol(1) = g(0);
sol(N+2) = g(1);

xVals = linspace(0,1,N+2)';

L = makeLaplace(N);

fVals = f(xVals);
fVals = h2*fVals;

fVals(2) = fVals(2) + g(0);
fVals(N+1) = fVals(N+1) + g(1);

sol(2:N+1) = L\fVals(2:N+1);


plot(xVals,sol)

L=makeLaplace(N) returns the tridiagonal matrix $L=\operatorname{tridiag}(-1,2,-1) \in \mathbb{R}^{n \times n}$

This code yields me a solution, but only shows an error reduction of $\mathcal{O}(h)$ instead of $\mathcal{O}(h^2)$.

I have studied the answers at Testing 1D Poisson Solver, but that didn't get me any closer. That question also only considers homogeneous BC's, so I think this question has a right of existence.

It is probably something small and stupid, but for some reason I can't figure it out.

Thanks in advance.

Answer 1

Update: the advice below is necessary, but was written just to expand on the comments above and isn't the actual issue with the code above. See the other answers.

You have a discretized system $Ax=b$, and certain entries of $x$ are specified by your boundary conditions. Let $x_1$ denote the vector of "free" (unknown) values, and let $x_2$ denote the vector of "fixed" (known) values. Similarly partition $A$ and $b$, obtaining the block system

$\left[ \begin{array}{cc} A_{11} & A_{12} \\ A_{21} & A_{22} \end{array} \right]\left[ \begin{array}{c} x_1 \\ x_2 \end{array} \right] = \left[ \begin{array}{c} b_1 \\ b_2 \end{array} \right] $

The values of $x_2$ are specified by your problem, so your task is to compute $x_1 = A_{11}^{-1}(b_1 - A_{12}x_2$).

The values of $b_2$ should be consistent (if they are specified).

Answer 2

The answers so far seem to me to be completely off target. The boundary conditions you have are consistent with the r.h.s. of the equation and the matrix $L$, so they are not the problem.

Instead, there are two things wrong here. First, the usual finite second-order difference approximation to $u''$ is $$ u''(x_0) = h^{-2}\big(u(x_{-1}) - 2u(x_0) + u(x_1)\big), $$ which has the opposite sign to what you have. In effect, if you use $\mathrm{diag}(-1,2,-1)$ you end up solving $\Delta u=-f$, which is a bug your testing doesn't seem to have caught.

Second, when you have $N$ interior points and $N+2$ points in total, with $x_1 = 0$ and $x_{N+2} = 1$, the $j$-th point is $$ x_j = \frac{j-1}{N+1}, $$ and so $$ h = \frac{1}{N+1}, $$ whereas you have $h=1/(N+2)$. This is the reason for incorrect order of convergence. The easiest way here to test for this kind of mistake is to use the exact solution $u=4x(x-1)$ to the equation $\Delta u=8$ with $g(x)=0$, so that $u(\frac12)=-1$. Because the error of the FD formula behaves as $u^{(4)}$, the scheme is exact in this case, so gettting any error other than round-off error is a bug.

Third, for how to set boundary conditions see the answer by Patrick Sanan, but that doesn't seem to be a problem here as far as I can see. The matrix $A_{12}$ has first row $(-1,0)$ and last row $(0,-1)$, with zeros everywhere else, so you are correct to add the boundary conditions to the rhs at $x_2$ and $x_{N+1}$ (although the sign would still be wrong as above). The units of $h^2f$ and $g$ are also correct ($f$ has units $u/x^2$ (like $u''$), $g$ has units $u$, $h$ has units $x$ and $h^2f$ has units $u$, same as $g$).

Fourth, for testing, it is much easier to pick a known exact solution, say $v(x)=e^x$, go backwards to work out $f(x)=e^x$ and $g(x)=e^x$, and then apply the scheme and compare $u$ with $v$. This is commonly known as the method of manufactured solutions.

Here is a modified version of your code

function [sol, xVals, L, fVals] = poissonSolveFD2(f, g, N)

h = 1/(N+1);
h2= h*h;

sol = zeros(N+2,1);
sol(1) = g(0);
sol(N+2) = g(1);

xVals = linspace(0,1,N+2)';

L = spdiags(repmat([1,-2,1], N, 1), -1:1, N, N);

fVals = f(xVals);
fVals = h2*fVals;

fVals(2) = fVals(2) - g(0);
fVals(N+1) = fVals(N+1) - g(1);

sol(2:N+1) = L\fVals(2:N+1);

plot(xVals,sol)

Simple check that it works:

>> u = @(x) exp(x);
>> e1 = norm(poissonSolveFD2(u, u, 100) - u(linspace(0,1,102)'), 'inf')
e1 =
   1.7307e-06
>> e2 = norm(poissonSolveFD2(u, u, 201) - u(linspace(0,1,203)'), 'inf')
e2 =
   4.3269e-07
>> e1 / e2
ans =
    3.9999