Why does my implementation of the galerkin method yield the correct values, but with inverted sign

Question

I tried to implement the galerkin-method for the one dimensional poisson equation. I chose for the boundary conditions $$ u(0) = u(1) = 0 $$ and for the "source" term $f(x)$ $$ \frac{\partial^2 u(x)}{\partial x^2}=f(x) $$ $$f(x)=1$$ Now the problem is the solution, that I expected was $\frac{1}{2}*(x-0.5)^2-0.125$, instead I got the mirrored function around the x-axis, as if had picked $f(x)=-1$.

Below is my code where I generate the solution:

size = 1; % length of the interval
N = 10; % number of elements
K = N+1; % number of knots
h = size / N; % lenght of one element

A = zeros(K);
b = zeros(K, 1);
%initialize "stiffness" matrix and rhs 
for i = 1:K
    for j = 1:K
        if(i==j)
            A(i,j) = 2 / h;
        elseif(abs(i-j)==1)
            A(i,j) = -1 / h;
        end
    end
end
for i = 1:K
    b(i) = h;
end
%remove first and last rows and columns, due to boundary condititons = 0
A(:,1) = [];
A(1,:) = [];
A(:,K-1) = [];
A(K-1,:) = [];
b = b(2:end-1)

%solve and add boundary values
y = [0;linsolve(A, b);0]
x = 0:h:1;

plot(x,y);

The values generated are right apart from the sign, is it an error produced by me or is it maybe even the correct solution, given my input conditions? I would be happy to hear any suggestions.

Answer 1

The equation which you attempt to solve in the code is $-u''=f$. Take care of the right signs when you do integration by parts to obtain the stiffness matrix.

In fact, the weak formulation of your PDE would read: $$a(u,v)=(f,v)$$

where $a(u,v)=-\int_0^1 u'v'$. If you put $u=v$ you see that the diagonal of your stiffness matrix should be negative, whereas yours is positive.

Summing up, your stiffness matrix has the wrong sign.

Answer 2

An easier test: The finite difference approximation of the second derivative is $$ u''(x)=\frac{u(x-h)-2u(x)+u(x+h)}{h^2}, $$ which also indicates that your matrix elements have the wrong sign.

---

With piecewise linear test and kernel functions, the derivatives are the same as (some) finite difference approximation. The difference is in the right side, for Galerkin you get $\frac16(f(x-h)+4f(x)+f(x+h))$ instead of $f(x)$. Here with a constant $f$ both amount to the same.