0
$\begingroup$

I've written a code in Python to solve the 1D wave equation with the finite difference method (the explicit and the implicit methods). I'm trying to perform a mesh convergence study to estimate the empirical convergence rate of the method but I'm always getting an error of 2. I don't know what I'm doing wrong, whether my code is wrong or what. I'd appreciate the help.

# Équation des ondes
# Modélisation de la déformation u(x, t) d'une corde vibrante, de longueur L = 1, fixée à ses extremités
# Tester deux schémas aux différences finies pour approcher une solution de ce problème

from numpy import *
from numpy.linalg import *
from matplotlib.pyplot import *

def u0(x):
    return sin(pi*x) + 0.25*sin(10*pi*x)
def u1(x):
    return 0

print('Choix de la vitesse de transport c : ')
c = float(input('c = '))

def solex(x, t):
    return sin(pi*x) + 0.25*sin(10*pi*x)*cos(10*c*pi*t)

print('Choix de pas h : ')
h = float(input('h = '))

print('Choix du pas dt et du temps final T : ')
dt = float(input('dt = '))
T = float(input('T = '))

print("Schéma choisi pour le calcul : ")
print("1.- Schéma explicite")
print("2.- Schéma implicite")
meth = int(input('Choix = '))

# Maillage

N = int((1/h) - 1)
x = linspace(0, 1, N + 2)
M = int((T/dt) - 1)
t = linspace(0, T, M + 2)

# Matrice du système linéaire

u = ones(N - 1)*(-1)
v = ones(N - 1)*(-((c**2)*(dt**2)/(h**2)))
A = 2*eye(N)
B = (1 + 2*((c**2)*(dt**2)/(h**2)))*eye(N)
C = diag(u, 1)
D = diag(u, -1)
F = diag(v, 1)
G = diag(v, -1)

A = (1/(h**2))*(A + C + D)
B = B + F + G

# Itération

sol = zeros((N, M + 2))

for i in range(1, N + 1):
    sol[i - 1, 0] = u0(x[i])
for i in range(1, N + 1):
    sol[i - 1, 1] = sol[i - 1, 0] + dt*u1(x[i])

if (meth == 1):
    if (c*dt/h < 1):
        print("La condition CFL est satisfaite, donc le schéma converge")
    else:
        print("La condition CFL n'est pas satisfaite, donc le schéma ne converge pas")
    for j in range(1, size(t) - 1):
        sol[:, j + 1] = 2*sol[:, j] - sol[:, j - 1] - (c**2)*(dt**2)*dot(A, sol[:, j])
elif (meth == 2):
    for j in range(1, size(t) - 1):
        sol[:, j + 1] = dot(inv(B), 2*sol[:, j] - sol[:, j - 1])
                
sol = vstack([zeros(M + 2), sol, zeros(M + 2)])

# Calcul de l'erreur max

maxim1 = 0
erreur = 0

for j in range(0, size(t)):
        for i in range(0, size(x)):
            maxim1 = abs(solex(x[i], t[j]) - sol[i, j])
            if (maxim1 > erreur):
                erreur = maxim1
            else:
                erreur = erreur
        
print("L'erreur max est: ", erreur)

# Traçage de la solution approchée et exacte avec suffisamment de points (50-100 p.e)

plot(x, solex(x, T), label = 'Sol exacte')
plot(x, sol[:, size(t) - 1], label = 'Sol approché')
axis('equal')

title("Déformation d'une corde vibrante, de longueur L = 1, fixée à ses extremités")

# on ajoute les labels sur les axes
xlabel('x')
ylabel("Déformation à l'instant t")

legend()
# savefig('image.png')
show()

I've tried using dt = 0.002, and h = 0.1, 0.01, 0.001 to refine the mesh but I'm still getting an error of 2.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.