Matlab & Python (avec numpy), qui n'est pas un langage de compilation mais qui est largement accepté pour les calculs numériques en raison de la facilité de manipulation des tableaux, est lent lorsqu'il est calculé avec une instruction for, alors évitez-le autant que possible! Je pense que vous entendez souvent le discours. Quel est l'impact? J'ai essayé de le vérifier en utilisant la méthode explicite de l'équation d'onde bidimensionnelle.
Jetons un coup d'œil à l'équation dominante, l'équation de la vague et sa solution numérique.
L'équation d'onde bidimensionnelle
\frac{\partial^2 u}{\partial t^2} - c^2\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)=0
Et la dispersion du système central de ceci est
\frac{u^{n+1}_{i,j}-2u^n_{i,j}+u^{n-1}_{i,j}}{\Delta t^2} - c^2\left(\frac{u^n_{i+1,j}-2u^n_{i,j}+u^n_{i-1,j}}{\Delta x^2}+\frac{u^n_{i,j+1}-2u^n_{i,j}+u^n_{i,j-1}}{\Delta y^2}\right)=0
Ce sera. Si vous organisez cela de manière à ce qu'il soit facile à calculer,
u^{n+1}_{i,j}=2u^n_{i,j}-u^{n-1}_{i,j} + c^2 \Delta t^2 \left(\frac{u^n_{i+1,j}-2u^n_{i,j}+u^n_{i-1,j}}{\Delta x^2}+\frac{u^n_{i,j+1}-2u^n_{i,j}+u^n_{i,j-1}}{\Delta y^2}\right)
Ce sera.
Si la condition aux limites est une condition Diricle, c'est facile car il vous suffit de substituer la valeur limite.
Matlab Tout d'abord, écrivons en utilisant une boucle selon le cas.
c = 1.0;
xmin = 0.;
ymin = 0.;
xmax = 1.;
ymax = 1.; %La zone de calcul est[0,1],[0,1]
xmesh = 400;
ymesh = 400; %Le nombre de divisions est x,400 pour les deux axes y
dx = (xmax-xmin)/xmesh;
dy = (ymax-ymin)/ymesh;
dt = 0.2*dx/c;
u0 = zeros(xmesh,ymesh); % u^{n-1}
u1 = zeros(xmesh,ymesh); % u^n
u2 = zeros(xmesh,ymesh); % u^{n+1}
u1(100:130,100:130)=1e-6;%La vitesse initiale est donnée à une certaine zone.
x = xmin+dx/2:dx:xmax-dx/2;
y = ymin+dy/2:dy:ymax-dy/2;
t=0.;
tic %Un gars pratique qui peut mesurer le temps avec tic toc
while t<1.0
for j = 2:ymesh-1
for i = 2:xmesh-1
u2(i,j) = 2*u1(i,j)-u0(i,j) + c*c*dt*dt*((u1(i+1,j)-2*u1(i,j)+u1(i-1,j))/(dx*dx) +(u1(i,j+1)-2*u1(i,j)+u1(i,j-1))/(dy*dy) );
end
end
u0=u1;
u1=u2;
t = t+dt;
%Donner des conditions Diricre
for i=1:xmesh
u1(i,1)=0.;
u1(i,ymesh)=0.;
end
for j=1:ymesh
u1(1,j)=0.;
u1(xmesh,j)=0.;
end
end
toc
[X,Y] = meshgrid(x,y);
mesh(X,Y,u1); %Enfin dessiner
Dans mon environnement, le temps écoulé était de 5 249 230 secondes. Écrivons ceci sans boucle, sauf pour la partie temporelle.
c = 1.0;
xmin = 0.;
ymin = 0.;
xmax = 1.;
ymax = 1.; %La zone de calcul est[0,1],[0,1]
xmesh = 400;
ymesh = 400; %Le nombre de divisions est x,400 pour les deux axes y
dx = (xmax-xmin)/xmesh;
dy = (ymax-ymin)/ymesh;
dt = 0.2*dx/c;
u0 = zeros(xmesh,ymesh); % u^{n-1}
u1 = zeros(xmesh,ymesh); % u^n
u2 = zeros(xmesh,ymesh); % u^{n+1}
u1(100:130,100:130)=1e-6;%La vitesse initiale est donnée à une certaine zone.
x = xmin+dx/2:dx:xmax-dx/2;
y = ymin+dy/2:dy:ymax-dy/2;
t=0.;
tic %Un gars pratique qui peut mesurer le temps avec tic toc
while t<1.0
u2(2:xmesh-1,2:ymesh-1) = 2*u1(2:xmesh-1,2:ymesh-1)-u0(2:xmesh-1,2:ymesh-1)+c*c*dt*dt*(diff(u1(:,2:ymesh-1),2,1)/(dx*dx)+diff(u1(2:xmesh-1,:),2,2)/(dy*dy));
u0=u1;
u1=u2;
t = t+dt;
%Donner des conditions Diricre
u1(:,1)=0.;
u1(:,ymesh)=0.;
u1(1,:)=0.;
u1(xmesh,:)=0.;
end
toc
[X,Y] = meshgrid(x,y);
mesh(X,Y,u1)
Utilisez une petite tête lorsque vous écrivez un code qui a l'air soigné. En effet, il est nécessaire de prendre en compte l'écart d'index lors de l'utilisation de diff. Le temps d'exécution était de 2,734515 secondes. La différence de temps d'exécution n'est qu'environ deux fois plus grande, et l'impression est que la différence est plus petite que prévu.
Python J'utilise souvent numpy parce que c'est pratique. Écrivons en utilisant une boucle.
import numpy as np
import matplotlib.pyplot as plt
from copy import deepcopy
import time
c = 1.0
xmin = 0.
ymin = 0.
xmax = 1.
ymax = 1. #La zone de calcul est[0,1],[0,1]
xmesh = 400
ymesh = 400 #Le nombre de divisions est x,400 pour les deux axes y
dx = (xmax-xmin)/xmesh
dy = (ymax-ymin)/ymesh
dt = 0.2*dx/c
u0 = np.zeros((xmesh, ymesh)) # u^{n-1}
u1 = np.zeros((xmesh, ymesh)) # u^n
u2 = np.zeros((xmesh, ymesh)) # u^{n+1}
u1[100:130, 100:130] = 1e-6 #La vitesse initiale est donnée à une certaine zone.
x = np.linspace(xmin+dx/2, xmax-dx/2, xmesh)
y = np.linspace(ymin+dy/2, ymax-dy/2, ymesh)
t = 0.
before = time.time()
while t < 1.:
for j in range(1, xmesh-1):
for i in range(1, ymesh-1):
u2[i, j] = 2*u1[i, j]-u0[i, j]+c*c*dt*dt * \
((u1[i+1, j]-2.*u1[i, j]+u1[i-1, j])/(dx*dx) +
(u1[i, j+1]-2.*u1[i, j]+u1[i, j-1])/(dy*dy))
u0 = deepcopy(u1)
u1 = deepcopy(u2)
t = t+dt
#Donner des conditions Diricre
for i in range(0,xmesh):
u1[i, 0] = 0.
u1[i, ymesh-1] = 0.
for j in range(0,ymesh):
u1[0, j] = 0.
u1[xmesh-1, j] = 0.
print(time.time()-before)
X, Y = np.meshgrid(x, y)
fig = plt.figure() #Créer une zone de tracé
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, u1)
plt.show()
Le calcul a été effectué en 880,909 secondes. Trop lent.
Vient ensuite le code sans boucles.
import numpy as np
import matplotlib.pyplot as plt
from copy import deepcopy
import time
c = 1.0
xmin = 0.
ymin = 0.
xmax = 1.
ymax = 1. #La zone de calcul est[0,1],[0,1]
xmesh = 400
ymesh = 400 #Le nombre de divisions est x,400 pour les deux axes y
dx = (xmax-xmin)/xmesh
dy = (ymax-ymin)/ymesh
dt = 0.2*dx/c
u0 = np.zeros((xmesh, ymesh)) # u^{n-1}
u1 = np.zeros((xmesh, ymesh)) # u^n
u2 = np.zeros((xmesh, ymesh)) # u^{n+1}
u1[100:130, 100:130] = 1e-6 #La vitesse initiale est donnée à une certaine zone.
x = np.linspace(xmin+dx/2, xmax-dx/2, xmesh)
y = np.linspace(ymin+dy/2, ymax-dy/2, ymesh)
t = 0.
before = time.time()
while t < 1.0:
u2[1:xmesh-1, 1:ymesh-1] = 2*u1[1:xmesh-1, 1:ymesh-1]-u0[1:xmesh-1, 1:ymesh-1]+c*c*dt*dt * \
(np.diff(u1[:, 1:ymesh-1], n=2, axis=0)/(dx*dx) +
np.diff(u1[1: xmesh-1, :], n=2, axis=1)/(dy*dy))
u0 = deepcopy(u1)
u1 = deepcopy(u2)
t = t+dt
#Donner des conditions Diricre
u1[:, 0] = 0.
u1[:, ymesh-1] = 0.
u1[0, :] = 0.
u1[xmesh-1, :] = 0.
print(time.time()-before)
X, Y = np.meshgrid(x, y)
fig = plt.figure() #Créer une zone de tracé
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, u1)
plt.show()
Le résultat est de 5,539646625518799 secondes. C'est plus de 100 fois plus rapide que la version de la boucle python, mais cela prend beaucoup plus de temps que Matlab. Est-ce mal écrit?
Évitez autant que possible de balancer les doubles boucles en Python. Cela prend un temps de calcul inhabituel. Avec Matlab, c'était un peu plus lent.
2020/10/14 Addendum: Une personne familière avec Matlab m'a appris. Matlab semble être automatiquement compilé JIT. Merveilleux. https://jp.mathworks.com/products/matlab/matlab-execution-engine.html
Le reflet des vagues est magnifique
c = 1.0;
xmin = 0.;
ymin = 0.;
xmax = 1.;
ymax = 1.; %La zone de calcul est[0,1],[0,1]
xmesh = 100;
ymesh = 100; %Le nombre de divisions est x,400 pour les deux axes y
dx = (xmax-xmin)/xmesh;
dy = (ymax-ymin)/ymesh;
dt = 0.1*dx/c;
u0 = zeros(xmesh,ymesh); % u^{n-1}
u1 = zeros(xmesh,ymesh); % u^n
u2 = zeros(xmesh,ymesh); % u^{n+1}
x = xmin+dx/2:dx:xmax-dx/2;
y = ymin+dy/2:dy:ymax-dy/2;
t=0.;
[X,Y] = meshgrid(x,y);
[smallX,smallY] = meshgrid(linspace(0,pi,21),linspace(0,pi,21));
u1(20:40,30:50)=2e-7*sin(smallX).*sin(smallY);%La vitesse initiale est donnée à une certaine zone.
u1(60:80,40:60)=1e-7*sin(smallX).*sin(smallY);
mesh(X,Y,u1);
zlim([-1e-5 1e-5]);
tic %Un gars pratique qui peut mesurer le temps avec tic toc
while t<10.0
u2(2:xmesh-1,2:ymesh-1) = 2*u1(2:xmesh-1,2:ymesh-1)-u0(2:xmesh-1,2:ymesh-1)+c*c*dt*dt*(diff(u1(:,2:ymesh-1),2,1)/(dx*dx)+diff(u1(2:xmesh-1,:),2,2)/(dy*dy));
u0=u1;
u1=u2;
t = t+dt;
u1(:,1)=0.;
u1(:,ymesh)=0.;
u1(1,:)=0.;
u1(xmesh,:)=0.;
mesh(X,Y,u1);
zlim([-1e-5 1e-5]);
pause(0.01);
end
toc
[X,Y] = meshgrid(x,y);
mesh(X,Y,u1)
Une fois exécuté, le reflet des vagues peut être vu magnifiquement comme une animation.
Recommended Posts