How to solve a boundary value problem with Neumann Boundary conditions using the Finite Element Method

114 Views Asked by At

I have been given the question:

Consider the boundary value problem: $$-u'' + u = f(x) \forall x \in (0,1), u(0) = 0, u'(1) = 0.$$ Using continuous piecewise linear basis functions on a uniform mesh of size h = $\frac{1}{N}$ , N $\geq$ 2 write down the finite element formulation of this problem and show that it has a unique solution $u_h$. Expand $u_h$ in terms of the finite element basis functions $\phi_i$,where $\phi_i(x) = (1 - \frac{|x - x_i|}{h})_+$, i = 1,...,N, by writing $$u_h(x) = \sum_{i=1}^N U_i\phi_i(x)$$ to obtain a system of linear equations for the vector unknowns ($U_1$,...,$U_N$ )$^T$. [Hint: For the Neumann boundary condition u'(1) = 0 you do not need to enforce anything on the function space]

So far I have managed to introudce a test function v satisfying the same boundary conditions as u to obtain:

$$\int_0^1{(-u'' + u)v}dx = \int_0^1{fv}dx$$

which given the boundary conditions and through integration by parts I have shown to be

$$\int_0^1{(u'v' + uv)}dx = \int_0^1{fv}dx$$

but I'm not really sure how to continue

Edit: I have since defined the $L^2$-inner product $\langle u,v \rangle = \int_{I} uvdx$ and introduced the space $$V_h = \{v_h \in C[0,1]:v_h(0) = v_h'(1) = 0, v_h|I_j \: is \: linear \: \forall j = 1,...,M+1\}$$ so I know I need to find $u_h \in V_h$ satisfying: $\langle u_h',v_h'\rangle + \langle u_h,v_h\rangle = \langle f,v_h\rangle$. Then using $u_h(x) = \sum_{i=1}^N U_i\phi_i(x)$, and setting $v_h = \phi_j$ the problem can be rewritten as $$\sum_{i=1}^N U_i(\ \langle \phi_i',\phi_j'\rangle + \langle \phi_i,\phi_j\rangle \ ) = \langle f,\phi_j\rangle$$ which is equivalent to a linear system AU = F with $$[A]_{ij} = \langle \phi_i',\phi_j'\rangle + \langle \phi_i,\phi_j\rangle$$ $$[F]_{j} = \langle f,\phi_j\rangle$$ I am now stuck trying to construct $[A]_{ij}$ and thus proving it is invertible to show the existence of U

1

There are 1 best solutions below

1
On

The basic issue you have here pertains to the weak formulation of the original problem, not really FEM specifically. The way that goes is to start with a candidate solution function $u$ and a test function $v$ without pinning down what function space each is in, multiply both sides by $v$, and integrate over the whole domain. Then you integrate by parts, and keep the boundary terms (temporarily), before analyzing how to be rid of them afterwards.

So in this problem you start with

$$\int_0^1 -u'' v + u v dx = \int_0^1 fv dx$$

and IBP gives

$$-u'(1) v(1) + u'(0) v(0) + \int_0^1 u' v' + uv dx = \int_0^1 fv dx.$$

The weird trick with Neumann boundary conditions is to simply replace $u'(1)$ in the equation with what the BVP wants it to be (so $0$ in your case), rather than forcing $u'(1)$ to actually be $0$. There are a bunch of ways to support why this is the right way to do things, but as a newcomer it is easier to just think of it as a weird trick.

On the flip side, to get rid of the second boundary term, we force $v(0)=0$. This is OK because we are going to also force $u(0)=0$, so we lose one degree of freedom in both spaces.

So the resulting weak formulation is

$$\int_0^1 u' v' + uv dx = \int_0^1 fv dx$$

with the boundary condition $u(0)=v(0)=0$.

Now as for the function spaces, we pose the original problem with $u$ and $v$ in $H^1$ (so that the bilinear form will be continuous) and subject to the above boundary condition. Just defining $H^1$ and proving that the Dirichlet boundary condition makes sense for elements of it requires a detour into functional analysis. Fortunately, you don't need to pay too much attention to this part when you're just starting out with FEM.

The beauty of the FEM theory in this setting is that the Lax-Milgram theorem applies with exactly the same math in the finite dimensional FEM system as it does in the original problem. You check that the bilinear form $B(u,v)= \int_0^1 u' v' + uv dx$ is bounded and coercive, then the FEM version of that bilinear form is automatically also bounded and coercive (because you're just doing a domain restriction). So now, given a subspace $U$ of $H^1$ with this boundary condition, Lax-Milgram tells you that there is a unique $u \in U$ such that $B(u,v)=\int_0^1 P(f) v \, dx$ for all $v \in U$, where $P(f)$ is the projection of $f$ onto $U$.

As an aside, it may help a little bit with new vocabulary overload to know that "coercive" in this context just means that $A$ is positive definite.