Matlab和python同时集成代码，matlab稳定，python炸了
Same time integration code in Matlab and python, matlab is stable, python blows up
这里是一个指数时间差分法的算法，使用原来的
来自牛津的 Matlab 代码
clc
clear
% viscosity
nu = 0.5;
% Spatial grid and initial condition:
N = 128;
x = 2*pi*(0:N-1)'/N;
u = cos(x).*(1+sin(x));
v = fft(u);
% Precompute various ETDRK4 scalar quantities:
h = 0.01;
tot = 5;
% time step
k = [0:N/2-1 0 -N/2+1:-1]';
% wave numbers
L = 1/16*k.^2 - 1/16^3*nu*k.^4;
% Fourier multipliers
E = exp(h*L);
E2 = exp(h*L/2);
M =16;
% no. of points for complex means
r = exp(1i*pi*((1:M)-.5)/M); % roots of unity
% construct things on the
LR = h*L(:,ones(M,1)) + r(ones(N,1),:);
Q = h*real(mean( (exp(LR/2)-1)./LR ,2) );
f1 = h*real(mean( (-4-LR+exp(LR).*(4-3*LR+LR.^2))./LR.^3 ,2));
f2 = h*real(mean( (2+LR+exp(LR).*(-2+LR))./LR.^3 ,2));
f3 = h*real(mean( (-4-3*LR-LR.^2+exp(LR).*(4-LR))./LR.^3 ,2));
% Main time-stepping loop:
uu=u;tt=0;
vv=v;
NvNV = [];
aa = [];
bb = [];
cc = [];
NvNv = [];
NaNa = [];
NbNb = [];
NcNc = [];
nmax = floor(tot/h);
g = -0.5i*k;
%
for n = 1:nmax
t = n*h;
Nv = g.*fft(real(ifft(v)).^2);
a = E2.*v + Q.*Nv;
Na = g.*fft(real(ifft(a)).^2);
b = E2.*v + Q.*Na;
Nb = g.*fft(real(ifft(b)).^2);
c = E2.*a + Q.*(2*Nb-Nv);
Nc = g.*fft(real(ifft(c)).^2);
v = E.*v + Nv.*f1 + 2*(Na+Nb).*f2 + Nc.*f3;
u = real(ifft(v));
uu = [uu,u];
vv = [vv,v];
tt = [tt,t];
aa = [aa,a];
bb = [bb,b];
cc = [cc,c];
NvNv = [NvNv,Nv];
NaNa = [NaNa,Na];
NbNb = [NbNb,Nb];
NcNc = [NcNc,Nc];
end
Python代码
# spatial domain: 0,2pi
# time domain: 0,2
# init
import numpy as np
from matplotlib import pyplot as plt
np.seterr(all='raise')
plt.style.use('siads')
np.random.seed(0)
pi = np.pi
# viscosity
nu = 0.5
# mesh
mesh = 128
# time restriction
tot = 5
# distributed x point
x = np.linspace(0, 2.0 * pi, mesh, endpoint=False)
# force IC
u0 = np.cos(x)*(1.0 + np.sin(x))
N = mesh
k_array_noshift = np.fft.fftfreq(mesh)*mesh
u = u0
v = np.fft.fft(u)
h = 0.01
## follow them, set -N/2 to zer
k_array_noshift[mesh / 2] = 0
## Linear part in fourier space
L = 1.0/16.0*k_array_noshift**2 - nu*1.0/16.0**3 * k_array_noshift**4
## Fourier mulitplier
E = np.exp(h * L)
E2 = np.exp(h * L / 2.0)
## select number of points on the circle
M = 16
## choose radius 1, since k^2-k^4 ranges a lot... so 1 is enough to make sure only one singular point
r = np.exp(1j * np.pi * (np.arange(1, M + 1) - 0.5) / M)
r = r.reshape(1, -1)
r_on_circle = np.repeat(r, mesh, axis=0)
## define hL on M copies
LR = h * L
LR = LR.reshape(-1, 1)
LR = np.repeat(LR, M, axis=1)
## obtain hL on circle
LR = LR + r_on_circle
## important quantites used in ETDRK4
# g = -0.5*i*k
g = -0.5j * k_array_noshift
# averaged Q, f1,f2,f3
Q = h*np.real(np.mean( (np.exp(LR/2.0)-1)/LR, axis=1 ))
f1 = h*np.real(np.mean( (-4.0-LR + np.exp(LR)*(4.0-3.0*LR+LR**2))/LR**3, axis=1 ))
f2 = h*np.real(np.mean( (2.0+LR + np.exp(LR)*(-2.0 + LR))/LR**3, axis=1 ))
f3 = h*np.real(np.mean( (-4.0-3.0*LR - LR**2 + np.exp(LR)*(4.0 - LR))/LR**3, axis=1 ))
def compute_u2k_zeropad_dealiased(uk_):
u2k = np.fft.fft(np.real(np.fft.ifft(uk_))**2)
# three over two law
# N = uk.size
#
# # map uk to uk_fine
#
# uk_fine = np.hstack((uk[0:N / 2], np.zeros((N / 2)), uk[-N / 2:])) * 3.0 / 2.0
#
# # convert uk_fine to physical space
# u_fine = np.real(np.fft.ifft(uk_fine))
#
# # compute square
# u2_fine = np.square(u_fine)
#
# # compute fft on u2_fine
# u2k_fine = np.fft.fft(u2_fine)
#
# # convert u2k_fine to u2k
# u2k = np.hstack((u2k_fine[0:N / 2], u2k_fine[-N / 2:])) / 3.0 * 2.0
return u2k
print 'dt =', h
# np.linalg.norm(np.fft.ifft(uk_0)-u0) # very good
ntsnap = int(tot/h)
isnap = 0
tsnap = np.linspace(0, tot, ntsnap)
usnap = np.zeros((mesh, ntsnap))
uksnap = np.zeros((mesh, ntsnap),dtype=complex)
aasnap = np.zeros((mesh, ntsnap),dtype=complex)
bbsnap = np.zeros((mesh, ntsnap),dtype=complex)
ccsnap = np.zeros((mesh, ntsnap),dtype=complex)
Nvsnap = np.zeros((mesh, ntsnap),dtype=complex)
Nasnap = np.zeros((mesh, ntsnap),dtype=complex)
Nbsnap = np.zeros((mesh, ntsnap),dtype=complex)
Ncsnap = np.zeros((mesh, ntsnap),dtype=complex)
tcur = 0.0
## main loop time stepping
while tcur <= tot and isnap < ntsnap:
# print tcur
# record snap
# if abs(tcur - tsnap[isnap]) < 1e-2:
print ' current progress =', tcur / tsnap[-1] * 100, ' % '
u = np.real(np.fft.ifft(v))
usnap[:, isnap] = u
uksnap[:, isnap] = v
Nv = g * np.fft.fft(np.real(np.fft.ifft(v))**2)
a = E2 * v + Q * Nv
Na = g * np.fft.fft(np.real(np.fft.ifft(a))**2)
b = E2 * v + Q * Na
Nb = g * np.fft.fft(np.real(np.fft.ifft(b))**2)
c = E2 * a + Q * (2.0*Nb - Nv)
Nc = g * np.fft.fft(np.real(np.fft.ifft(c))**2)
v = E*v + Nv*f1 + 2.0*(Na + Nb)*f2 + Nc*f3
aasnap[:, isnap] = a
bbsnap[:, isnap] = b
ccsnap[:, isnap] = c
Nvsnap[:, isnap] = Nv
Nasnap[:, isnap] = Na
Nbsnap[:, isnap] = Nb
Ncsnap[:, isnap] = Nc
# algo: ETDRK4
# # 1. compute nonlinear part in first fractional step
# u2k = compute_u2k_zeropad_dealiased(uk)
# Nuk = g * u2k
#
# # 2. update fractional step
# uk_a = E2 * uk + Q * Nuk
#
# # 3. compute nonlinear part in second fractional step
# Nuk_a = g * compute_u2k_zeropad_dealiased(uk_a)
#
# # 4. update fractional step
# uk_b = E2 * uk + Q * Nuk_a
#
# # 5. compute nonlinear part in third fractional step
# Nuk_b = g * compute_u2k_zeropad_dealiased(uk_b)
#
# # 6. update fractional step
# uk_c = E2 * uk_a + Q * (2.0 * Nuk_b - Nuk)
#
# # 7 compute nonlinear part in the fourth fractional step
# Nuk_c = g * compute_u2k_zeropad_dealiased(uk_c)
#
# # final, update uk
# uk = E * uk + Nuk * f1 + 2.0 * (Nuk_a + Nuk_b) * f2 + Nuk_c * f3
# record time
tcur = tcur + h
isnap = isnap + 1
# save uksnap
np.savez('output_uk',uksnap=uksnap,f1=f1,f2=f2,f3=f3,Q=Q,LR=LR,L=L,g=g, Nasnap=Nasnap, Nvsnap=Nvsnap, Nbsnap=Nbsnap, Ncsnap=Ncsnap,
aasnap=aasnap, bbsnap=bbsnap, ccsnap=ccsnap, usnap=usnap)
# plot snapshots
plt.figure(figsize=(16,16))
for isnap in xrange(ntsnap):
plt.plot(x, usnap[:,isnap])
plt.xlabel('x')
plt.ylabel('u')
plt.savefig('snapshots.png')
# plot contours
from matplotlib import cm
fig = plt.figure(figsize=(8, 8))
X, Y = np.meshgrid(x, tsnap)
Z = usnap.transpose()
V = np.linspace(-10, 10, 500)
surf = plt.contourf(X, Y, Z, V, cmap=cm.seismic)
plt.savefig('contour.png')
Python 代码 运行s 在某种程度上表示某些值进入 NaN，但不适用于 matlab。相同的代码，正好
所以我想知道为什么以及发生了什么，你可以运行自己的代码。
更让我感到奇怪的是，第一次迭代，两者的匹配度非常高！
有关详细信息，请查看此 linkedIn link。
My LinkedIn Post that explains this problem
终于！！！！！！！！！一切正常!!!!!!!!!!!!!!!!!!!
我发现我可以让 python 和 matlab 一样好用（注意 python 总是在特定的物理时间爆炸到 inf），只需
将 numpy.fft 更改为 scipy.fft
请注意，我测试过，在我的代码中评估的 numpy.fft 和 scipy.fft 几乎相同，但在 1e-16..[=11= 附近有所不同]
但这对于这样一个混乱的系统模拟来说很重要。
如果我想让我的代码崩溃，只需将 scipy.fft 更改为 numpy.fft。
以我个人用户的观点，numpy.fft一定有什么玄机！由于 matlab fft 和 scipy fft 在我的代码中完全没有问题。
这里是一个指数时间差分法的算法，使用原来的 来自牛津的 Matlab 代码
clc
clear
% viscosity
nu = 0.5;
% Spatial grid and initial condition:
N = 128;
x = 2*pi*(0:N-1)'/N;
u = cos(x).*(1+sin(x));
v = fft(u);
% Precompute various ETDRK4 scalar quantities:
h = 0.01;
tot = 5;
% time step
k = [0:N/2-1 0 -N/2+1:-1]';
% wave numbers
L = 1/16*k.^2 - 1/16^3*nu*k.^4;
% Fourier multipliers
E = exp(h*L);
E2 = exp(h*L/2);
M =16;
% no. of points for complex means
r = exp(1i*pi*((1:M)-.5)/M); % roots of unity
% construct things on the
LR = h*L(:,ones(M,1)) + r(ones(N,1),:);
Q = h*real(mean( (exp(LR/2)-1)./LR ,2) );
f1 = h*real(mean( (-4-LR+exp(LR).*(4-3*LR+LR.^2))./LR.^3 ,2));
f2 = h*real(mean( (2+LR+exp(LR).*(-2+LR))./LR.^3 ,2));
f3 = h*real(mean( (-4-3*LR-LR.^2+exp(LR).*(4-LR))./LR.^3 ,2));
% Main time-stepping loop:
uu=u;tt=0;
vv=v;
NvNV = [];
aa = [];
bb = [];
cc = [];
NvNv = [];
NaNa = [];
NbNb = [];
NcNc = [];
nmax = floor(tot/h);
g = -0.5i*k;
%
for n = 1:nmax
t = n*h;
Nv = g.*fft(real(ifft(v)).^2);
a = E2.*v + Q.*Nv;
Na = g.*fft(real(ifft(a)).^2);
b = E2.*v + Q.*Na;
Nb = g.*fft(real(ifft(b)).^2);
c = E2.*a + Q.*(2*Nb-Nv);
Nc = g.*fft(real(ifft(c)).^2);
v = E.*v + Nv.*f1 + 2*(Na+Nb).*f2 + Nc.*f3;
u = real(ifft(v));
uu = [uu,u];
vv = [vv,v];
tt = [tt,t];
aa = [aa,a];
bb = [bb,b];
cc = [cc,c];
NvNv = [NvNv,Nv];
NaNa = [NaNa,Na];
NbNb = [NbNb,Nb];
NcNc = [NcNc,Nc];
end
Python代码
# spatial domain: 0,2pi
# time domain: 0,2
# init
import numpy as np
from matplotlib import pyplot as plt
np.seterr(all='raise')
plt.style.use('siads')
np.random.seed(0)
pi = np.pi
# viscosity
nu = 0.5
# mesh
mesh = 128
# time restriction
tot = 5
# distributed x point
x = np.linspace(0, 2.0 * pi, mesh, endpoint=False)
# force IC
u0 = np.cos(x)*(1.0 + np.sin(x))
N = mesh
k_array_noshift = np.fft.fftfreq(mesh)*mesh
u = u0
v = np.fft.fft(u)
h = 0.01
## follow them, set -N/2 to zer
k_array_noshift[mesh / 2] = 0
## Linear part in fourier space
L = 1.0/16.0*k_array_noshift**2 - nu*1.0/16.0**3 * k_array_noshift**4
## Fourier mulitplier
E = np.exp(h * L)
E2 = np.exp(h * L / 2.0)
## select number of points on the circle
M = 16
## choose radius 1, since k^2-k^4 ranges a lot... so 1 is enough to make sure only one singular point
r = np.exp(1j * np.pi * (np.arange(1, M + 1) - 0.5) / M)
r = r.reshape(1, -1)
r_on_circle = np.repeat(r, mesh, axis=0)
## define hL on M copies
LR = h * L
LR = LR.reshape(-1, 1)
LR = np.repeat(LR, M, axis=1)
## obtain hL on circle
LR = LR + r_on_circle
## important quantites used in ETDRK4
# g = -0.5*i*k
g = -0.5j * k_array_noshift
# averaged Q, f1,f2,f3
Q = h*np.real(np.mean( (np.exp(LR/2.0)-1)/LR, axis=1 ))
f1 = h*np.real(np.mean( (-4.0-LR + np.exp(LR)*(4.0-3.0*LR+LR**2))/LR**3, axis=1 ))
f2 = h*np.real(np.mean( (2.0+LR + np.exp(LR)*(-2.0 + LR))/LR**3, axis=1 ))
f3 = h*np.real(np.mean( (-4.0-3.0*LR - LR**2 + np.exp(LR)*(4.0 - LR))/LR**3, axis=1 ))
def compute_u2k_zeropad_dealiased(uk_):
u2k = np.fft.fft(np.real(np.fft.ifft(uk_))**2)
# three over two law
# N = uk.size
#
# # map uk to uk_fine
#
# uk_fine = np.hstack((uk[0:N / 2], np.zeros((N / 2)), uk[-N / 2:])) * 3.0 / 2.0
#
# # convert uk_fine to physical space
# u_fine = np.real(np.fft.ifft(uk_fine))
#
# # compute square
# u2_fine = np.square(u_fine)
#
# # compute fft on u2_fine
# u2k_fine = np.fft.fft(u2_fine)
#
# # convert u2k_fine to u2k
# u2k = np.hstack((u2k_fine[0:N / 2], u2k_fine[-N / 2:])) / 3.0 * 2.0
return u2k
print 'dt =', h
# np.linalg.norm(np.fft.ifft(uk_0)-u0) # very good
ntsnap = int(tot/h)
isnap = 0
tsnap = np.linspace(0, tot, ntsnap)
usnap = np.zeros((mesh, ntsnap))
uksnap = np.zeros((mesh, ntsnap),dtype=complex)
aasnap = np.zeros((mesh, ntsnap),dtype=complex)
bbsnap = np.zeros((mesh, ntsnap),dtype=complex)
ccsnap = np.zeros((mesh, ntsnap),dtype=complex)
Nvsnap = np.zeros((mesh, ntsnap),dtype=complex)
Nasnap = np.zeros((mesh, ntsnap),dtype=complex)
Nbsnap = np.zeros((mesh, ntsnap),dtype=complex)
Ncsnap = np.zeros((mesh, ntsnap),dtype=complex)
tcur = 0.0
## main loop time stepping
while tcur <= tot and isnap < ntsnap:
# print tcur
# record snap
# if abs(tcur - tsnap[isnap]) < 1e-2:
print ' current progress =', tcur / tsnap[-1] * 100, ' % '
u = np.real(np.fft.ifft(v))
usnap[:, isnap] = u
uksnap[:, isnap] = v
Nv = g * np.fft.fft(np.real(np.fft.ifft(v))**2)
a = E2 * v + Q * Nv
Na = g * np.fft.fft(np.real(np.fft.ifft(a))**2)
b = E2 * v + Q * Na
Nb = g * np.fft.fft(np.real(np.fft.ifft(b))**2)
c = E2 * a + Q * (2.0*Nb - Nv)
Nc = g * np.fft.fft(np.real(np.fft.ifft(c))**2)
v = E*v + Nv*f1 + 2.0*(Na + Nb)*f2 + Nc*f3
aasnap[:, isnap] = a
bbsnap[:, isnap] = b
ccsnap[:, isnap] = c
Nvsnap[:, isnap] = Nv
Nasnap[:, isnap] = Na
Nbsnap[:, isnap] = Nb
Ncsnap[:, isnap] = Nc
# algo: ETDRK4
# # 1. compute nonlinear part in first fractional step
# u2k = compute_u2k_zeropad_dealiased(uk)
# Nuk = g * u2k
#
# # 2. update fractional step
# uk_a = E2 * uk + Q * Nuk
#
# # 3. compute nonlinear part in second fractional step
# Nuk_a = g * compute_u2k_zeropad_dealiased(uk_a)
#
# # 4. update fractional step
# uk_b = E2 * uk + Q * Nuk_a
#
# # 5. compute nonlinear part in third fractional step
# Nuk_b = g * compute_u2k_zeropad_dealiased(uk_b)
#
# # 6. update fractional step
# uk_c = E2 * uk_a + Q * (2.0 * Nuk_b - Nuk)
#
# # 7 compute nonlinear part in the fourth fractional step
# Nuk_c = g * compute_u2k_zeropad_dealiased(uk_c)
#
# # final, update uk
# uk = E * uk + Nuk * f1 + 2.0 * (Nuk_a + Nuk_b) * f2 + Nuk_c * f3
# record time
tcur = tcur + h
isnap = isnap + 1
# save uksnap
np.savez('output_uk',uksnap=uksnap,f1=f1,f2=f2,f3=f3,Q=Q,LR=LR,L=L,g=g, Nasnap=Nasnap, Nvsnap=Nvsnap, Nbsnap=Nbsnap, Ncsnap=Ncsnap,
aasnap=aasnap, bbsnap=bbsnap, ccsnap=ccsnap, usnap=usnap)
# plot snapshots
plt.figure(figsize=(16,16))
for isnap in xrange(ntsnap):
plt.plot(x, usnap[:,isnap])
plt.xlabel('x')
plt.ylabel('u')
plt.savefig('snapshots.png')
# plot contours
from matplotlib import cm
fig = plt.figure(figsize=(8, 8))
X, Y = np.meshgrid(x, tsnap)
Z = usnap.transpose()
V = np.linspace(-10, 10, 500)
surf = plt.contourf(X, Y, Z, V, cmap=cm.seismic)
plt.savefig('contour.png')
Python 代码 运行s 在某种程度上表示某些值进入 NaN，但不适用于 matlab。相同的代码，正好
所以我想知道为什么以及发生了什么，你可以运行自己的代码。
更让我感到奇怪的是，第一次迭代，两者的匹配度非常高！
有关详细信息，请查看此 linkedIn link。
My LinkedIn Post that explains this problem
终于！！！！！！！！！一切正常!!!!!!!!!!!!!!!!!!!
我发现我可以让 python 和 matlab 一样好用（注意 python 总是在特定的物理时间爆炸到 inf），只需
将 numpy.fft 更改为 scipy.fft
请注意，我测试过，在我的代码中评估的 numpy.fft 和 scipy.fft 几乎相同，但在 1e-16..[=11= 附近有所不同]
但这对于这样一个混乱的系统模拟来说很重要。
如果我想让我的代码崩溃，只需将 scipy.fft 更改为 numpy.fft。
以我个人用户的观点，numpy.fft一定有什么玄机！由于 matlab fft 和 scipy fft 在我的代码中完全没有问题。