Finite element method: what are bilinear and linear forms?

Question

Sample problem

Consider this ODE:

$$ u : \Omega \to \Re $$ $$\nabla u-u=0 \qquad \mathrm{in} \qquad \Omega$$ $$\Omega\subset\Re$$ $$ u = 1 \qquad \mathrm{on} \qquad \Gamma$$

• $\Omega\subset\Re$ is the domain of interest where $0\le x \le 1$
• $\Gamma$ is domain's boundary where $x=0$
• $\nabla u-u=0$ would become $\frac{du}{dx}-u=0$

Finite element as I know it

I would solve it by using the finite element method. First, I assume the solution to be a polynomial which satisfies the boundary condition:

$$ u := \begin{bmatrix} 1 & x & x^2 & x^3 & x^4 \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix}$$

So, the equation $\frac{du}{dx}-u=0$ would become:

$$ \frac{d}{dx}(\begin{bmatrix} 1 & x & x^2 & x^3 & x^4 \end{bmatrix}) \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} - \begin{bmatrix} 1 & x & x^2 & x^3 & x^4 \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} = 0$$

I multiply both sides of the question by transpose of $\begin{bmatrix} 1 & x & x^2 & x^3 & x^4 \end{bmatrix}$ like:

$$ \begin{bmatrix} 1 \\ x \\ x^2 \\ x^3 \\ x^4 \end{bmatrix} \left( \frac{d}{dx}(\begin{bmatrix} 1 & x & x^2 & x^3 & x^4 \end{bmatrix}) \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} - \begin{bmatrix} 1 & x & x^2 & x^3 & x^4 \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} = 0 \right)$$

The equation would become:

$$ \begin{bmatrix} 1 \\ x \\ x^2 \\ x^3 \\ x^4 \end{bmatrix} \frac{d}{dx}(\begin{bmatrix} 1 & x & x^2 & x^3 & x^4 \end{bmatrix}) \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} - \begin{bmatrix} 1 \\ x \\ x^2 \\ x^3 \\ x^4 \end{bmatrix} \begin{bmatrix} 1 & x & x^2 & x^3 & x^4 \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

The equation would be:

$$ \begin{bmatrix} 1 \\ x \\ x^2 \\ x^3 \\ x^4 \end{bmatrix} \begin{bmatrix} 0 & 1 & 2x & 3x^2 & 4x^3 \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} - \begin{bmatrix} 1 \\ x \\ x^2 \\ x^3 \\ x^4 \end{bmatrix} \begin{bmatrix} 1 & x & x^2 & x^3 & x^4 \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

The equation would be:

$$ \begin{bmatrix} 0 & 1 & 2x & 3x^2 & 4x^3 \\ 0 & x & 2x^2 & 3x^3 & 4x^4 \\ 0 & x^2 & 2x^3 & 3x^4 & 4x^5 \\ 0 & x^3 & 2x^4 & 3x^5 & 4x^6 \\ 0 & x^4 & 2x^5 & 3x^6 & 4x^7 \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} - \begin{bmatrix} 1 & x & x^2 & x^3 & x^4 \\ x & x^2 & x^3 & x^4 & x^5 \\ x^2 & x^3 & x^4 & x^5 & x^6 \\ x^3 & x^4 & x^5 & x^6 & x^7 \\ x^4 & x^5 & x^6 & x^7 & x^8 \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

Now we do the integral of both sides of the equation over the $\Omega$ domain:

$$ \int_0^1 \left( \begin{bmatrix} 0 & 1 & 2x & 3x^2 & 4x^3 \\ 0 & x & 2x^2 & 3x^3 & 4x^4 \\ 0 & x^2 & 2x^3 & 3x^4 & 4x^5 \\ 0 & x^3 & 2x^4 & 3x^5 & 4x^6 \\ 0 & x^4 & 2x^5 & 3x^6 & 4x^7 \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} - \begin{bmatrix} 1 & x & x^2 & x^3 & x^4 \\ x & x^2 & x^3 & x^4 & x^5 \\ x^2 & x^3 & x^4 & x^5 & x^6 \\ x^3 & x^4 & x^5 & x^6 & x^7 \\ x^4 & x^5 & x^6 & x^7 & x^8 \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} \right) \mathrm{d}x$$

