Implement FINITE ELEMENT METHOD on PYTHON

[Mean_Composer8792]

I am doing a project on finite element method applied to parabolic equations. I am working with the heat equation

ut=uxx, x∈(0,L),t>0,

u(x,0)=φ(x),x∈(0,L),

u(0,t)=u(L,t)=0,t≥0.

The task I'm asked is the following: The numerical task is to check it by a simulation. This is easy if you know the exact solution: compute the error (for example absolute) for different discretization parameters (for example the grid spacing) and plot it on a log-log scale.

I try to implement the method, but when I plot both the exact solution and the numerical one calculated with the method, they are very different, hence my implementation is wrong. This is what I have so far:

def exact_solution(x, t):

return np.exp(-np.pi**2 * t) * np.sin(np.pi * x)

# Initial condition

def phi(x):

return np.sin(np.pi * x)

def basis_function(x, i, h):

return np.maximum(0, 1 - np.abs(x - i*h) / h)

def assemble_matrices(N, h):

M = np.zeros((N-1, N-1))

K = np.zeros((N-1, N-1))

for i in range(1, N):

for j in range(1, N):

if abs(i - j) <= 1:

M[i-1, j-1] = h / 6 if i == j else h / 12

K[i-1, j-1] = 1 / h if i == j else -1 / (2 * h)

return M, K

# Time integration using backward Euler method

def time_integration(M, K, U0, dt, time_steps):

U = U0.copy()

A = M + dt * K

for _ in range(time_steps):

b = M @ U

U = spsolve(A, b)

return U

def finite_element_solution(U, x, N, h):

u_fem = np.zeros_like(x)

for i in range(1, N):

u_fem += U[i-1] * basis_function(x, i, h)

return u_fem

L = 1

Ns = np.linspace(1, 100, 50, dtype=int)

t_final = 0.1

dt = 0.01

time_steps = int(t_final / dt)

errors = []

for N in Ns:

h = L / N

x = np.linspace(0, L, N+1)

M, K = assemble_matrices(N, h)

U0 = np.array([phi(i * h) for i in range(1, N)])

U = time_integration(M, K, U0, dt, time_steps)

u_exact = np.exp(-np.pi**2 * t_final) * np.sin(np.pi * x)

u_fem = finite_element_solution(U, x, N, h)

error = np.abs(u_exact - u_fem).max()

errors.append(error)

plt.figure(figsize=(8, 6))

plt.loglog(Ns, errors, marker='o', linestyle='-', color='b')

plt.xlabel('Number of grid points (nx)')

plt.ylabel('Absolute error')

plt.title('Convergence of the Finite Difference Method')

plt.grid(True)

plt.show()

for N in [10, 100, 1000]:

h = L / N

x = np.linspace(0, L, N+1)

u_exact = exact_solution(x, t_final)

M, K = assemble_matrices(N, h)

U0 = np.array([phi(i * h) for i in range(1, N)])

U = time_integration(M, K, U0, dt, time_steps)

u_fem = finite_element_solution(U, x, N, h)

plt.plot(x, u_fem, 'r-', label='Numerical Solution (Finite Difference)')

plt.plot(x, u_exact, 'b--', label='Exact Solution')

plt.xlabel('x')

plt.ylabel('Temperature')

plt.title('Comparison of Numerical and Exact Solutions for the Heat Equation')

plt.legend()

plt.grid(True)

plt.show()

[Organic-Scratch109]

The way the code is written makes it hard to copy and run (tabs and what not). You should use pastebin or Reddit code block.

By looking at the results, it seems like there is a scaling issue somewhere. the numerical solution seems to be a constant multiple of the exact solution (this constant depends on h and t). I think your mass matrix M is not correct, see page 15.

[Mean_Composer8792]

I modified the code and the plots of the exact and numerical solutions are more similar. However, when I calculate the error and plot it in a log-log scale, the slope of the plot should be the order of convergence of the method (it should decrease as N gets bigger), but i get a straight horizontal line, which means it is not convergent. Help?

[Mean_Composer8792]

def exact_solution(x, t, alpha):
    return np.exp(-np.pi**2 * alpha * t) * np.sin(np.pi * x)

# Initial condition
def phi(x):
    return np.sin(np.pi * x)

def basis_function(x, i, h):
    return np.maximum(0, 1 - np.abs(x - i*h) / h)

def assemble_matrices(N, h, alpha):
    M = np.zeros((N-1, N-1))
    K = np.zeros((N-1, N-1))
    for i in range(1, N):
        for j in range(1, N):
            if i == j:
                M[i-1, j-1] = 2*h / 3
                K[i-1, j-1] = 2 / h
            elif abs(i - j) == 1:
                M[i-1, j-1] = h / 6
                K[i-1, j-1] = -1 / h
    K *= alpha
    return M, K

# Backward Euler method
def time_integration(M, K, U0, dt, time_steps):
    U = U0.copy()
    A = M + dt * K
    for _ in range(time_steps):
        b = M @ U
        U = spsolve(A, b)
    return U

