Choses à faire

Je voulais simuler une équation différentielle incluant le laplacien avec python, mais l'autre jour, j'ai proposé une méthode pour traiter la différence du laplacien autre que l'instruction for (j'ai réinventé la roue), j'ai donc écrit une animation de l'équation de diffusion. C'était.

Quand je l'ai cherché à nouveau après l'avoir trouvé, il y avait une personne qui était compatible avec ce que j'ai fait il y a quelques années. L'article a été très utile: Accélérer les calculs numériques à l'aide de NumPy / SciPy: Application 1

L'équation différentielle à résoudre est une équation de diffusion unidimensionnelle exprimée par l'équation suivante: $ \frac{\partial }{\partial t}\psi(x,t)= D\frac{\partial^2}{\partial x^2}\psi(x,t) $ La condition initiale était donnée par Gaussian, et le développement temporel était la méthode d'Euler pour le moment.

Code que j'ai écrit

`Diffusion_1D_ani.ipynb`


%matplotlib nbagg
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import ArtistAnimation

# calculate Laplacian with periodic boundary condition
def Lap(L, M):
    # expand vector
    VL = [L[0]]
    VR = [L[-1]]
    Le = np.hstack([VR ,L, VL])
    Lc = np.dot(M, Le) # central diff
    return Lc[1:N+1]

# generate 1-D lattice and initialize
N = 64
x = np.linspace(-10, 10, N)
psi = np.exp(-x**2/7)

# time and diffusion const
dt = 0.05
T = np.arange(0.0, 37, dt)
D = 3

# make (N+2)x(N+2) matrix to calc Laplacian
I = np.eye(N+2,N+2)
S = np.vstack((I[1:],np.array([0]*(N+2))))
M = S + S.T -2*I 

# make animation
fig = plt.figure()
ims = []
for i in range(len(T)):
    psi += dt*(D*Lap(psi, M))
    line = plt.plot(x, psi, 'b')
    ims.append(line)

ani = ArtistAnimation(fig, ims, interval=5, blit=True, repeat=True)

plt.show()

Exemple d'exécution

Ingéniosité

Lors du calcul du laplacien, les valeurs aux extrémités droite et gauche du tableau dimensionnel $ N $ original $ \ psi $ étaient attachées aux extrémités gauche et droite avec np.hstack pour créer un tableau dimensionnel $ N + 2 $. En conséquence, le terme de diffusion incluant la condition aux limites périodique pourrait être calculé automatiquement en prenant le produit interne avec la $ (N + 2) \ times (N + 2) $ matrice $ M $ préparée à l'avance. La pièce ajoutée pour la condition aux limites n'est pas nécessaire, elle est donc ignorée à la fin de la fonction.

Résumé

Problème: je voulais le rendre visible pour le moment, donc je n'ai pas pris en compte la stabilité numérique et l'erreur de 1 mm. Je pense qu'il y a quelque chose dans NumPy qui peut effectuer l'intégration numérique rapidement et en toute sécurité, alors je vais enquêter et faire quelque chose à ce sujet. J'ajouterai à ces questions lorsque des progrès seront réalisés à l'avenir.

Futur: Je pense que j'aimerais essayer une simulation du système de diffusion de réaction car il semble que je puisse faire une version 2D en étendant cette version 1D de manière appropriée. (Si ce n'est pas une équation qui inclut un terme non linéaire très laid, ce sera en quelque sorte ...)

Post-scriptum:

(11/23) Si nous faisons un compromis sur la différence à terme, nous pourrions mettre un terme non linéaire en 3 lignes. Il semble qu'une équation de Navier-Stokes unidimensionnelle puisse être créée.

def nl(L, N):
    Lf = np.hstack((L,L[0]))
    return L * np.diff(Lf)

Dans le cas de la différence directe, le calcul a divergé immédiatement si non seulement le nombre de Reynolds, mais aussi les conditions initiales et la taille du pas de temps n'ont pas été soigneusement sélectionnés. Que devrais-je faire

Animation de l'équation de diffusion avec NumPy