3
$\begingroup$

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.

$\endgroup$
  • 4
    $\begingroup$ You haven't applied the boundary conditions correctly. You can't just add them in to the RHS on those nodes. You actually have to correct the equations to yield what you want. $\endgroup$ – Bill Barth Jun 24 '15 at 2:43
  • $\begingroup$ So what would you suggest? Because I really can't see the problem... $\endgroup$ – Nigel Overmars Jun 24 '15 at 7:22
  • 2
    $\begingroup$ As Bill Barth mentioned, your BCs are wrong. Notice that they don't share the same units, so you can't add them together. To compute that sum you need to take into account the influence of the matrix L in the fixed components of u. One way of doing this is considering that the first column will influence $u(a)$ and the last column $u(b)$. $\endgroup$ – nicoguaro Jun 24 '15 at 13:53
  • 1
    $\begingroup$ @BillBarth I think he applied the boundary conditions correctly. $\endgroup$ – Kirill Jun 25 '15 at 0:22
2
$\begingroup$

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).

$\endgroup$
  • $\begingroup$ A (much longer) discussion of the same procedure is given in lectures 21.6 and 21.65 here: math.tamu.edu/~bangerth/videos.html $\endgroup$ – Wolfgang Bangerth Jun 24 '15 at 21:10
  • 2
    $\begingroup$ Doesn't this lead to just the code in question? The way the question was written, $A_{12}$ has just two $-1$ elements, so $b_1-A_{12}x_2$ just adds the boundary conditions to the rhs just like in the code. Or am I wrong? $\endgroup$ – Kirill Jun 25 '15 at 0:18
  • $\begingroup$ I admit that I didn't read the code carefully - I was just attempting to clarify what was in the comments. $\endgroup$ – Patrick Sanan Jun 25 '15 at 12:36
5
$\begingroup$

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
$\endgroup$
  • $\begingroup$ It's not the first column $(2,-1,0)$ and the last one $(0,-1,2)$? $\endgroup$ – nicoguaro Jun 25 '15 at 15:31
  • $\begingroup$ @nicoguaro No, the "2" matches with the value $u_2$ at interior point $x_2$, so it goes into the $A_{11}$ matrix. The contents of vector $x_2$ are the values $u_1$ and $u_{N+2}$, which get multiplied with $(1,0)$ and $(0,1)$ at points $x_2$ and $x_{N+1}$. Apart from those two rows, $A_{12}$ is all zeros. This is maybe easier to see from the original equations by just moving the knowns $u_1$ and $u_{N+2}$ from the l.h.s. to the r.h.s. $\endgroup$ – Kirill Jun 25 '15 at 21:21
  • $\begingroup$ Is not the diagonal term 2? like this \begin{bmatrix} -2 & 1 & \cdots & 0 & 0 \\ 1 & -2 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & -2 & 1 \\ 0 & 0 & \cdots & 1 & -2 \end{bmatrix} $\endgroup$ – nicoguaro Jun 25 '15 at 21:28
  • $\begingroup$ @nicoguaro Yes, but in the notation of the other answer, that's the $A_{11}$ matrix, whereas $A_{12}$ is $$ \begin{pmatrix} 1&0\\0&0\\\vdots&\vdots\\0&0\\0&1 \end{pmatrix}, $$ and the ordering of points is $x_1=(u_{2:N+1})$, $x_2=(u_1, u_{N+2})$. I was talking about $A_{12}$, which is what determines how the boundary conditions go to the r.h.s. In this question the equations are simple enough that this is kind of overkill. My point is the original code handles this correctly. $\endgroup$ – Kirill Jun 25 '15 at 21:32
  • $\begingroup$ It seems that we are saying the same in different ways... or I don't get it. But it's ok, this is not a chat anyway. $\endgroup$ – nicoguaro Jun 25 '15 at 21:35

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.