The equation would be:

$$ \int_0^1 \left( \begin{bmatrix} 0 & 1 & 2x & 3x^2 & 4x^3 \\ 0 & x & 2x^2 & 3x^3 & 4x^4 \\ 0 & x^2 & 2x^3 & 3x^4 & 4x^5 \\ 0 & x^3 & 2x^4 & 3x^5 & 4x^6 \\ 0 & x^4 & 2x^5 & 3x^6 & 4x^7 \end{bmatrix} \right) \mathrm{d}x \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} - \int_0^1 \left( \begin{bmatrix} 1 & x & x^2 & x^3 & x^4 \\ x & x^2 & x^3 & x^4 & x^5 \\ x^2 & x^3 & x^4 & x^5 & x^6 \\ x^3 & x^4 & x^5 & x^6 & x^7 \\ x^4 & x^5 & x^6 & x^7 & x^8 \end{bmatrix} \right) \mathrm{d}x \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

Having computed the integrals, the equation would be:

$$ \begin{bmatrix} 0 & 1 & 1 & 1 & 1 \\ 0 & \frac{1}{2} & \frac{2}{3} & \frac{3}{4} & \frac{4}{5} \\ 0 & \frac{1}{3} & \frac{1}{2} & \frac{3}{5} & \frac{2}{3} \\ 0 & \frac{1}{4} & \frac{2}{5} & \frac{1}{2} & \frac{4}{7} \\ 0 & \frac{1}{5} & \frac{1}{3} & \frac{3}{7} & \frac{1}{2} \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} - \begin{bmatrix} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} \\ \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} \\ \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

Having done the subtraction, the equation would be:

$$ \begin{bmatrix} -1 & \frac{1}{2} & \frac{2}{3} & \frac{3}{4} & \frac{4}{5} \\ \frac{-1}{2} & \frac{1}{6} & \frac{5}{12} & \frac{11}{20} & \frac{19}{30} \\ \frac{-1}{3} & \frac{1}{12} & \frac{3}{10} & \frac{13}{30} & \frac{11}{21} \\ \frac{-1}{4} & \frac{1}{20} & \frac{7}{30} & \frac{5}{14} & \frac{25}{56} \\ \frac{-1}{5} & \frac{1}{30} & \frac{4}{21} & \frac{17}{56} & \frac{7}{18} \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

Getting rid of the first row, the equation would become:

$$ \begin{bmatrix} \frac{-1}{2} \\ \frac{-1}{3} \\ \frac{-1}{4} \\ \frac{-1}{5} \end{bmatrix} + \begin{bmatrix} \frac{1}{6} & \frac{5}{12} & \frac{11}{20} & \frac{19}{30} \\ \frac{1}{12} & \frac{3}{10} & \frac{13}{30} & \frac{11}{21} \\ \frac{1}{20} & \frac{7}{30} & \frac{5}{14} & \frac{25}{56} \\ \frac{1}{30} & \frac{4}{21} & \frac{17}{56} & \frac{7}{18} \end{bmatrix} \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

The solution of the above linear algebra problem would be:

$$ \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} \frac{1704}{1709} \\ \frac{882}{1709} \\ \frac{224}{1709} \\ \frac{126}{1709} \end{bmatrix} $$

The 1st row which we got rid of giving:

$$ \begin{bmatrix} -1 & \frac{1}{2} & \frac{2}{3} & \frac{3}{4} & \frac{4}{5} \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} -1 & \frac{1}{2} & \frac{2}{3} & \frac{3}{4} & \frac{4}{5} \end{bmatrix} \begin{bmatrix} 1 \\ \frac{1704}{1709} \\ \frac{882}{1709} \\ \frac{224}{1709} \\ \frac{126}{1709} \end{bmatrix} = \frac{-1}{8545} = -0.000117 \approx 0 $$

