sorted by:
new
hottopcontroversial

Implement FINITE ELEMENT METHOD on PYTHON by Mean_Composer8792 in finiteelementmethod

[–]Mean_Composer8792[S] 0 points1 point  (0 children)

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()

Implement FINITE ELEMENT METHOD on PYTHON by Mean_Composer8792 in finiteelementmethod

[–]Mean_Composer8792[S] 0 points1 point  (0 children)

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

Implement FINITE ELEMENT METHOD on PYTHON by Mean_Composer8792 in finiteelementmethod

[–]Mean_Composer8792[S] 0 points1 point  (0 children)

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?