Finite Element Method for a Two-Point Problem

Question

Given the following two-point problem: $$-y''(x)+(by)'(x)=f(x), \forall x \in [0,1]$$ $$y(0)=0, y'(1)=my(1)$$ where $ b \in C^1([0,1];\Bbb{R}), f \in C([0,1];\Bbb{R})$ and $ m \in \Bbb{R}$ a constant. What could be a finite element method for the construction of the approximation of the solution $y$ of the problem above, where the finite element space consists of continuous and partially linear functions?

Answer 1

As usual, we start by finding the variational form of the problem.

So, let's multiply the ODE by a test function $v$ with $v(0)=0$ (because we have a homogeneous boundary condition at $x=0$) and integrate over $[0,1]$ to get $$\int\limits_0^1 { - y''v} + \int\limits_0^1 {by'v} = \int\limits_0^1 {fv}.$$ Using integration by parts on the first integral and bearing in mind that $v(0)=0$, we get $$\int\limits_0^1 { - y''v} = \left[ { - y'v} \right]_{x = 0}^{x = 1} + \int\limits_0^1 {y'v'} = - y'\left( 1 \right)v\left( 1 \right) + \int\limits_0^1 {y'v'} = - my\left( 1 \right)v\left( 1 \right) + \int\limits_0^1 {y'v'}$$ Substituting this result in the previous expression leads to $$\int\limits_0^1 {y'v'} + \int\limits_0^1 {by'v} = \int\limits_0^1 {fv} + my\left( 1 \right)v\left( 1 \right).$$ Thus, the continuous variational formulation reads as follows:

Find $u \in {V^0} = \left\{ {\omega :\left\| \omega \right\| + \left\| {\omega '} \right\| < \infty ,\omega \left( 0 \right) = 0} \right\}$ such that $$\int\limits_0^1 {y'v'} + \int\limits_0^1 {by'v} = \int\limits_0^1 {fv} + my\left( 1 \right)v\left( 1 \right),{\rm{ }}\forall v \in {V^0}.$$

To get the FEM formulation, let's consider a uniform partition $T_h$ of $[0,1]$: $0=x_0 < x_1 < \ldots < x_{N+1}=1$ into the subintervals $I_n=[x_n,x_{n-1}]$, such that $x_n-x_{n-1}=h$, $n=1,\ldots,N+1$.

The finite element method is now formulated as follows:

Find $u \in V_h^0 = \left\{ {{\omega _h}:{\omega _h}{\textrm{ is piecewise linear, continuous on }}{T_h},{\omega _h}\left( 0 \right) = 0} \right\}$ such that $$\int\limits_0^1 {{{y_h}^\prime}{{v'}_h}} + \int\limits_0^1 {b{{y_h}^\prime}{v_h}} = \int\limits_0^1 {f{v_h}} + m{y_h}\left( 1 \right){v_h}\left( 1 \right),{\rm{ }}\forall {v_h} \in V_h^0$$

Let's consider the 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 ,N$$ plus the semi-hat function: $${\varphi _{N + 1}}\left( x \right) = \left\{ {\begin{array}{*{20}{c}}{\frac{{x - {x_{j - 1}}}}{h}}&{,x \in [{x_N},{x_{N + 1}}]}\\0&{,x \notin [{x_N},{x_{N + 1}}]}\end{array}} \right. .$$ Now, using the ansatz ${y_h} = \sum\limits_{j = 1}^{N + 1} {{\xi _j}{\varphi _j}}$ and inserting this expression on the FEM formulation we get $$\int\limits_0^1 {\left( {\sum\limits_{j = 1}^{N + 1} {{\xi _j}{\varphi _j}^\prime } } \right){\varphi _i}^\prime } + \int\limits_0^1 {b\left( {\sum\limits_{j = 1}^{N + 1} {{\xi _j}{\varphi _j}^\prime } } \right){\varphi _i}} = \int\limits_0^1 {f{\varphi _i}} + m{\varphi _{N + 1}}\left( 1 \right){\varphi _i}\left( 1 \right),{\rm{ }}i = 1, \ldots ,N + 1$$ or, more specifically, we can formulate the FEM as follows:

For $i=1,\ldots,N+1$, find $\xi _j$ from the following system of linear equations $$\sum\limits_{j = 1}^{N + 1} {\left( {\int\limits_0^1 {{\varphi _j}^\prime {\varphi _i}^\prime } } \right){\xi _j}} + \sum\limits_{j = 1}^{N + 1} {\left( {\int\limits_0^1 {b{\varphi _j}^\prime {\varphi _i}} } \right){\xi _j}} = \int\limits_0^1 {f{\varphi _i}} + m{\varphi _{N + 1}}\left( 1 \right){\varphi _i}\left( 1 \right)$$ or, equivalently, $$A\xi = b,$$ where $A$ is a tridiagonal matrix (due to the local support of the hat functions) with entries:

${\left( A \right)_{i,i}} = \int\limits_{{x_{i - 1}}}^{{x_{i + 1}}} {{{\left( {{\varphi _i}^\prime } \right)}^2}} + \int\limits_{{x_{i - 1}}}^{{x_{i + 1}}} b {\varphi _i}^\prime {\varphi _i} = \frac{2}{h} + \frac{1}{{{h^2}}}\int\limits_{{x_{i - 1}}}^{{x_i}} b \left( {x - {x_{i - 1}}} \right) - \frac{1}{{{h^2}}}\int\limits_{{x_i}}^{{x_{i + 1}}} b \left( {{x_{i + 1}} - x} \right)$

${\left( A \right)_{i,i + 1}} = \int\limits_{{x_{i - 1}}}^{{x_{i + 1}}} {{\varphi _{i + 1}}^\prime {\varphi _i}^\prime } + \int\limits_{{x_{i - 1}}}^{{x_{i + 1}}} b {\varphi _{i + 1}}^\prime {\varphi _i} = - \frac{1}{h} + \frac{1}{{{h^2}}}\int\limits_{{x_i}}^{{x_{i + 1}}} b \left( {{x_i} - x} \right)$

${\left( A \right)_{i-1,i}} = \int\limits_{{x_{i - 1}}}^{{x_{i + 1}}} {{\varphi _{i - 1}}^\prime {\varphi _i}^\prime } + \int\limits_{{x_{i - 1}}}^{{x_{i + 1}}} b {\varphi _{i - 1}}^\prime {\varphi _i} = - \frac{1}{h} - \frac{1}{{{h^2}}}\int\limits_{{x_i}}^{{x_{i + 1}}} b \left( {x - {x_{i - 1}}} \right)$

and the vector $b$ is as follows: $$b = \left[ {\begin{array}{*{20}{c}}{\int\limits_0^1 {f{\varphi _1}} }\\{\int\limits_0^1 {f{\varphi _2}} }\\ \vdots \\{\int\limits_0^1 {f{\varphi _N}} }\\{\int\limits_0^1 {f{\varphi _{N + 1}}} + m}\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{\frac{1}{h}\int\limits_{{x_0}}^{{x_1}} {f\left( {x - {x_1}} \right)} + \frac{1}{h}\int\limits_{{x_1}}^{{x_2}} {f\left( {{x_1} - x} \right)} }\\{\frac{1}{h}\int\limits_{{x_1}}^{{x_2}} {f\left( {x - {x_2}} \right)} + \frac{1}{h}\int\limits_{{x_2}}^{{x_3}} {f\left( {{x_2} - x} \right)} }\\ \vdots \\{\frac{1}{h}\int\limits_{{x_{N - 1}}}^{{x_N}} {f\left( {x - {x_N}} \right)} + \frac{1}{h}\int\limits_{{x_N}}^{{x_{N + 1}}} {f\left( {{x_N} - x} \right)} }\\{\frac{1}{h}\int\limits_{{x_N}}^{{x_{N + 1}}} {f\left( {x - {x_{N + 1}}} \right)} + m}\end{array}} \right]$$

Please note that, in order to implement this method, we would have to approximate all the integrals numerically, most certainly using some quadrature rule, e.g. the Simpson's formula which states that $$\int\limits_a^b {g\left( x \right)} dx \approx \frac{{g\left( a \right) + 4g\left( {\frac{{a + b}}{2}} \right) + g\left( b \right)}}{6} \cdot \left( {b - a} \right).$$