Considérant le Schéma Leapfrog suivant utilisé pour discrétiser un équation d'onde vectorielle avec des conditions initiales et des conditions aux limites périodiques données. J'ai implémenté le schéma et maintenant je veux faire des tests de convergence numérique pour montrer que le schéma est du second ordre dans l'espace et dans le temps.

Je suis principalement aux prises avec deux points ici:

  1. Je ne suis pas sûr à 100% d'avoir correctement mis en œuvre le schéma. Je voulais vraiment utiliser le slicing car c'est tellement plus rapide que d'utiliser des boucles.
  2. Je ne sais pas vraiment comment obtenir le bon tracé d'erreur, car je ne sais pas quelle norme utiliser. Dans les exemples que j'ai trouvés (ils étaient en 1D), nous avons toujours utilisé la norme L2.
import numpy as np
import matplotlib.pyplot as plt


# Initial conditions
def p0(x):
    return np.cos(2 * np.pi * x)


def u0(x):
    return -np.cos(2 * np.pi * x)


# exact solution
def p_exact(x, t):
    # return np.cos(2 * np.pi * (x + t))
    return p0(x + t)


def u_exact(x, t):
    # return -np.cos(2 * np.pi * (x + t))
    return u0(x + t)


# function for doing one time step, considering the periodic boundary conditions
def leapfrog_step(p, u):
    p[1:] += CFL * (u[:-1] - u[1:])
    p[0] = p[-1]
    u[:-1] += CFL * (p[:-1] - p[1:])
    u[-1] = u[0]
    return p, u


# Parameters
CFL = 0.3
LX = 1  # space length
NX = 100  # number of space steps

T = 2  # end time

NN = np.array(range(50, 1000, 50))  # list of discretizations
Ep = []
Eu = []
for NX in NN:
    print(NX)
    errorsp = []
    errorsu = []
    x = np.linspace(0, LX, NX)    # space grid
    dx = x[1] - x[0]  # spatial step
    dt = CFL * dx  # time step
    t = np.arange(0, T, dt)  # time grid

    # TEST

    # time loop
    for time in t:
        if time == 0:
            p = p0(x)
            u = u0(x)
        else:
            p, u = leapfrog_step(p, u)
            errorsp.append(np.linalg.norm((p - p_exact(x, time)), 2))
            errorsu.append(np.linalg.norm((u - u_exact(x, time)), 2))
    errorsp = np.array(errorsp) * dx ** (1 / 2)
    errorsu = np.array(errorsu) * dx ** (1 / 2)
    Ep.append(errorsp[-1])
    Eu.append(errorsu[-1])

# plot the error
plt.figure(figsize=(8, 5))
plt.xlabel("$Nx$")
plt.ylabel(r'$\Vert p-\bar{p}\Vert_{L_2}$')
plt.loglog(NN, 15 / NN ** 2, "green", label=r'$O(\Delta x^{2})$')
plt.loglog(NN, Ep, "o", label=r'$E_p$')
plt.loglog(NN, Eu, "o", label=r'$E_u$')
plt.legend()
plt.show()


J'apprécierais vraiment que quelqu'un puisse vérifier rapidement la mise en œuvre du schéma et une indication sur la façon d'obtenir le tracé d'erreur.

0
jerremaier 30 oct. 2020 à 21:26

1 réponse

Meilleure réponse

En dehors de l'initialisation, je ne vois aucune erreur dans votre code.

En ce qui concerne l'initialisation, considérons la première étape. Là, vous devez calculer, selon la description de la méthode, des approximations pour p(dt,j*dx) à partir des valeurs de p(0,j*dx) et u(0.5*dt, (j+0.5)*dx). Cela signifie que vous devez initialiser à time==0

u = u_exact(x+0.5*dx, 0.5*dt).

Et doivent également comparer la solution alors obtenue avec u_exact(x+0.5*dx, time+0.5*dt).

Le fait que vous ayez obtenu le bon ordre est pour l'OMI plus un artefact du problème de test qu'un algorithme accidentellement toujours correct.

Si aucune solution exacte n'est connue, ou si vous souhaitez utiliser un algorithme plus réaliste dans le test, vous devrez calculer les valeurs initiales u à partir de p(0,x) et u(0,x) via des développements de Taylor

u(t,x) = u(0,x) + t*u_t(0,x) + 0.5*t^2*u_tt(0,x) + ...

u(0.5*dt,x) = u(0,x) - 0.5*dt*p_x(0,x) + 0.125*dt^2*u_xx(0,x) + ...
            = u(0,x) - 0.5*CFL*(p(0,x+0.5*dx)-p(0,x-0.5*dx)) 
                     + 0.125*CFL^2*(u(0,x+dx)-2*u(0,x)+u(0,x-dx)) + ...

Il pourrait être suffisant de prendre juste l'expansion linéaire,

u[j] = u0(x[j]+0.5*dx) - 0.5*CFL*(p0(x[j]+dx)-p0(x[j]) 

Ou avec des opérations de tableau

p = p0(x)
u = u0(x+0.5*dx)
u[:-1] -= 0.5*CFL*(p[1:]-p[:-1])
u[-1]=u[0] 

Car alors l'erreur de second ordre dans les données initiales s'ajoute simplement à l'erreur générale de second ordre.


Vous voudrez peut-être changer la grille d'espace en x = np.linspace(0, LX, NX+1) pour avoir dx = LX/NX.


Je définirais la solution exacte et la condition initiale dans l'autre sens, car cela permet plus de flexibilité dans les problèmes de test.

# components of the solution
def f(x): return np.cos(2 * np.pi * x)
def g(x): return 2*np.sin(6 * np.pi * x)
# exact solution
def u_exact(x,t): return  f(x+t)+g(x-t)
def p_exact(x,t): return -f(x+t)+g(x-t)
# Initial conditions
def u0(x): return u_exact(x,0)
def p0(x): return p_exact(x,0)
0
Lutz Lehmann 5 nov. 2020 à 13:11