function kdv2a()
N = 128
x = collect(linspace(-10,10,2N))
delta_x = x[2] - x[1]
delta_k = pi/(N*delta_x)
delta_t = 1/(2N)^2
k = delta_k*[0:N; -(N-1):-1]
c=16
u = 1/2*c*(sech.(sqrt(c)/2*(x.+8.0))).^2
c=9
u += 1/2*c*(sech.(sqrt(c)/2*(x.+4.0))).^2
tstep = 0.005
nstep = Integer(floor(tstep/delta_t))
#T = Integer(floor(1/tstep))
T = 180
t = 10*tstep*collect(0:T)
kdv = Array{Float64,2}(length(u), T+1)
kdv[:,1] = u
for s=1:T
U=fft(kdv[:,s])
for n = 1:nstep
U = U .* exp.(im*k.^3*delta_t)
U = U - delta_t*(3*im*k.*fft(real(ifft(U)).^2))
end
kdv[:,s+1]=real(ifft(U))
end
return kdv, x, t
end
kdv2a (generic function with 1 method)
@time u, x, t = kdv2a()
5.590454 seconds (8.56 M allocations: 3.105 GiB, 9.60% gc time)
([0.0107279 0.00662186 … -0.00189289 -0.00147245; 0.0146776 0.00944989 … -0.00117761 -0.0012299; … ; 1.30946e-17 0.00191254 … -0.000701666 0.000129507; 1.03491e-17 0.0035073 … -0.000308399 -0.000774489], [-10.0, -9.92157, -9.84314, -9.76471, -9.68627, -9.60784, -9.52941, -9.45098, -9.37255, -9.29412 … 9.29412, 9.37255, 9.45098, 9.52941, 9.60784, 9.68627, 9.76471, 9.84314, 9.92157, 10.0], [0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45 … 8.55, 8.6, 8.65, 8.7, 8.75, 8.8, 8.85, 8.9, 8.95, 9.0])
using PyPlot
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T
m, n = length(vy), length(vx)
vx = reshape(vx, 1, n)
vy = reshape(vy, m, 1)
(repmat(vx, m, 1), repmat(vy, 1, n))
end
tt, xx = meshgrid(t,x)
pcolormesh(xx, tt, u, cmap="rainbow")
grid(ls=":")
filename = "kdv2a"
withfig(figure(figsize=(6.4,3.6))) do
k = 1
n = size(u,2)
is = [1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1; 1:n; n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;]
maxu = maximum(u)
minu = minimum(u)
for i in is
clf()
plot(x, u[:,i])
grid(ls=":")
xlim(-10,10)
ylim(minu-0.2, maxu+0.2)
tstr = @sprintf("%.2f", t[i])
title("KdV: t = $tstr")
savefig(filename * @sprintf("_%04d",k), bbox_inches="tight")
k += 1
end
end
run(`convert -delay 5 -loop 0 $(filename)_\*.png $filename.gif`)
open(filename*".gif") do f
base64_video = base64encode(f)
display("text/html", """<img src="data:image/gif;base64,$base64_video">""")
end
run(`rm $(filename)_\*.png`)
error in running finalizer: Base.ArgumentError(msg="ref of NULL PyObject")
function kdv2b()
N = 128
x = collect(linspace(-10,10,2N))
delta_x = x[2] - x[1]
delta_k = pi/(N*delta_x)
delta_t = 1/(2N)^2
k = delta_k*[0:N; -(N-1):-1]
c=16
u = 1/2*c*(sech.(sqrt(c)/2*(x.+8.0))).^2
c=4
u += 1/2*c*(sech.(sqrt(c)/2*(x.+2.0))).^2
tstep = 0.005
nstep = Integer(floor(tstep/delta_t))
#T = Integer(floor(1/tstep))
T = 180
t = 10*tstep*collect(0:T)
kdv = Array{Float64,2}(length(u), T+1)
kdv[:,1] = u
for s=1:T
U=fft(kdv[:,s])
for n = 1:nstep
U = U .* exp.(im*k.^3*delta_t)
U = U - delta_t*(3*im*k.*fft(real(ifft(U)).^2))
end
kdv[:,s+1]=real(ifft(U))
end
return kdv, x, t
end
kdv2b (generic function with 1 method)
@time u, x, t = kdv2b()
4.149392 seconds (7.70 M allocations: 3.066 GiB, 13.60% gc time)
([0.0107285 0.00663015 … -0.0014681 -0.000545773; 0.0146783 0.00943505 … -0.00126715 -0.00191005; … ; 3.53303e-10 0.001905 … -0.00181044 -0.00132111; 3.02011e-10 0.00350786 … -0.00176343 -0.00106325], [-10.0, -9.92157, -9.84314, -9.76471, -9.68627, -9.60784, -9.52941, -9.45098, -9.37255, -9.29412 … 9.29412, 9.37255, 9.45098, 9.52941, 9.60784, 9.68627, 9.76471, 9.84314, 9.92157, 10.0], [0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45 … 8.55, 8.6, 8.65, 8.7, 8.75, 8.8, 8.85, 8.9, 8.95, 9.0])
using PyPlot
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T
m, n = length(vy), length(vx)
vx = reshape(vx, 1, n)
vy = reshape(vy, m, 1)
(repmat(vx, m, 1), repmat(vy, 1, n))
end
tt, xx = meshgrid(t,x)
pcolormesh(xx, tt, u, cmap="rainbow")
grid(ls=":")
filename = "kdv2b"
withfig(figure(figsize=(6.4,3.6))) do
k = 1
n = size(u,2)
is = [1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1; 1:n; n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;]
maxu = maximum(u)
minu = minimum(u)
for i in is
clf()
plot(x, u[:,i])
grid(ls=":")
xlim(-10,10)
ylim(minu-0.2, maxu+0.2)
tstr = @sprintf("%.2f", t[i])
title("KdV: t = $tstr")
savefig(filename * @sprintf("_%04d",k), bbox_inches="tight")
k += 1
end
end
run(`convert -delay 5 -loop 0 $(filename)_\*.png $filename.gif`)
open(filename*".gif") do f
base64_video = base64encode(f)
display("text/html", """<img src="data:image/gif;base64,$base64_video">""")
end
run(`rm $(filename)_\*.png`)
function kdvsin_nopf()
N = 128
x = collect(linspace(-10,10,2N))
delta_x = x[2] - x[1]
delta_k = pi/(N*delta_x)
delta_t = 1/(2N)^2
k = delta_k*[0:N; -(N-1):-1]
u0 = -5*sin.(pi*x/10)
tstep = 0.005
nstep = Integer(floor(tstep/delta_t))
#T = Integer(floor(1/tstep))
T = 180
t = 10*tstep*collect(0:T)
u = Array{Float64,2}(length(u0), T+1)
u[:,1] = u0
U = similar(u0)
for s in 1:T
U = fft(@view u[:,s])
for n in 1:nstep
U = U .* exp.(im .* k.^3 .* delta_t)
U = U .- delta_t .*(3im .* k.*(fft(real(ifft(U).^2))))
end
u[:,s+1] = real(ifft(U))
end
return u, x, t
end
kdvsin_nopf (generic function with 1 method)
@time u, x, t = kdvsin_nopf()
3.459570 seconds (7.46 M allocations: 1.906 GiB, 9.35% gc time)
([6.12323e-16 0.0300603 … 1.18984 2.8248; 0.123187 0.144892 … 0.578729 1.89663; … ; -0.123187 -0.174981 … 4.6908 6.62513; -6.12323e-16 -0.0683861 … 2.65886 4.37568], [-10.0, -9.92157, -9.84314, -9.76471, -9.68627, -9.60784, -9.52941, -9.45098, -9.37255, -9.29412 … 9.29412, 9.37255, 9.45098, 9.52941, 9.60784, 9.68627, 9.76471, 9.84314, 9.92157, 10.0], [0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45 … 8.55, 8.6, 8.65, 8.7, 8.75, 8.8, 8.85, 8.9, 8.95, 9.0])
using PyPlot
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T
m, n = length(vy), length(vx)
vx = reshape(vx, 1, n)
vy = reshape(vy, m, 1)
(repmat(vx, m, 1), repmat(vy, 1, n))
end
tt, xx = meshgrid(t,x)
pcolormesh(xx, tt, u, cmap="rainbow")
grid(ls=":")
function kdvsin()
N = 128
x = collect(linspace(-10,10,2N))
delta_x = x[2] - x[1]
delta_k = pi/(N*delta_x)
delta_t = 1/(2N)^2
k = delta_k*[0:N; -(N-1):-1]
u0 = -5*sin.(pi*x/10)
tstep = 0.005
nstep = Integer(floor(tstep/delta_t))
#T = Integer(floor(1/tstep))
T = 180
t = 10*tstep*collect(0:T)
u = Array{Float64,2}(length(u0), T+1)
u[:,1] = u0
P = plan_fft(u0)
U = similar(u0)
for s in 1:T
U = P*(@view u[:,s])
for n in 1:nstep
U = U .* exp.(im .* k.^3 .* delta_t)
U = U .- delta_t .*(3im .* k.*(P*real(P\U).^2))
end
u[:,s+1] = real(P\U)
end
return u, x, t
end
kdvsin (generic function with 1 method)
@time u, x, t = kdvsin()
1.465380 seconds (894.55 k allocations: 1.424 GiB, 8.54% gc time)
([6.12323e-16 0.0300603 … 1.18984 2.8248; 0.123187 0.144892 … 0.578729 1.89663; … ; -0.123187 -0.174981 … 4.6908 6.62513; -6.12323e-16 -0.0683861 … 2.65886 4.37568], [-10.0, -9.92157, -9.84314, -9.76471, -9.68627, -9.60784, -9.52941, -9.45098, -9.37255, -9.29412 … 9.29412, 9.37255, 9.45098, 9.52941, 9.60784, 9.68627, 9.76471, 9.84314, 9.92157, 10.0], [0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45 … 8.55, 8.6, 8.65, 8.7, 8.75, 8.8, 8.85, 8.9, 8.95, 9.0])
using PyPlot
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T
m, n = length(vy), length(vx)
vx = reshape(vx, 1, n)
vy = reshape(vy, m, 1)
(repmat(vx, m, 1), repmat(vy, 1, n))
end
tt, xx = meshgrid(t,x)
pcolormesh(xx, tt, u, cmap="rainbow")
grid(ls=":")
filename = "kdvsin"
withfig(figure(figsize=(6.4,3.6))) do
k = 1
n = size(u,2)
is = [1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1; 1:n; n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;n;]
maxu = maximum(u)
minu = minimum(u)
for i in is
clf()
plot(x, u[:,i])
grid(ls=":")
xlim(-10,10)
ylim(minu-0.2, maxu+0.2)
tstr = @sprintf("%.2f", t[i])
title("KdV: t = $tstr")
savefig(filename * @sprintf("_%04d",k), bbox_inches="tight")
k += 1
end
end
run(`convert -delay 5 -loop 0 $(filename)_\*.png $filename.gif`)
open(filename*".gif") do f
base64_video = base64encode(f)
display("text/html", """<img src="data:image/gif;base64,$base64_video">""")
end
run(`rm $(filename)_\*.png`)
function kdv(x::Array{Float64}, t::Array{Float64}, u0::Array{Float64})
N = length(x)
@assert iseven(N) "length(x) is not even."
@assert N == length(u0) "length(u0) is not equal to length(x)."
delta_x = x[2] - x[1]
delta_k = 2*pi/(N*delta_x)
k = delta_k*[0:div(N,2); -(div(N,2)-1):-1]
T = length(t)
delta_t = 1/N^2
tstep = (t[2] - t[1])/(N*delta_x/2)
nstep = Integer(floor(N^2*tstep))
u = Array{Float64,2}(length(u0), T)
u[:,1] = u0
P = plan_fft(u0)
U = similar(u0)
for s in 1:T-1
U = P*(@view u[:,s])
for n in 1:nstep
U = U .* exp.(im .* k.^3 .* delta_t)
U = U .- delta_t .*(3im .* k.*(P*real(P\U).^2))
end
u[:,s+1] = real(P\U)
end
return u
end
kdv (generic function with 1 method)
N = 2^8
x = collect(linspace(-10,10,N))
t = collect(linspace(0,9,180))
u0 = -5*sin.(pi*x/10)
@time u = kdv(x,t,u0);
1.118376 seconds (689.65 k allocations: 1.410 GiB, 8.58% gc time)
using PyPlot
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T
m, n = length(vy), length(vx)
vx = reshape(vx, 1, n)
vy = reshape(vy, m, 1)
(repmat(vx, m, 1), repmat(vy, 1, n))
end
tt, xx = meshgrid(t,x)
figure(figsize=(6.4,4))
pcolormesh(xx, tt, u, cmap="rainbow")
colorbar()
grid(ls=":")
N = 2^9
x = collect(linspace(-10,10,N))
t = collect(linspace(0,9,360))
u0 = -5*sin.(pi*x/10)
@time u = kdv(x,t,u0);
7.698346 seconds (2.35 M allocations: 11.010 GiB, 7.69% gc time)
using PyPlot
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T
m, n = length(vy), length(vx)
vx = reshape(vx, 1, n)
vy = reshape(vy, m, 1)
(repmat(vx, m, 1), repmat(vy, 1, n))
end
tt, xx = meshgrid(t,x)
figure(figsize=(6.4,4))
pcolormesh(xx, tt, u, cmap="rainbow")
colorbar()
grid(ls=":")