simualtion d'une onde qui se propage en 1D

[artchi01]

Bonjour, j'essaie de simuler la propagation d'une onde qui se propage en 1D dans un milieu limité à l'aide de la méthode des différences finies (Cf équation de d'Alembert, approximation de la dérivée seconde en espace et en temps). Cependant, toutes mes tentatives se sont avérées infructueuses. Je vous communique le code commenté. Merci par avance.
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 10 15:03:29 2024

@author: arthu
"""

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML

# Parameters
c = 1.0 # Wave speed
L = 10.0 # Length of the domain
T = 2.0 # Total time
dx = 0.01 # Spatial step size
dt = 0.005 # Temporal step size
x = np.arange(0, L, dx) # Spatial grid
t = np.arange(0, T, dt) # Temporal grid
Nx = len(x) # Number of spatial points
Nt = len(t) # Number of time points

U_0 = 5
w = 5

def matriceC(dt, dx, Nx, L):
C = np.zeros((Nx, Nx))
dx = L / (Nx - 1)

# Remplissage des valeurs internes de la matrice
for i in range(1, Nx - 1):
C[i, i - 1] = (dt * c / dx)**2
C[i, i] = 2 * (dt * c / dx)**2 + 1
C[i, i + 1] = (dt * c / dx)**2
# Conditions aux limites de Dirichlet
C[0, 0] = 1
C[0, 1] = -2
C[0, 2] = 1

C[Nx - 1, Nx - 3] = 1
C[Nx - 1, Nx - 2] = -2
C[Nx - 1, Nx - 1] = 1

return C

def matriceD(Nx, w, U_0, U, L):
dx = L / (Nx - 1)
matD = np.zeros(Nx)
# On applique les conditions aux limites trouvées plus haut
matD[0] = -U_0 * (w)**2
matD[-1] = 0
# Le corps de matD prend l'onde au temps précédent
matD[1:-1] = U[1:-1] - 2 * U[1:-1]
return matD

# Initial condition
U = np.zeros((Nt, Nx))
U[0] = np.sin(np.pi * x / L) # Condition initiale: sinusoïdale

C = matriceC(dt, dx, Nx, L)

# Time evolution
for i in range(1, Nt - 1):
D = matriceD(Nx, w, U_0, U[i - 1], L)
U[i] = np.linalg.solve(C, D) # Utilisation de np.linalg.solve au lieu de np.linalg.inv

# Animation function
fig, ax = plt.subplots()

def animate(i):
ax.clear()
ax.plot(x, U[i])
ax.set_ylim(-1, 1)
ax.set_title(f"Time: {i * dt:.2f}")
ax.set_xlabel('X')
ax.set_ylabel('Amplitude')

# Create the animation
ani = animation.FuncAnimation(fig, animate, frames=Nt, interval=50)

# Display the animation
plt.close()
HTML(ani.to_jshtml())
-----

[gts2]

Bonjour,

Avant de passer au codage, pourriez-vous indiquer la méthode mathématique utilisée (précisément pas uniquement différences finies), j'ai du mal à suivre.

[artchi01]

Bien sûr voici ci-dessus l’explication de ma méthode numérique

[gts2]

Avec votre méthode, pourquoi transformer en un système linéaire : on trouve donc l'état (n+1) est connu en fonction de l'état n et (n-1), soit en Python

C[i,j] = 2*C[i-1,j] - C[i-2,j]+k*(C[i-1,j+1]+C[i-1,j-1]-2*C[i-1,j])

et cela donne : (je n'ai pas gérer les CL)

Comme c'est une équation d'ordre 2, il faut des conditions initiales sur la fonction et sa dérivée
Les CL sont à définir à tout t, donc les C[0,i] pour tout i

[A voir en vidéo sur Futura]