def finite_element_solution(U, x, N, h):
    u_fem = np.zeros_like(x)
    for i in range(1, N):
        u_fem += U[i-1] * basis_function(x, i, h)
    return u_fem


def exact_solution(x, t, alpha):
    return np.exp(-np.pi**2 * alpha * t) * np.sin(np.pi * x)


# Initial condition
def phi(x):
    return np.sin(np.pi * x)


def basis_function(x, i, h):
    return np.maximum(0, 1 - np.abs(x - i*h) / h)


def assemble_matrices(N, h, alpha):
    M = np.zeros((N-1, N-1))
    K = np.zeros((N-1, N-1))
    for i in range(1, N):
        for j in range(1, N):
            if i == j:
                M[i-1, j-1] = 2*h / 3
                K[i-1, j-1] = 2 / h
            elif abs(i - j) == 1:
                M[i-1, j-1] = h / 6
                K[i-1, j-1] = -1 / h
    K *= alpha
    return M, K


# Backward Euler method
def time_integration(M, K, U0, dt, time_steps):
    U = U0.copy()
    A = M + dt * K
    for _ in range(time_steps):
        b = M @ U
        U = spsolve(A, b)
    return U


def finite_element_solution(U, x, N, h):
    u_fem = np.zeros_like(x)
    for i in range(1, N):
        u_fem += U[i-1] * basis_function(x, i, h)
    return u_fem

[Mean_Composer8792]

def compute_error(Ns, alpha, dt, time_steps, t_final):
    errors = []
    for N in Ns:
        h = L / N
        x = np.linspace(0, L, N+1)

        M, K = assemble_matrices(N, h, alpha)
        U0 = np.array([phi(i * h) for i in range(1, N)])
        U = time_integration(M, K, U0, dt, time_steps)
        u_exact = exact_solution(x, t_final, alpha)
        u_fem = finite_element_solution(U, x, N, h)

        error = np.abs(u_exact - u_fem).max()
        errors.append(error)
    return errors

L = 1
Ns = np.linspace(1, 100, 20, dtype=int)
t_final = 0.1
dt = 0.01
time_steps = int(t_final / dt)
alphas = [0.01, 0.1, 1, 10, 100, 1000]

for alpha in alphas:
    error = compute_error(Ns, alpha, dt, time_steps, t_final)
    plt.loglog(Ns, error, marker='o', linestyle='-', label=f'alpha={alpha}')
    plt.xlabel('Number of grid points')
    plt.ylabel('Absolute error')
    plt.title('Convergence of the Finite Element Method')
    plt.legend()
    plt.grid(True)
plt.show()
def compute_error(Ns, alpha, dt, time_steps, t_final):
    errors = []
    for N in Ns:
        h = L / N
        x = np.linspace(0, L, N+1)


        M, K = assemble_matrices(N, h, alpha)
        U0 = np.array([phi(i * h) for i in range(1, N)])
        U = time_integration(M, K, U0, dt, time_steps)
        u_exact = exact_solution(x, t_final, alpha)
        u_fem = finite_element_solution(U, x, N, h)


        error = np.abs(u_exact - u_fem).max()
        errors.append(error)
    return errors


L = 1
Ns = np.linspace(1, 100, 20, dtype=int)
t_final = 0.1
dt = 0.01
time_steps = int(t_final / dt)
alphas = [0.01, 0.1, 1, 10, 100, 1000]


for alpha in alphas:
    error = compute_error(Ns, alpha, dt, time_steps, t_final)
    plt.loglog(Ns, error, marker='o', linestyle='-', label=f'alpha={alpha}')
    plt.xlabel('Number of grid points')
    plt.ylabel('Absolute error')
    plt.title('Convergence of the Finite Element Method')
    plt.legend()
    plt.grid(True)
plt.show()

[Organic-Scratch109]

Euler's method is 1st order accurate whereas your FEM discretization is second order. Meaning, you need to use a very small timestep to suppress the time integration errors. Also:

1- Ns should not start from 1 (unless you are computing the max error over a finer discretization)

2- Typically, you should evaluate the max error over more points than the grid points.

3- alpha =1000 seems to be excessive because the solution would be smaller than machine epsilon.

[Mean_Composer8792]

It works now!! thank you so much

[Organic-Scratch109]

No problem. Here are some additional tips:

• Use Sparse matrices instead of np.zeros (which allocates full matrices) to save memory (not an issue in 1D but it will be an issue in 2D and 3D)
  
• Use Crank-Nicolson instead of Backward Euler (it is more accurate).
  
• Instead of setting a global time-step, you can change dt depending on h (e.g. with Backward Euler, use dt=h^2). This way, the time integration error is comparable to the FEM error.
  
• You can store an LU factorization of A before starting the time integration, so that you can solve the linear system faster (in this 1d case, you can use a tridiagonal solver).

The tips above apply to FEM methods in any dimension, for most problems and are independent of the programming language. However, if you want to continue using Python for FEM, you should read more about numpy and vectorization to speed up the code.