4
$\begingroup$

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:

Sample FEM Solution - Odd number Elements


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

            Ie(ii,jj) = -a*lgquad(fun,[],-1,1); % Calling Guass-Legendre Quad routine

            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
        b((kk-1)+ii) = b((kk-1)+ii) + lgquad(fun,[],-1,1);
    end
end

Can someone give me an idea of what my issue might be?

$\endgroup$
  • 2
    $\begingroup$ I believe your fundamental problem is that you are trying to prescribe two boundary conditions for a first-order ODE. I modified your code to prescribe a BC only at the left end and it seems to give correct results for both even and odd numbers of elements. However, I don't understand why your code apparently works as-is for an odd number of elements. $\endgroup$ – Bill Greene Apr 10 '16 at 18:13
  • $\begingroup$ @BillGreene, you were right, removing the second Dirichlet BC solved the problem. I'm guessing that I created an overdetermined system by including BCs - I'll remember that for next time. I'm not sure why an odd number of elements works either, but I appreciate that you reminded me of that issue. Thanks again! $\endgroup$ – cbcoutinho Apr 11 '16 at 12:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.