Один и тот же код интеграции в Matlab и python, Matlab стабилен, python взрывается

Вот алгоритм метода экспоненциальной разницы во времени с использованием исходного кода 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 в какой-то степени говорит, что какое-то значение переходит в NaN, но не для Matlab. Точно такой же код

Поэтому я хочу знать, почему и что происходит, вы можете запустить код самостоятельно.

Что заставляет меня чувствовать себя более странно, так это то, что первая итерация очень хорошо сочетается между ними!