Finite Element Method in 1D: Poisson Equation

In this tutorial, we demonstrate how to solve the 1D Poisson equation using the Finite Element Method (FEM). We will cover the following topics:

1. Creating a mesh
2. Defining the source term
3. Constructing local stiffness and mass matrices
4. Assembling the global stiffness and mass matrices
5. Assembling the global load vector
6. Solving the linear system
7. Plotting the solution
8. Calculating the L² error

Problem Statement

We aim to solve the Poisson equation:

$$
-\Delta U + U = 2x
$$

with homogeneous boundary conditions:

$$
U(0) = U(1) = 0
$$

Creating a Mesh

We first define the interval (a, b) and the number of elements in the mesh. We calculate the mesh size (element length) and create the mesh (node coordinates).

# Define the interval (a, b)
a = 0
b = 1

# Define the number of elements
num_elements = 64

# Define the number of nodes
num_nodes = num_elements + 1

# Calculate the mesh size (element length)
h = (b - a) / num_elements

# Create the mesh (node coordinates)
mesh = np.linspace(a, b, num_nodes)

Defining the Source Term

The source term for our problem is defined as:

$$
f(x)=2 * x
$$

# Define the source term function
def f(x):
    return 2 * x

Local Stiffness and Mass Matrices

We create local stiffness and mass matrices for each element. The local stiffness matrix is computed using the element's mesh size, while the local mass matrix uses both the mesh size and the element's index.

# Local stiffness matrix K_e for element e
def local_stiffness_matrix(e):
    assert 0 <= e < num_elements, "e must be between 0 and num_elements - 1"

    K_e = np.zeros((2, 2))

    # Diagonal entries
    K_e[0, 0] = 1 / h
    K_e[1, 1] = 1 / h

    # Off-diagonal entries
    K_e[0, 1] = -1 / h
    K_e[1, 0] = -1 / h

    return K_e

# Local mass matrix M_e for element e
def local_mass_matrix(e):
    assert 0 <= e < num_elements, "e must be between 0 and num_elements - 1"

    M_e = np.zeros((2, 2))

    # Diagonal entries
    M_e[0, 0] = h / 3
    M_e[1, 1] = h / 3

    # Off-diagonal entries
    M_e[0, 1] = h / 6
    M_e[1, 0] = h / 6

    return M_e

Assembling the Global Stiffness and Mass Matrices

We then assemble the global stiffness and mass matrices by iterating through all the elements and adding the corresponding local matrices.

# Assemble the global stiffness matrix
def assemble_global_stiffness_matrix():
    K_global = np.zeros((num_nodes, num_nodes))

    for e in range(num_elements):
    K_e = local_stiffness_matrix(e)
    K_global[e:e+2, e:e+2] += K_e
    # Apply homogeneous boundary conditions by eliminating the first and last rows and columns
    K_global = K_global[1:-1, 1:-1]

    return K_global

Assemble the global mass matrix

def assemble_global_mass_matrix():
M_global = np.zeros((num_nodes, num_nodes))

for e in range(num_elements):
    M_e = local_mass_matrix(e)
    M_global[e:e+2, e:e+2] += M_e

# Apply homogeneous boundary conditions by eliminating the first and last rows and columns
M_global = M_global[1:-1, 1:-1]

return M_global

Assembling the Global Load Vector

We create a local load vector for each element and assemble the global load vector by iterating through all the elements and adding the corresponding local load vectors.

# Local load vector f_e for element e

def local_load_vector(element):
    assert 0 <= element < num_elements, "element must be between 0 and num_elements - 1"

    f_e = np.zeros(2)

    # Load vector entries
    f_e[0] = f(mesh[element]) * h / 2
    f_e[1] = f(mesh[element + 1]) * h / 2

    return f_e

# Assemble the global load vector

def assemble_global_load_vector():
    f_global = np.zeros(num_nodes)

    for element in range(num_elements):
        f_e = local_load_vector(element)
        f_global[element:element+2] += f_e

    # Apply homogeneous boundary conditions by eliminating the first and last entries
    f_global = f_global[1:-1]

    return f_global

Solving the Linear System

We solve the linear system to obtain the solution U.

# Assemble the load vector
f_global = assemble_global_load_vector()

# Solve the linear system
U = np.zeros(num_nodes)
U[1:-1] = np.linalg.solve(K_global + M_global, f_global)

print("Solution U:")
print(U)

Plotting the Solution

We plot the finite element solution alongside the exact solution.

# the exact solution
def u(x):
    C2 = -2 / (np.exp(1) - np.exp(-1))
    C1 = -C2
    return C1 * np.exp(-x) + C2 * np.exp(x) + 2 * x

X = np.linspace(a, b, 100)
exact_solution = u(X)
u = u(mesh)

# Plot the solution
plt.plot(mesh, U, marker='o', label='FEM solution')
plt.plot(X, exact_solution,label='Exact solution')
plt.xlabel('x')
plt.ylabel('U(x)')
plt.title('Finite Element Solution')
plt.legend()
plt.grid()
plt.show()

Calculating the L² Error

We compute the L² error using the mass matrix.

def mass_matrix_l2_error(U, u, M_global):
    # Calculate the difference between the exact solution and the finite element solution
    diff = U - u

    # Apply homogeneous boundary conditions by eliminating the first and last entries
    diff = diff[1:-1]

    # Compute the L² error using the mass matrix
    error = np.sqrt(np.dot(diff, np.dot(M_global, diff)))

    return error

# Compute the L² error using the mass matrix
L2_error = mass_matrix_l2_error(U, u, M_global)

print(f"L² error: {L2_error}")

This concludes our tutorial on solving the 1D Poisson equation using the Finite Element Method. We have gone through the process of creating a mesh, defining the source term, constructing local stiffness and mass matrices, assembling global stiffness and mass matrices, assembling the global load vector, solving the linear system, plotting the solution, and calculating the L² error.

Remark

num_elements	L² error
8	9.79e-5
16	2.47e-5
32	6.19e-6
64	1.55e-6

위의 출력 결과에서 확인할 수 있듯이, 원하는 convergence order에 따라 num_elements 값을 증가시킬 때마다 L² error가 감소하고 있습니다.

a priori error estimate와 같은 convergence order를 얻기 위해서는, 원래 error를 h^p (h는 mesh size, p는 convergence order)에 비례하는 형태로 만들어야 합니다. 따라서, num_elements를 두 배로 증가시키면, mesh size인 h가 1/2로 줄어듭니다. 이에 따라 error는 h^p에 비례하므로, error는 h^p의 제곱근 형태로 감소합니다. 예를 들어, convergence order가 2일 경우, mesh size를 두 배로 줄이면 error는 1/4로 감소합니다.

따라서, 위의 출력 결과를 보면, num_elements를 두 배로 증가시킬 때마다 L² error가 약 4배씩 줄어들고 있습니다. 따라서, convergence order가 2임을 알 수 있습니다.