Matlab & Python (mit Numpy), das keine Kompilierungssprache ist, aber aufgrund der einfachen Handhabung von Arrays für numerische Berechnungen allgemein akzeptiert wird, ist langsam, wenn es mit einer for-Anweisung berechnet wird. Vermeiden Sie es daher so weit wie möglich! Ich denke, Sie hören oft den Diskurs. Was ist die Auswirkung? Ich habe versucht, es mit der expliziten Methode der zweidimensionalen Wellengleichung zu verifizieren.
Schauen wir uns die dominante Gleichung, die Wellengleichung und ihre numerische Lösung an.
Die zweidimensionale Wellengleichung
\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
Und die Zerstreuung des zentralen Systems davon ist
\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
Es wird sein. Wenn Sie dies so organisieren, dass es leicht zu berechnen ist,
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)
Es wird sein.
Wenn die Randbedingung eine Richtungsbedingung ist, ist dies einfach, da Sie nur den Randwert ersetzen müssen.
Matlab Schreiben wir zunächst mit einer geeigneten Schleife.
c = 1.0;
xmin = 0.;
ymin = 0.;
xmax = 1.;
ymax = 1.; %Der Berechnungsbereich ist[0,1],[0,1]
xmesh = 400;
ymesh = 400; %Die Anzahl der Unterteilungen beträgt x,400 für beide y-Achsen
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;%Die Anfangsgeschwindigkeit wird einem bestimmten Bereich gegeben.
x = xmin+dx/2:dx:xmax-dx/2;
y = ymin+dy/2:dy:ymax-dy/2;
t=0.;
tic %Ein praktischer Typ, der die Zeit mit Tic Toc messen kann
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;
%Geben Sie Diricre Bedingungen
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); %Zum Schluss zeichnen
In meiner Umgebung betrug die verstrichene Zeit 5,249230 Sekunden. Schreiben wir dies ohne Schleife mit Ausnahme des Zeitteils.
c = 1.0;
xmin = 0.;
ymin = 0.;
xmax = 1.;
ymax = 1.; %Der Berechnungsbereich ist[0,1],[0,1]
xmesh = 400;
ymesh = 400; %Die Anzahl der Unterteilungen beträgt x,400 für beide y-Achsen
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;%Die Anfangsgeschwindigkeit wird einem bestimmten Bereich gegeben.
x = xmin+dx/2:dx:xmax-dx/2;
y = ymin+dy/2:dy:ymax-dy/2;
t=0.;
tic %Ein praktischer Typ, der die Zeit mit Tic Toc messen kann
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;
%Geben Sie Diricre Bedingungen
u1(:,1)=0.;
u1(:,ymesh)=0.;
u1(1,:)=0.;
u1(xmesh,:)=0.;
end
toc
[X,Y] = meshgrid(x,y);
mesh(X,Y,u1)
Verwenden Sie einen kleinen Kopf, wenn Sie Code schreiben, der ordentlich aussieht. Dies liegt daran, dass bei Verwendung von diff die Indexabweichung berücksichtigt werden muss. Die Ausführungszeit betrug 2,734515 Sekunden. Der Unterschied in der Ausführungszeit ist nur etwa doppelt so groß, und der Eindruck ist, dass der Unterschied kleiner als erwartet ist.
Python Ich benutze oft Numpy, weil es bequem ist. Schreiben wir mit einer Schleife.
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. #Der Berechnungsbereich ist[0,1],[0,1]
xmesh = 400
ymesh = 400 #Die Anzahl der Unterteilungen beträgt x,400 für beide y-Achsen
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 #Die Anfangsgeschwindigkeit wird einem bestimmten Bereich gegeben.
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
#Geben Sie Diricre Bedingungen
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() #Erstellen eines Grundstücksbereichs
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, u1)
plt.show()
Die Berechnung wurde in 880.909 Sekunden abgeschlossen. Zu langsam.
Als nächstes kommt der Code ohne Schleifen.
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. #Der Berechnungsbereich ist[0,1],[0,1]
xmesh = 400
ymesh = 400 #Die Anzahl der Unterteilungen beträgt x,400 für beide y-Achsen
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 #Die Anfangsgeschwindigkeit wird einem bestimmten Bereich gegeben.
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
#Geben Sie Diricre Bedingungen
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() #Erstellen eines Grundstücksbereichs
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, u1)
plt.show()
Das Ergebnis ist 5,539646625518799 Sekunden. Es ist mehr als 100 Mal schneller als die Python-Loop-Version, aber es dauert viel länger als Matlab. Ist es schlecht geschrieben?
Vermeiden Sie es, Doppelschleifen in Python so weit wie möglich zu schwingen. Es dauert ungewöhnlich viel Rechenzeit. Mit Matlab war es etwas langsamer.
Nachtrag 2020/10/14: Eine Person, die mit Matlab vertraut ist, hat mich unterrichtet. Matlab scheint automatisch JIT-kompiliert zu werden. Wunderbar. https://jp.mathworks.com/products/matlab/matlab-execution-engine.html
Das Spiegelbild der Wellen ist wunderschön
c = 1.0;
xmin = 0.;
ymin = 0.;
xmax = 1.;
ymax = 1.; %Der Berechnungsbereich ist[0,1],[0,1]
xmesh = 100;
ymesh = 100; %Die Anzahl der Unterteilungen beträgt x,400 für beide y-Achsen
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);%Die Anfangsgeschwindigkeit wird einem bestimmten Bereich gegeben.
u1(60:80,40:60)=1e-7*sin(smallX).*sin(smallY);
mesh(X,Y,u1);
zlim([-1e-5 1e-5]);
tic %Ein praktischer Typ, der die Zeit mit Tic Toc messen kann
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)
Bei der Ausführung kann die Reflexion der Wellen wunderschön wie eine Animation gesehen werden.
Recommended Posts