黒木玄
2018-01-23~2018-09-08
このノートブックは
の続き. このノートブックでは in-place でFFTを使うことによってメモリ効率を最適化することを行いたい. 例としてKdV方程式を扱う.
using Base64
using FFTW
using LinearAlgebra
const e = ℯ
linspace(start, stop, length) = range(start, stop=stop, length=length)
using PyPlot
default_figsize = (6.4, 4.8) # このデフォルトサイズは大き過ぎるので
rc("figure", figsize=(3,2.4)) # このように小さくしておく.
using PyCall
const animation = pyimport("matplotlib.animation")
function displayfile(mimetype, file)
open(file) do f
base64text = base64encode(f)
display("text/html", """<img src="data:$mimetype;base64,$base64text">""")
end
end
displayfile (generic function with 1 method)
@show Sys.CPU_THREADS
FFTW.set_num_threads(Sys.CPU_THREADS)
Sys.CPU_THREADS = 8
周期 $L=2\pi K$ を持つ函数を
$$ f(x) = \sum_{k=-N/2+1}^{N/2} (-1)^{k-1} a_k\;e^{2\pi i(k-1)x/L} = \sum_{k=-N/2+1}^{N/2} (-1)^{k-1} a_k\;e^{i(k-1)x/K} $$の形の有限Fourier級数でよく近似できていると仮定する. これの導函数は
$$ f'(x) = \sum_{k=-N/2+1}^{N/2} (-1)^{k-1} \frac{i(k-1)}{K}a_k\;e^{i(k-1)x/K} $$になる. すなわち, $f(x)$ を微分する操作は, $-N/2< k \leqq N/2$ のとき $a_k$ に $i(k-1)/K$ をかける操作に対応している. $-N/2< k \leqq N/2$ のとき $a_k$ に $i(k-1)/K$ をかける操作は(丸め誤差)を除いて, 有限Fourier級数の正確な微分を与えることに注意せよ. この計算方法は小さな $h$ について $(f(x+h/2)-f(x-h/2))/h$ のような方法で微分を求める方法よりも誤差が小さい.
さらに $x$ を
$$ \Delta x = \frac{L}{N} = \frac{2\pi K}{N}, \quad x_j = (j-1)\Delta x -\frac{L}{2} = L\left(\frac{j-1}{N} - \frac{1}{2}\right) $$で離散化したとする. このとき,
$$ f(x_j) = \sum_{k=-N/2-1}^{N/2} a_k\;e^{i(k-1)(j-1)/N}. $$これは $a_{k+N}=a_k$ という仮定のもとで
$$ f(x_j) = \sum_{k=1}^N a_k\; e^{i(k-1)(j-1)/N} $$と書き直される. これはJulia言語の離散Fourier変換のスタイルに一致している. この表示のもとで, $f(x)$ を微分する操作は, $1\leqq k\leqq N/2$ のとき $i(k-1)/N$ を $a_k$ にかけ, $N/2+1\leqq k\leqq N$ のとき $i(k-1-N)/N$ を $a_k$ にかける操作に対応している. FFT_Data
における D
を成分ごとにかける操作はまさにこの操作になっている.
?ifft
search: ifft ifft! ifftshift irfft plan_ifft plan_ifft! Cptrdiff_t plan_irfft
ifft(A [, dims])
Multidimensional inverse FFT.
A one-dimensional inverse FFT computes
$$ \operatorname{IDFT}(A)[k] = \frac{1}{\operatorname{length}(A)} \sum_{n=1}^{\operatorname{length}(A)} \exp\left(+i\frac{2\pi (n-1)(k-1)} {\operatorname{length}(A)} \right) A[n]. $$A multidimensional inverse FFT simply performs this operation along each transformed dimension of A
.
# right-open interval [a,b)
#
roi(a, b, N) = collect(linspace(a, b-(b-a)/N, N))
# x axis (Assume the priodic boundary condition with period L)
#
x_axis(L, N) = roi(-L/2, L/2, N)
# k axis (k = the wave number)
#
function k_axis(N)
@assert iseven(N)
Ndiv2 = div(N,2)
vcat(roi(0,Ndiv2,Ndiv2), roi(-Ndiv2,0,Ndiv2))
end
# FFT Data
#
mutable struct FFT_Data{T<:Real, S}
K::T
L::T
N::Int64
x::Array{T,1}
k::Array{T,1}
D::Array{Complex{T},1}
D2::Array{Complex{T},1}
D3::Array{Complex{T},1}
D4::Array{Complex{T},1}
FFT::S
end
# FFT \ (D.^m .* (FFT * u)) = (d/dx)^m u
#
function FFT_Data(K, N)
L = 2π*K
x = x_axis(L, N)
k = k_axis(N)
D = im .* k ./ K
FFT = plan_fft(Array{Complex{Float64},1}(undef, N))
FFT_Data(Float64(K), L, N, x, k, D, D.^2, D.^3, D.^4, FFT)
end
FFT_Data
# nearly zero
#
function imagisnealyzero(z)
maximum(abs, imag(z))/(maximum(abs, real(z)) + 1e-15) < 1e-3
end
# plot complex number array u on x
#
function plot_u(x, u)
plot(x, real(u), label="Re", lw=0.8)
imagisnealyzero(u) || plot(x, imag(u), label="Im", lw=0.8, ls="--")
grid(ls=":")
xticks(fontsize=8)
yticks(fontsize=8)
end
plot_u (generic function with 1 method)
# Application - the solution of the heat equation u_t = u_{xx}
#
# u(t) = exp(t (∂/∂x)^2) u(0)
#
sol_of_heat_eq(o::FFT_Data, u0, t) = o.FFT \ (exp.(t .* o.D2) .* (o.FFT * u0))
N = 2^8
K = 3
o = FFT_Data(K, N);
f0(x) = exp(-4(x+4)^2) + exp(-2(x-3)^2)
u0 = f0.(o.x)
t = roi(0,2,10)
@time u = [sol_of_heat_eq(o, u0, t) for t in t]
figure(figsize=(10, 3))
for j in 1:10
subplot(2,5,j)
plot_u(o.x, u[j])
ylim(-0.05,1.05)
end
tight_layout()
0.937872 seconds (2.18 M allocations: 113.085 MiB, 4.17% gc time)
function plot_u(o, u, nrows, j, ymin, ymax)
subplot(nrows, 5, j+1)
plot_u(o.x, u)
ylim(ymin, ymax)
xticks(fontsize=8)
yticks(fontsize=8)
title("j = $j", fontsize=8)
end
function plot_KdV(o::FFT_Data, u, t; thin=10, ymin=-7.0, ymax=15.0)
nplots = div(length(t), thin)
nrows = div(nplots+4, 5)
figure(figsize=(10, 1.75nrows))
j = 0
plot_u(o, @view(u[:,j+1]), nrows, j, ymin, ymax)
for j in 1:nplots-1
plot_u(o, @view(u[:,thin*j+1]), nrows, j, ymin, ymax)
end
tight_layout()
end
function anim_KdV(gifname, o::FFT_Data, u, t; thin=2, interval=100, ymin=-7.0, ymax=15.0)
function plotframe(j)
clf()
plot(o.x, real(u[:, (j-1)*thin + 1]), lw=1.0)
grid(ls=":")
ylim(ymin, ymax)
xticks(fontsize=8)
yticks(fontsize=8)
title("$(gifname): j = $(j-1)", fontsize=10)
plot()
end
fig = figure(figsize=(4, 3))
n = (length(t)-1)÷thin + 1
frames = [1;1;1;1;1;1;1;1;1; 1:n; n;n;n;n;n;n;n;n;n;]
ani=animation[:FuncAnimation](fig, plotframe, frames=frames, interval=interval)
ani[:save]("$(gifname).gif", writer="imagemagick")
clf()
displayfile("image/gif", "$(gifname).gif")
end
anim_KdV (generic function with 1 method)
FFT * v
, FFT \ v
の形式で利用した場合¶# FFT*v and FFT\v version
function update_KdV(o::FFT_Data, u, Δt, niters)
v = o.FFT * u
for iter in 1:niters
v .= exp.(-Δt .* o.D3) .* v
v .-= 3.0 .* Δt .* o.D .*(o.FFT * (o.FFT \ v).^2)
end
o.FFT \ v
end
function solve_KdV(o::FFT_Data{T}, f0, tmax) where T
Δt = 1/o.N^2
skip = floor(Int, 0.005/Δt)
t = 0:skip*Δt:tmax
M = length(t)
u = Array{Complex{T},2}(undef, o.N, M)
u[:,1] = Complex.(f0.(o.x))
for j in 2:M
u[:,j] = update_KdV(o, @view(u[:,j-1]), Δt, skip)
end
return real(u), t
end
N = 2^8
K = 20/(2π)
o = FFT_Data(K, N);
f0(x) = -5*sin(π*x/10)
Δt = 1/N^2
skip = 10*floor(Int, 0.005/Δt)
tmax = 1.0
@time u, t = solve_KdV(o, f0, tmax)
sleep(0.1)
plot_KdV(o, u, t)
8.667141 seconds (2.12 M allocations: 888.221 MiB, 1.89% gc time)
function update_KdV!(o::FFT_Data, u1, u0, Δt, niters)
v = similar(u0)
v_tmp = similar(u0)
u_tmp = similar(u0)
mul!(v, o.FFT, u0)
for iter in 1:niters
# v .= exp.(-Δt .* o.D3) .* v
v .= exp.(-Δt .* o.D3) .* v
# v .-= 3.0 .* Δt .* o.D .*(o.FFT * (o.FFT \ v).^2)
ldiv!(u_tmp, o.FFT, v)
u_tmp .*= u_tmp
mul!(v_tmp, o.FFT, u_tmp)
v .-= 3.0 .* Δt .* o.D .* v_tmp
end
ldiv!(u1, o.FFT, v)
end
function solve_KdV_inplace(o::FFT_Data{T}, f0, tmax) where T
Δt = 1/o.N^2
skip = floor(Int, 0.005/Δt)
t = 0:skip*Δt:tmax
M = length(t)
u = Array{Complex{T},2}(undef, o.N, M)
u[:,1] = Complex.(f0.(o.x))
for j in 2:M
update_KdV!(o, @view(u[:,j]), @view(u[:,j-1]), Δt, skip)
end
return real(u), t
end
N = 2^8
K = 20/(2π)
o = FFT_Data(K, N);
f0(x) = -5*sin(π*x/10)
Δt = 1/N^2
skip = 10*floor(Int, 0.005/Δt)
tmax = 1.0
@time u, t = solve_KdV_inplace(o, f0, tmax)
sleep(0.1)
plot_KdV(o, u, t)
7.794829 seconds (840.88 k allocations: 45.956 MiB, 0.20% gc time)
gifname = "KdV sin"
@time anim_KdV(gifname, o, u, t)
26.332842 seconds (3.50 M allocations: 181.681 MiB, 0.24% gc time)
using BenchmarkTools
N = 2^8
K = 20/(2π)
o = FFT_Data(K, N)
f0(x) = -5*sin(π*x/10)
Δt = 1/N^2
skip = 10*floor(Int, 0.005/Δt)
tmax = 1.4
1.4
@benchmark u, t = solve_KdV(o, f0, tmax)
BenchmarkTools.Trial: memory estimate: 1.10 GiB allocs estimate: 642611 -------------- minimum time: 11.433 s (2.25% GC) median time: 11.433 s (2.25% GC) mean time: 11.433 s (2.25% GC) maximum time: 11.433 s (2.25% GC) -------------- samples: 1 evals/sample: 1
@benchmark u, t = solve_KdV_inplace(o, f0, tmax)
BenchmarkTools.Trial: memory estimate: 20.48 MiB allocs estimate: 276931 -------------- minimum time: 10.609 s (0.00% GC) median time: 10.609 s (0.00% GC) mean time: 10.609 s (0.00% GC) maximum time: 10.609 s (0.00% GC) -------------- samples: 1 evals/sample: 1
以上のように in-place 版は非 in-place 版よりも計算時間が1割程度短く, in-place版の消費メモリは非 in-pace 版の50分の1である.
N = 2^8
K = 20/(2π)
o = FFT_Data(K, N);
KdVsoliton(c, a, x) = c/2*(sech(sqrt(c)/2*(x-a)))^2
f0(x) = KdVsoliton(16.0, -8.0, x) + KdVsoliton(4.0, -2.0, x)
Δt = 1/N^2
skip = 10*floor(Int, 0.005/Δt)
tmax = 1.0
@time u, t = solve_KdV_inplace(o, f0, tmax)
sleep(0.1)
plot_KdV(o, u, t, ymin=-0.5, ymax=10.0)
7.884885 seconds (384.44 k allocations: 24.193 MiB, 0.81% gc time)
gifname = "KdV 2-soliton"
@time anim_KdV(gifname, o, u, t, ymin=-0.5, ymax=10.0)
26.929848 seconds (89.47 k allocations: 9.704 MiB)
# V はポテンシャル函数
#
mutable struct Schroedinger_Data{T, S, R}
o::FFT_Data{T,S}
V::R
v::Array{T,1}
end
# FFT_Data o とポテンシャル函数 V から Schroedinger_Data を作成する函数
#
Schroedinger_Data(o, V) = Schroedinger_Data(o, V, V.(o.x))
function update_Schroedinger!(s::Schroedinger_Data, u1, u0, Δt, skip)
o = s.o
U_tmp = similar(u0)
u1 .= u0
for iter in 1:skip
# exp.(-im .* Δt .* s.v) .* (o.FFT\(exp.(0.5im .* Δt .* o.D2).*(o.FFT * u)))
mul!(U_tmp, o.FFT, u1)
U_tmp .*= exp.(0.5im .* Δt .* o.D2)
ldiv!(u1, o.FFT, U_tmp)
u1 .*= exp.(-im .* Δt .* s.v)
end
end
function solve_Schroedinger(s::Schroedinger_Data{T,S,R}, u0, tmax; Δt=2π/1000, skip=10) where {T,S,R}
o = s.o
t = 0:skip*Δt:tmax+10eps()
M = length(t)
u = Array{Complex{T},2}(undef, o.N, M)
u[:,1] = u0
for j in 2:M
update_Schroedinger!(s, @view(u[:,j]), @view(u[:,j-1]), Δt, skip)
end
return u, t
end
function plot_complexvalued(x, u, ymin, ymax)
plot(o.x, real(u), lw=0.8)
plot(o.x, imag(u), lw=0.8, ls="--")
ylim(ymin, ymax)
grid(ls=":")
xticks(fontsize=8)
yticks(fontsize=8)
end
function plot_Schroedinger(s::Schroedinger_Data, u, t; thin=10, ymin=-1.1, ymax=1.1)
nplots = (length(t)-1)÷thin + 1
o = s.o
nrows = div(nplots+3, 5)
figure(figsize=(4, 1.75))
subplot(121)
plot_complexvalued(o.x, u[:,1], ymin, ymax)
title("initial state", fontsize=9)
subplot(122)
plot_u(o.x, s.v)
title("potential", fontsize=9)
tight_layout()
figure(figsize=(10, 1.75nrows))
for j in 1:nplots-1
subplot(nrows, 5, j)
plot_complexvalued(o.x, u[:,(j-1)*thin+2], ymin, ymax)
title("j = $j", fontsize=9)
end
tight_layout()
end
function anim_Schroedinger(gifname, s::Schroedinger_Data, u, t; thin=2, interval=100, ymin=-1.1, ymax=1.1, giftitle=gifname)
o = s.o
function plotframe(j)
clf()
plot(o.x, real(u[:, (j-1)*thin+1]), lw=1.0)
plot(o.x, imag(u[:, (j-1)*thin+1]), lw=1.0, ls="--")
grid(ls=":")
ylim(ymin, ymax)
xticks(fontsize=8)
yticks(fontsize=8)
title("$(giftitle): j = $(j-1)", fontsize=10)
plot()
end
fig = figure(figsize=(4, 3))
n = (length(t)-1)÷thin + 1
frames = [1;1;1;1;1;1;1;1;1; 1:n; n;n;n;n;n;n;n;n;n;]
ani=animation[:FuncAnimation](fig, plotframe, frames=frames, interval=interval)
ani[:save]("$(gifname).gif", writer="imagemagick")
clf()
displayfile("image/gif", "$(gifname).gif")
end
anim_Schroedinger (generic function with 1 method)
mutable struct GaussianPacket{T<:Real}
x₀::T
p₀::T
σ::T
Δx::T
end
function GaussianPacket(x₀, p₀, σ, o::FFT_Data)
Δx = o.x[2] - o.x[1]
GaussianPacket(x₀, p₀, σ, Δx)
end
function (f::GaussianPacket)(x)
sqrt(f.Δx)/(π^(1/4)*sqrt(f.σ)) *
exp(im*f.p₀*(x-f.x₀/2) - (x-f.x₀)^2/(2f.σ^2))
end
N = 2^10
K = 10
o = FFT_Data(K, N);
g0 = GaussianPacket(-20.0, 5.0, 2.0, o)
plot_u(o.x, g0.(o.x));
# 自由粒子
N = 2^10
K = 10
o = FFT_Data(K, N);
@show typeof(o)
V(x) = 0.0
s = Schroedinger_Data(o, V)
g0 = GaussianPacket(-20.0, 5.0, 2.0, o)
u0 = g0.(o.x)
T = 2.0
tmax = 2π*T
@time u, t = solve_Schroedinger(s, u0, tmax)
sleep(0.1)
plot_Schroedinger(s, u, t, ymin=-0.2, ymax=0.2)
typeof(o) = FFT_Data{Float64,FFTW.cFFTWPlan{Complex{Float64},-1,false,1}} 0.809061 seconds (1.43 M allocations: 75.963 MiB, 3.53% gc time)
gifname = "free particle"
@time anim_Schroedinger(gifname, s, u, t, ymin=-0.2, ymax=0.2)
28.971086 seconds (405.89 k allocations: 36.768 MiB, 0.05% gc time)
# 壁=(高さ10、幅10)、粒子=運動量5)
N = 2^10
K = 10
o = FFT_Data(K, N);
@show typeof(o)
V(x) = ifelse(0 ≤ x ≤ 10, 10.0, 0.0)
s = Schroedinger_Data(o, V)
g0 = GaussianPacket(-20.0, 5.0, 1.0, o)
u0 = g0.(o.x)
T = 1.5
tmax = 2π*T
@time u, t = solve_Schroedinger(s, u0, tmax)
sleep(0.1)
plot_Schroedinger(s, u, t, ymin=-0.2, ymax=0.2)
typeof(o) = FFT_Data{Float64,FFTW.cFFTWPlan{Complex{Float64},-1,false,1}} 0.261462 seconds (5.03 k allocations: 5.007 MiB)
gifname = "low barrier"
@time anim_Schroedinger(gifname, s, u, t, ymin=-0.2, ymax=0.2)
22.670839 seconds (80.45 k allocations: 18.642 MiB, 0.03% gc time)
# 壁=(高さ20、幅10)、粒子=運動量5)
N = 2^10
K = 10
o = FFT_Data(K, N);
@show typeof(o)
V(x) = ifelse(0 ≤ x ≤ 10, 20.0, 0.0)
s = Schroedinger_Data(o, V)
g0 = GaussianPacket(-20.0, 5.0, 1.0, o)
u0 = g0.(o.x)
T = 1.5
tmax = 2π*T
@time u, t = solve_Schroedinger(s, u0, tmax)
sleep(0.1)
plot_Schroedinger(s, u, t, ymin=-0.2, ymax=0.2)
typeof(o) = FFT_Data{Float64,FFTW.cFFTWPlan{Complex{Float64},-1,false,1}} 0.298520 seconds (5.03 k allocations: 5.007 MiB, 2.52% gc time)
gifname = "high barrier"
@time anim_Schroedinger(gifname, s, u, t, ymin=-0.2, ymax=0.2)
23.436036 seconds (80.00 k allocations: 14.509 MiB, 0.03% gc time)
# 調和振動子 (周期 2π)
N = 2^10
K = 10
o = FFT_Data(K, N);
@show typeof(o)
V(x) = x^2/2
s = Schroedinger_Data(o, V)
g0 = GaussianPacket(-20.0, 0.0, 0.5, o)
u0 = g0.(o.x)
T = 2.0
tmax = 2π*T
@time u, t = solve_Schroedinger(s, u0, tmax)
sleep(0.1)
plot_Schroedinger(s, u, t, ymin=-0.3, ymax=0.3)
typeof(o) = FFT_Data{Float64,FFTW.cFFTWPlan{Complex{Float64},-1,false,1}} 0.415711 seconds (6.68 k allocations: 6.664 MiB)
gifname = "harmonic oscillator"
@time anim_Schroedinger(gifname, s, u, t, ymin=-0.3, ymax=0.3, thin=1)
52.132494 seconds (191.10 k allocations: 30.885 MiB, 0.06% gc time)
# V(x) = -20 exp(-x^2/25)
N = 2^10
K = 10
o = FFT_Data(K, N);
@show typeof(o)
V(x) = -100 * exp(-x^2/100)
s = Schroedinger_Data(o, V)
g0 = GaussianPacket(-15.0, 0.0, 1.0, o)
u0 = g0.(o.x)
T = 4.0
tmax = 2π*T
@time u, t = solve_Schroedinger(s, u0, tmax, skip=20)
sleep(0.1)
plot_Schroedinger(s, u, t, ymin=-0.2, ymax=0.2)
typeof(o) = FFT_Data{Float64,FFTW.cFFTWPlan{Complex{Float64},-1,false,1}} 0.955456 seconds (115.11 k allocations: 12.158 MiB, 6.97% gc time)
gifname = "V(x)=-100exp(-x2/100)"
giftitle = "\$V(x)=-100\\exp(-100x^2)\$"
@time anim_Schroedinger(gifname, s, u, t, ymin=-0.2, ymax=0.2, giftitle=giftitle)
30.274456 seconds (115.35 k allocations: 22.111 MiB, 0.02% gc time)
KdV equation: $$ u_t = -\partial_x^3 u - 3\partial_x(u^2). $$
Smith's equation: $$ u_t = 2a^{-2}\left( \sqrt{1-a^2\partial_x^2} - 1 \right)\partial_x u - 3\partial_x(u^2). $$
以下では $a^2=0.15$ と仮定する.
function plot_Smith(o::FFT_Data, u, t; thin=10, ymin=-7.0, ymax=25.0)
nplots = div(length(t), thin)
nrows = div(nplots+4, 5)
figure(figsize=(10, 1.75nrows))
j = 0
plot_u(o, @view(u[:,j+1]), nrows, j, ymin, ymax)
for j in 1:nplots-1
plot_u(o, @view(u[:,thin*j+1]), nrows, j, ymin, ymax)
end
tight_layout()
end
function anim_Smith(gifname, o::FFT_Data, u, t; thin=2, interval=100, ymin=-7.0, ymax=25.0)
function plotframe(j)
clf()
plot(o.x, real(u[:, (j-1)*thin + 1]), lw=1.0)
grid(ls=":")
ylim(ymin, ymax)
xticks(fontsize=8)
yticks(fontsize=8)
title("$(gifname): j = $(j-1)", fontsize=10)
plot()
end
fig = figure(figsize=(4, 3))
n = (length(t)-1)÷thin + 1
frames = [1;1;1;1;1;1;1;1;1; 1:n; n;n;n;n;n;n;n;n;n;]
ani=animation[:FuncAnimation](fig, plotframe, frames=frames, interval=interval)
ani[:save]("$(gifname).gif", writer="imagemagick")
clf()
displayfile("image/gif", "$(gifname).gif")
end
anim_Smith (generic function with 1 method)
function update_Smith!(o::FFT_Data, u1, u0, Δt, niters, a²)
v = similar(u0)
v_tmp = similar(u0)
u_tmp = similar(u0)
mul!(v, o.FFT, u0)
for iter in 1:niters
v .= exp.(2 ./ a² .* Δt .* (sqrt.(1 .- a² .* o.D2) .- 1).*o.D) .* v
# v .-= 3.0 .* Δt .* o.D .*(o.FFT * (o.FFT \ v).^2)
ldiv!(u_tmp, o.FFT, v)
u_tmp .*= u_tmp
mul!(v_tmp, o.FFT, u_tmp)
v .-= 3.0 .* Δt .* o.D .* v_tmp
end
ldiv!(u1, o.FFT, v)
end
function solve_Smith_inplace(o::FFT_Data{T}, f0, tmax; a²=0.15) where T
Δt = 1/o.N^2
skip = floor(Int, 0.005/Δt)
t = 0:skip*Δt:tmax
M = length(t)
u = Array{Complex{T},2}(undef, o.N, M)
u[:,1] = Complex.(f0.(o.x))
for j in 2:M
update_Smith!(o, @view(u[:,j]), @view(u[:,j-1]), Δt, skip, a²)
end
return real(u), t
end
solve_Smith_inplace (generic function with 1 method)
N = 2^9
K = 20/(2π)
o = FFT_Data(K, N);
f0(x) = -5*sin(π*x/10)
Δt = 1/N^2
skip = 10*floor(Int, 0.005/Δt)
tmax = 1.0
@time u, t = solve_Smith_inplace(o, f0, tmax)
sleep(0.1)
plot_Smith(o, u, t, thin=10)
40.363118 seconds (2.11 M allocations: 115.233 MiB, 0.12% gc time)
gifname = "Smith sin"
@time anim_Smith(gifname, o, u, t)
27.640708 seconds (360.26 k allocations: 24.151 MiB, 0.03% gc time)
N = 2^9
K = 20/(2π)
o = FFT_Data(K, N);
KdVsoliton(c, a, x) = c/2*(sech(sqrt(c)/2*(x-a)))^2
f0(x) = KdVsoliton(16.0, -8.0, x) + KdVsoliton(4.0, -2.0, x)
Δt = 1/N^2
skip = 10*floor(Int, 0.005/Δt)
tmax = 1.0
@time u, t = solve_Smith_inplace(o, f0, tmax)
sleep(0.1)
plot_Smith(o, u, t, ymin=-0.5, ymax=12.0)
38.722846 seconds (912.37 k allocations: 57.235 MiB, 0.04% gc time)
gifname = "Smith 2-soliton"
@time anim_Smith(gifname, o, u, t, ymin=-0.5, ymax=12.0)
26.687429 seconds (88.84 k allocations: 10.016 MiB)