Intro:

I am developing a Galerkin FEM code in matlab - starting small with a simple 1D ODE. The equation I'm trying to solve is: $$a u_x = cos(x), x \in [0, 2\pi]$$ Which has a known exact solution $$u = sin(x)/a$$

When I solve the equation using my code on an odd number of linear elements (nodes = N+1), the solution is as expected:

Problem

I'm running into a problem when I set the number of elements to an even number, in which case the solution vector is just filled with NaNs and I receive the following warning/error:

Warning: Matrix is singular to working precision

Most of the questions I've looked up on mathworks answers have to do with incorrectly constructed matrices - that's not the case here. I'm not sure if I'm running into a numerical issue, where a zero is too small, or if this is an implementation issue.

I construct both the elemental stiffness matrices and the RHS vector in the same way, and I can't tell a difference between the two systems from a quick comparison of their respective Ax=b structures.

Methodology

I'm generating the elemental stiffness matrices individually, where N is a basis function. The weak form of the ODE above for an elemental stiffness matrix is: $$I_e = \int_\Omega \underline{N} a u_x dx = \int_\Omega \underline{N} cos(x) dx$$

Using integration by parts (i.e. magic), this equation turns into this: $$(au) \underline{N} |_{\partial \Omega L}^{\partial \Omega R} - \int_\Omega (au) \underline{N}_x dx$$

Then I substitute a continuous value of u for the weighted sum of node values ($\hat{u}$), as such: $$u = N_1 \hat{u}_1 + N_2 \hat{u}_2 = \underline{N}^\intercal \underline{\hat{u}}$$

Which is finally equal to the following:

$$I_e = a \underline{N}^\intercal \underline{\hat{u}} \underline{N} |_{\partial \Omega L}^{\partial \Omega R} - a \int_\Omega \underline{N}^\intercal \underline{\hat{u}} \underline{N}_x dx$$

Since the first term in $I_e$ cancels out except for the boundaries, in which case $u$ is already 0, I'm assuming that term can be safely ignored. The resulting elemental stiffness takes the form:

$$I_e = \begin{bmatrix} -a \int_\Omega N_1 N_{1x} dx & -a \int_\Omega N_1 N_{2x} dx \\ -a \int_\Omega N_2 N_{1x} dx & -a \int_\Omega N_2 N_{2x} dx \end{bmatrix}$$

Code

The entire code is in a repository on my github here, but the most important part is here below. All of the functions called from this script have been checked and verified to work properly, namely basis_1D(), lgquad(), and lgwt(). Feel free to check the repo if you're curious

The portion of the code that constructs A and b is:

for kk = 1:numElem
% Initialize elemental stiffness matrix
Ie  = zeros(2);

% Anonymous function that maps natural coordinate to global coordinate
X   = @(s)basis_1D(1,1,s,0)*coordPairs(kk,1) + ...
basis_1D(1,2,s,0)*coordPairs(kk,2); % X = x(s)

% Anonymous function that calculates the Jacobian as a function of the natural coordinate
J   = @(s)basis_1D(1,1,s,1)*coordPairs(kk,1) + ...
basis_1D(1,2,s,1)*coordPairs(kk,2); % J = dx/ds

for ii = 1:2
for jj = 1:2
% Anonymous function to what I denote above as N_x
fun1    = @(s) basis_1D(1,ii,s,1)./J(s);

% Anonymous function to what I denote above as N
fun2    = @(s) basis_1D(1,jj,s,0);

% The two functions combined together, as well as the jacobian for numerical quadrature
fun     = @(s) fun1(s)*fun2(s).*J(s);

II = (kk-1)+ii;
JJ = (kk-1)+jj;
A(II,JJ) = A(II,JJ) + Ie(ii,jj); % Assembling global matrix
end

fun1    = @(s) cos(X(s));
fun2    = @(s) basis_1D(1,ii,s,0);
fun     = @(s) fun1(s)*fun2(s).*J(s); % Combined function for RHS

% Calling Gauss-Legendre Quad routine and assembling RHS