# Boundary conditions in the Finite Element Method at only one side of computational domain

I want to solve a Sturm-Liouville problem in 1D, i.e., \begin{align} [p(x)\ u'(x)]'+q(x)u(x) = f(x) \end{align} with boundary conditions \begin{align} u(0)=a \hspace{1cm} u'(0)=b \end{align} How do I implement these conditions?

What I did:

First I multiply both sides by a test function $v\in\lbrace c\in H^1\big\vert v(0)=0, v'(0)=0 \rbrace$ and integrate over the computational domain, i.e. from 0 to 1 in my case. \begin{align} -\int\limits_0^1\left(p(x)u'(x)\right)'v(x)dx +\int\limits_0^1q(x)u(x)v(x)dx = \int\limits_0^1f (x)v(x)dx \end{align} The first term can be transformed by partial integration. \begin{align} -p(1)u'(1)v(1) + \underset{=0}{\underbrace{p(0)u'(0)v(0)}} +\int\limits_0^1p(x)u'(x)v'(x)dx \end{align} I want to use hermite cubic splines denoted by $\phi_k^l$ with the property \begin{align} \left(\phi_k^l\right)^{(m)}(x_j)=\delta_{kj}\delta_{lm} \end{align} where $(m)$ denotes the 0th or 1st derivative. Inserting this in the weak formulation above we finally arrive at the matrix equation \begin{align} \nonumber -p(1)u^1_{N +1}&\delta_{0l}\delta_{N+1j} +a\int\limits_0^1 dx p(x)\left[\phi_0^0\,'(x)\phi_j^l\,'(x) + q(x)\phi_0^0(x)\phi_j^l(x)\right]\\ & +b\int\limits_0^1 dx\, p(x)\left[\phi_0^1\,'(x)\phi_j^l\,'(x) + q(x)\phi_0^1(x)\phi_j^l(x)\right]\\ &+\sum\limits_{i=1}^{N+1}\sum\limits_{k=0}^1u_i^k \int\limits_0^1 dx\, p(x)\left[\phi_i^k\,'(x)\phi_j^l\,'(x) + q(x)\phi_i^k(x)\phi_j^l(x)\right]\\ =&\int\limits_0^1f (x)\phi_j^l (x)dx \hspace{2cm}j=1,\dots,N+1\ \ l=0,1 \end{align} Is this derivation correct up to here?

The matrix equation is set up by first mapping every element to the standard intervall $[-1,1]$ and then using Gauss Integration.

Afterwards the global matrix is assembled. Finally boundary conditions have to be implemented according to the above equation. This means we delete the lines and columns corresponding to the coeffiicients of $\phi_0^0$ and $\phi_0^1$. Subtracting the second and third term (the ones with a and b in front) from the left hand side of the above equation on the right hand side and subtraction $p(1)$ at the corresponding position of the matrix we arrive at the final matrix equation. This is visualized in the figure Implementation of boundary conditions.

But using the extremely simple case of $p(x)=q(x)=1$, $f(x)=0$ and using only one domain we can find an analytic solution. The comparison between the analytic and the numerical solution is shown in figure comparision of solutions. Thus there must be a mistake in the above implementation process (other boundary conditions provided at both ends of the domain give correct solutions.)

For an increasing number of domains the solution does also not converge. Therefore something is obviously wrong but I do not have any idea.

• You are right, there is a typo in the figure showing the implementation of the boundary conditions, it should be $K_{32}$ and $K_{42}$. For a larger number of elements I am also seeing convergence to the wrong solution. The numerical solution behaves as $u(1) \rightarrow 0$. As this is my first larger numerical project I have no expirience thus I am not sure if that helps, but I also used Spectral Element Methods (Chebyshev Polynomials in several subintervals) for this problem and again see the tendency towards $u(1) \rightarrow 0$. – NoIdea Nov 3 '15 at 11:22
• For a large number of elements the obtained solution (not the correct analytic one) seems to be 0 at the endpoint $x=1$, i.e. the coefficient for the last basis function seems to be tend to 0. That's what I meant by $u(1)\rightarrow 0$. So I do not make the solution to be zero there but the numerical outcome turns to be zero. – NoIdea Nov 3 '15 at 15:52