Stiffness matrix for Galerkin method

Question

I am working on a finite element approximation to solve the following ODE type:

$$-\dfrac{d}{dx}\bigg(a(x)\dfrac{du}{dx}\bigg)=f(x)$$

with $u(0) = u(1) = a(0)$, using a Galerkin method. We start by defining a mesh $x_j=jh$ over $N$ points. We take the piecewise function $\phi$ as a basis:

$$\phi_j(x) =\begin{cases} 0&\text{if}\, ~x\not\in(x_{j-1},x_{j+1})\\\\ \dfrac{x-x_{j-1}}{x_{j}-x_{j-1}}&\text{if}\, x\in (x_{j-1},x_{j})\\\\ \dfrac{x_{j+1}-x}{x_{j+1}-x_j}&\text{if}~x\in (x_j,x_{j+1}) \end{cases}$$

Specifically, for the problem:

We take the test method $v(x)=\phi_i(x)$ and the Galerkin method we are using is:

Find $u$ such that $\displaystyle\int_0^1 a(x)\dfrac{du}{dx}\dfrac{dv}{dx}~dx=\displaystyle\int_0^1 f(x)v(x) dx$, for approximation:

$u=\displaystyle\sum_{j=1}^M\xi_j~\phi_j(x)$

Substituting for u and v into the method, the left hand side gives the stiffness matrix with entries

$(a_{ij})=\displaystyle\int_0^1 a(x)\dfrac{d\phi_j}{dx}\dfrac{d\phi_i}{dx}~dx$

In the specific example I'm doing, $a(x)=1$ and we should find:

$(a_{ij})=\begin{cases} \frac{2}{h}&\text{if}\, ~i=j\\ \frac{-1}{h}&\text{if}\, i=j+1,~\text{or}~i=j-1 \end{cases}$

I don't understand how to get those values for $(a_{ij})$?

My attempt:

$$\phi_j'(x) =\begin{cases} 0&\text{if}\, ~x\not\in(x_{j-1},x_{j+1})\\ \dfrac{1}{h}&\text{if}\, x\in (x_{j-1},x_{j})\\ \dfrac{-1}{h}&\text{if}~x\in (x_j,x_{j+1}) \end{cases}$$

But computing $\phi_i'$ seems to be throwing me.

Answer 1

You are already on right track. Note that after substituting your ansatz for $u$ your bilinearform becomes $$\int_0^1 a(x) v' \sum_i^M \xi_j \phi_j' \mathrm d x$$ Where now we test with all $v \in \text{span} \{\phi_j\}_{j = 1, \dots , M}$ to obtain $$\int_0^1 a(x) \phi_i' \sum_i^M \xi_j \phi_j' \mathrm d x \: \forall \:i$$ Now for each $\phi_i'$ you just need to check which $\phi_j'$ are non-zero on shared intervals. More mathematically, you want to find $$ \phi_j' \colon \text{supp }\phi_i' \cap \text{supp }\phi_j' \neq \emptyset $$ This is done easily if you sketch your basis functions:

Thus, you can clearly see that only basis functions $\phi_{j-1}, \phi_j, \phi_{j+1}$ share support with $\phi_j$. The bilinearform is thus dramatically simplified as $$ \xi_{i-1} \int_{x_{i-1}}^{x_i} a(x) \phi_{i-1}' \phi_i' \mathrm d x + \xi_{i} \int_{x_{i-1}}^{x_{i+1}} a(x) \phi_{i}' \phi_i' \mathrm d x + \xi_{i+1} \int_{x_{i}}^{x_{i+1}} a(x) \phi_{i+1}' \phi_i' \mathrm d x \: \forall \:i$$ For equidistant $x_i - x_{i-1} = h \: \forall \: i$ grids and $a(x) \equiv 1$ this further simplifies to \begin{align} &\xi_{i-1} \int_{x_{i-1}}^{x_i} \frac{-1}{h} \frac{1}{h} \mathrm d x + \xi_{i} \bigg( \int_{x_{i-1}}^{x_{i}} \frac{1}{h} \frac{1}{h} \mathrm d x + \int_{x_{i}}^{x_{i+1}} \frac{-1}{h} \frac{1}{h} \mathrm d x \bigg)+ \xi_{i+1} \int_{x_{i}}^{x_{i+1}} \frac{1}{h} \frac{-1}{h} \mathrm d x \: \forall \:i \\ =& \xi_{i-1} \frac{-1}{h} + \xi_{i} \bigg( \frac{1}{h} + \frac{1}{h} \bigg)+ \xi_{i+1} \frac{-1}{h} \mathrm d x \: \forall \:i \\ =& \xi_{i-1} \frac{-1}{h} + \xi_{i} \frac{2}{h} + \xi_{i+1} \frac{-1}{h} \mathrm d x \: \forall \:i\end{align} The linear system is then obtained by grouping the coefficiencts $\xi_i$ into the $M$-dimensional vector. The condition $\forall \: i$ corresponds to one equation which has to be fulfilled, thus a row in matrix. Consequently, the system reads $$ \frac{1}{h} \begin{pmatrix} 2 & -1 & 0 & 0 & 0 & 0 & 0\\ -1 & 2 & 1 & 0 & 0 & 0 & 0\\ 0 & -1 & 2 & 1 & 0 & 0 & 0\\ 0 & 0 & \ddots & \ddots & \ddots & 0 & 0\\ 0 & 0 & 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & 0 & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & 0 & 0 & -1 & 2 \end{pmatrix} \begin{pmatrix} \xi_1 \\ \xi_2 \\ \xi_3 \\ \vdots \\ \xi_{M-2} \\ \xi_{M-1} \\ \xi_M\end{pmatrix} = \begin{pmatrix} \int_0^{x_1} f \phi_1 \mathrm dx \\ \int_0^{x_2} f \phi_2 \mathrm dx \\ \int_{x_1}^{x_3} f \phi_3 \mathrm dx \\ \vdots \\ \int_{x_{M - 3}}^{x_{M-1}} f \phi_{M-2} \mathrm dx \\ \int_{x_{M - 2}}^{x_M} f \phi_{M-1} \mathrm dx \\ \int_{x_{M - 1}}^{x_M} f \phi_M \mathrm dx \end{pmatrix}$$

where you can see that the siffness matrix fulfills the property you desire.