Comparison with the exact solution $e^x$ is done:

x	$e^x$	FEM approximation
0	1	1
0.1	1.10517	1.10501
0.2	1.2214	1.22123
0.3	1.34986	1.34971
0.4	1.49182	1.49168
0.5	1.64872	1.64855
0.6	1.82212	1.8219
0.7	2.01375	2.0135
0.8	2.22554	2.22526
0.9	2.4596	2.45932
1.0	2.71828	2.71796

Is my understanding correct?

My first question

Is my finite element method of solving the above ODE correct? Is it thorough?

Bilinear form and linear form

my second question

I'm studying this paper: https://arxiv.org/pdf/1911.09220.pdf

It contains this explanation about its Finite Element Abstractions:

I'm trying to find a correspondence between paper equations (3), (4), (5), and (6) and my finite element understanding. The paper is talking about bilinear and linear forms. I'm trying to figure out the corresponding bilinear and linear forms inside my finite element solution described above. I'm trying to wrap my head around the paper equations with my limited knowledge of the finite element method. I'd appreciate any help or hint that would make me understand the paper.

Answer 1

Basically starting with the equation:

$$-\Delta u = f$$

we multiply both sides by a test function $v$ to get

$$-v\Delta u = fv$$

and then we integrate over the domain $\Omega$

$$-\int_\Omega v\Delta u = \int_\Omega fv$$

and now we use integration by parts on the left hand side (there are regularity conditions on $u$ and $v$ to ensure you can do this)

$$\int_\Omega \nabla v \cdot \nabla u = -\int_\Omega v\Delta u$$

and now we recognize this as a bilinear form on the left hand side and a linear form on the right hand side.

How does integration by parts help?

The classical integration by parts formula is the following: $$\int_a^b u'(x)v(x)dx = u(b)v(b)-u(a)v(a) - \int_a^bu(x)v'(x)dx$$ If we apply it to a second derivative, we get: $$\int_a^b u''(x)v(x)dx = u'(b)v(b)-u'(a)v(a) - \int_a^bu'(x)v'(x)dx$$ Now, if $v(b)=v(a)=0$ then we have the one-dimensional version of what I wrote above $$\int_\Omega \nabla v\cdot \nabla u = \int_\Omega fv$$ $$\int_a^b v'u' = -\int_a^b v u''$$ The fact that we take $v$ to be a test function means that it vanishes on the boundary of $\Omega$ (in the one-dimensional version, this means $v(a)=v(b)=0$ because the boundary is $\{a,b\}$).

To apply this to a PDE you just need the higher-dimensional version of integration by parts, i.e., Green's first identity.

Answer 2

Basis functions

I tried to solve the same ODE with the help of basis functions which are like basis vectors. They are nonzero only on a small portion inside the domain and zero at the rest.

In 2D, the basis function would look like this taken from here:

So, in 1D, I created 5 piecewise basis functions like below to cover the domain of interest $\Omega$ i.e. $ 0\le x \le 1 $:

The general definition for a piecewise basis function is this:

$$ \Psi = \left\{ \begin{array}{rcl} 0 & \mbox{for} & x<x_0-\delta \\ \frac{1}{\delta}x-\frac{x_0-\delta}{\delta} & \mbox{for} & x_0 - \delta \leq x < x_0 \\ \frac{-1}{\delta}x+\frac{x_0+\delta}{\delta} & \mbox{for} & x_0 \leq x < x_0 + \delta \\ 0 & \mbox{for} & x_0 + \delta \leq x \end{array}\right. $$

I set $ \delta = 0.25 $ and $x_0 = 0, 0.25, 0.5, 0.75, 1.0 $ for five basis functions $ \Psi_1, \Psi_2, \Psi_3, \Psi_4, \Psi_5 $

Solving ODE

Then I set:

$$ u := \begin{bmatrix} \Psi_1 & \Psi_2 & \Psi_3 & \Psi_4 & \Psi_5 \end{bmatrix} \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} $$

