I am trying to use FEM to approximate the solution to the following BVP: $-\frac{d}{dx}[a(x)u'(x)]+b(x)u(x)=f(x)$, on [0,1] where $u(0)=0$ and $u(1)=1$.
I am using the Galerkin method with hat functions:
${\varphi _j}$$\left( x \right)$ = $\left\{ {\begin{array}{*{20}{c}}{\frac{{x - {x_{j - 1}}}}{h}}&{,x \in [{x_{j - 1}},{x_j}]}\\{\frac{{{x_{j + 1}} - x}}{h}}&{,x \in [{x_j},{x_{j + 1}}]}\\0&{,x \notin [{x_{j - 1}},{x_{j + 1}}]}\end{array}} \right.$ , j=1,$\cdots$ ,m.
To start I have set $u(x)=\sum_{j=1}^{m} \delta_j\varphi_j(x)$, $ $ $\rho(x)=-\frac{d}{dx}[a(x)u'(x)]+b(x)-f(x)$, used the product rule to evaluate $-\frac{d}{dx}[a(x)u'(x)]$, and required $\int_{0}^{1}\rho(x)\varphi_k(x)dx = 0, k=1,\cdots,m$.
After some simplifications this gives me the following integral: $\int_{0}^{1}\sum_{j=1}^{m}\delta_j[-a'(x)\varphi_j'(x)-a(x)\varphi_j''(x)+b(x)]$ $\varphi_k(x)dx-\int_{0}^{1}f(x)\varphi_k(x)dx = 0$
where the $\delta_j$ are the weights we are solving for.
I then simplified the RHS using integration by parts as follows: $-\int_{0}^{1}a(x)\varphi_j''(x)\varphi_k(x)dx = -a(x)\varphi_k(x)\varphi_j'(x)\Big|_{0}^{1}+\int_{0}^{1}a'(x)\varphi_j'(x)\varphi_k(x)dx + \int_{0}^{1}a(x)\varphi_j'(x)\varphi_k'(x)dx$ (***)
From the boundary conditions we know $-a(x)\varphi_k(x)\varphi_j'(x)\Big|_{0}^{1}=0$.
After plugging these two results into (***) above I obtained: $\int_{0}^{1}\sum_{j=1}^{m}\delta_j[a(x)\varphi_j'(x)\varphi_k'(x)+b(x)\varphi_j(x)\varphi_k(x)]dx-\int_{0}^{1}f(x)\varphi_k(x)dx = 0$
Is there a mistake somewhere? My professor was in a rush to cover FEM before the semester ended and he skipped the intermediate steps I have attempted to fill in. His result was not the same as mine and I do not see an error. The semester has ended so I can't go to his office hours for help but I was not able to implement the method and I think my derivation for the matrix A is the cause.
I apologize if my post is a little hard on the eyes. I fought with Latex for over an hour on this post! I've never used it before.
I think your way of arriving at the resulting matrix equation is... strange (read: non-standard and, hence, possibly leading to false conclusions). However, the resulting set of equations seems to be almost correct: the only mistake I could spot is that you are not handling the nonhomogeneous Dirichlet condition at $x=1$ in any way.
Using your notation, let us define (the symmetric) matrices $A, B \in \mathbb{R}^{m \times m}$ and vectors $f, u \in \mathbb{R}^m$ as $$ \begin{align*} A_{ij}&=\int_0^1 a(x)\varphi_i' \varphi_j'\,\mathrm{d}x, \\ B_{ij}&=\int_0^1 b(x) \varphi_i \varphi_j\,\mathrm{d}x, \\ f_i &= \int_0^1 f(x)\varphi_i \,\mathrm{d}x, \\ u_i &= \delta_i. \end{align*} $$ Suppose that you included all mesh points (also $x=0$ and $x=1$) in your matrices while assembling and let $i=m$ correspond to $x=1$ and $i=1$ correspond to $x=0$. Then your matrix equation reads $$ (A+B)u=f $$ with the extra requirement $u_1=0$ and $u_m=1$. Then you can write $$ \begin{bmatrix} u_1 \\ u_2 \\ \vdots \\ u_{m-1} \\ u_m \end{bmatrix} = \begin{bmatrix} 0 \\ u_2 \\ \vdots \\ u_{m-1} \\ 1 \end{bmatrix} = \begin{bmatrix} 0 \\ u_2 \\ \vdots \\ u_{m-1} \\ 0 \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ 1 \end{bmatrix} $$ and, hence, $$ (A+B)\begin{bmatrix} 0 \\ u_2 \\ \vdots \\ u_{m-1} \\ 0 \end{bmatrix} = f - (A+B)\begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ 1 \end{bmatrix}. $$ Finally you may simply drop the first and the last equations in the system since they are automatically satisfied.