This $ u $ function satisfies the boundary condition. We have $ u = 1 \: \mathrm{on} \: \Gamma $ i.e. $ x = 0 $

$ \begin{bmatrix} dof \end{bmatrix} $ are the degrees of freedom i.e. unknowns:

$$ \begin{bmatrix} \Psi \end{bmatrix} = \begin{bmatrix} \Psi_1 & \Psi_2 & \Psi_3 & \Psi_4 & \Psi_5 \end{bmatrix} \qquad \mbox{and} \qquad \begin{bmatrix} dof \end{bmatrix} = \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} $$

Let's simply write:

$$ u := \begin{bmatrix} \Psi \end{bmatrix} \begin{bmatrix} dof \end{bmatrix} $$

Then the equation would be:

$$ \frac{du}{dx}-u=0 $$

$$ \frac{d}{dx}(\begin{bmatrix} \Psi \end{bmatrix}) \begin{bmatrix} dof \end{bmatrix} - \begin{bmatrix} \Psi \end{bmatrix} \begin{bmatrix} dof \end{bmatrix} = 0 $$

Multiplied by the transpose of the basis functions:

$$ \begin{bmatrix} \Psi \end{bmatrix} ^T \frac{d}{dx}(\begin{bmatrix} \Psi \end{bmatrix}) \begin{bmatrix} dof \end{bmatrix} - \begin{bmatrix} \Psi \end{bmatrix} ^T \begin{bmatrix} \Psi \end{bmatrix} \begin{bmatrix} dof \end{bmatrix} = \begin{bmatrix} \Psi \end{bmatrix} ^T \; 0 $$

$$ (\begin{bmatrix} \Psi \end{bmatrix} ^T \frac{d}{dx}(\begin{bmatrix} \Psi \end{bmatrix}) - \begin{bmatrix} \Psi \end{bmatrix} ^T \begin{bmatrix} \Psi \end{bmatrix}) \; \begin{bmatrix} dof \end{bmatrix} = \begin{bmatrix} 0 \end{bmatrix} $$

Integrated over the domain of interest:

$$ \int_0^1 (\begin{bmatrix} \Psi \end{bmatrix} ^T \frac{d}{dx}(\begin{bmatrix} \Psi \end{bmatrix}) - \begin{bmatrix} \Psi \end{bmatrix} ^T \begin{bmatrix} \Psi \end{bmatrix}) dx \; \begin{bmatrix} dof \end{bmatrix} = \int_0^1 \begin{bmatrix} 0 \end{bmatrix} dx $$

Having integrated from 0 to 1, the linear algebra would become:

$$ ( \begin{bmatrix} -0.5 & 0.5 & 0 & 0 & 0 \\ -0.5 & 0 & 0.5 & 0 & 0 \\ 0 & -0.5 & 0 & 0.5 & 0 \\ 0 & 0 & -0.5 & 0 & 0.5 \\ 0 & 0 & 0 & -0.5 & 0.5 \end{bmatrix} - \begin{bmatrix} 0.083333 & 0.041667 & 0 & 0 & 0 \\ 0.041667 & 0.166667 & 0.041667 & 0 & 0 \\ 0 & 0.041667 & 0.166667 & 0.041667 & 0 \\ 0 & 0 & 0.041667 & 0.166667 & 0.041667 \\ 0 & 0 & 0 & 0.041667 & 0.083333 \end{bmatrix} ) \begin{bmatrix} 1 \\ a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

The final solution would be:

$$ \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} 1.2418 \\ 1.6334 \\ 2.0616 \\ 2.6800 \end{bmatrix} $$

Comparison

Comparison with the exact solution $e^x$ is done:

x	$e^x$	FEM approximation
0	1	1
0.1	1.10517	1.09673
0.2	1.2214	1.19345
0.3	1.34986	1.32013
0.4	1.49182	1.47676
0.5	1.64872	1.63339
0.6	1.82212	1.80466
0.7	2.01375	1.97593
0.8	2.22554	2.18526
0.9	2.4596	2.43265
1.0	2.71828	2.68004