2

我想自学如何用 Julia 求解 PDE，我现在正在尝试用 Julia 中的伪光谱方法求解复杂的 Ginzburg Landau 方程 (CGLE)。但是，我对此感到很挣扎，并且我对尝试什么有点想法。

CGLE 内容如下： cgle

通过傅里叶变换F及其逆变换F-1，我可以变换为谱形式：

cgle光谱

例如，在我发现的这个旧脚本中也给出了这一点（https://www.uni-muenster.de/Physik.TP/archive/fileadmin/lehre/NumMethoden/SoSe2009/Skript/script.pdf）来自同一个来源知道，alpha = 1，beta = 2 和在 0 附近具有 0.01 阶小噪声的初始条件应该导致平面波作为解决方案。这就是我想首先测试的。

遵循 Chris Rackauckas ( https://youtu.be/okGybBmihOE ) 的非常好的教程，我尝试通过以下方式使用 ApproxFun 和 DifferentialEquations 来解决这个问题：

编辑：我更正了原始帖子中的两个错误，一个缺失的点和减号，但代码仍然没有给出正确的结果

EDIT2：发现我计算的波数k完全错误 

using ApproxFun
using DifferentialEquations

F = Fourier()
n = 512
L = 100
T = ApproxFun.plan_transform(F, n)
Ti = ApproxFun.plan_itransform(F, n)
x = collect(range(-L/2,stop=L/2, length=n))
k = points(F, n)

alpha = 1im
beta = 2im
u0 = 0.01*(rand(ComplexF64, n) .- 0.5)
Fu0 = T*u0

function cgle!(du, u, p, t)
    a, b, k, T, Ti = p

    invu = Ti*u
    du .= (1.0 .- k.^2*(1.0 .+a)).*u .- T*( (1.0 .+b) .* (abs.(invu)).^2 .* invu)
end

pars = alpha, beta, k, T, Ti
prob = ODEProblem(cgle!, Fu0, (0.,50.), pars)
u = solve(prob, Rodas5(autodiff=false))

# plotting on a equidistant time stepping
t = collect(range(0, stop=50, length=1000))
sol = zeros(eltype(u),(n, length(t)))
for it in eachindex(t)
   sol[:,it] = Ti*u(t[it])
end

IM = PyPlot.imshow(real.(sol))
cb = PyPlot.colorbar(IM, orientation="horizontal")
gcf()

（已编辑）我尝试了不同的求解器，正如视频中所推荐的那样，有些显然不适用于复数，有些则可以，但是当我运行此代码时，它并没有给出预期的结果。该解决方案的价值仍然很小，并且不会产生实际上应该是结果的平面波。我还测试了其他会导致混乱的初始条件，但这些也会导致相同的非常小的解决方案。我还交替使用了 ApproxFun 的显式拉普拉斯算子，但结果是一样的。我的问题是，我既不是 PDE 数学的真正专家，也不是他们的数值处理专家，到目前为止，我主要使用 ODE。

EDIT2这现在似乎或多或少地起作用。我仍然想知道一些事情

  • 我如何在指定的域上计算它x\in[-L/2;L/2]，我对它如何与 ApproxFun 一起工作感到非常困惑，据我所知波数k应该是(2pi/L)*[-N/2+1 ; N/2 -1]，但我不太确定如何用 ApproxFun 做到这一点
  • https://codeinthehole.com/tutorial/coherent.html显示了方程的不同动态状态/相图。虽然我可以重现其中一些，但有些似乎不起作用，比如时空间歇性

编辑 3：我通过直接使用 FFTW 而不是 ApproxFun 解决了这些问题。如果有人知道如何使用 ApproxFun，我仍然会感兴趣。下面是 FFTW 的代码（它也对性能进行了一些优化）

begin
   using FFTW
   using DifferentialEquations
   using PyPlot
end

begin
   n = 512
   L = 200
   n2 = Int(n/2)
   alpha = 2im
   beta = 1im
   x = range(-L/2,stop=L/2,length=n)
   u0 = 0.01*(rand(ComplexF64, n) .- 0.5)

   k = [0:n/2-1; 0; -n/2+1:-1] .*(2*pi/L);
   k2 = k.*k
   k2[n2 + 1] = (n2*(2*pi/L))^2 

   T = plan_fft(u0)
   Ti = plan_ifft(T*u0)

   LinOp = (1.0 .- k2.*(1.0 .+alpha))
   Fu0 = T*u0
end

function cgle!(du, u, p, t)
    LinOp, b, T, Ti = p

    invu = Ti*u
    du .= LinOp.*u .- T*( (1.0 .+b) .* (abs.(invu)).^2 .* invu)
end

pars = LinOp, beta, T, Ti
prob = ODEProblem(cgle!, Fu0, (0.,100.), pars)
@time u = solve(prob)

t = collect(range(0, stop=50, length=1000))
sol = zeros(eltype(u),(n, length(t)))
for it in eachindex(t)
   sol[:,it] = Ti*u(t[it])
end

IM = PyPlot.imshow(abs.(sol))
cb = PyPlot.colorbar(IM, orientation="horizontal") 
gcf()

编辑 4：对于这种情况，Rodas 是一个极其缓慢的求解器，只使用默认值对我来说效果很好。

任何帮助表示赞赏。

4

1 回答 1

1 
du = (1. .- k.^2*(1. .+(im*a))).*u + T*( (1. .+(im*b)) .* abs.(invu).^2 .* invu)

请注意，这是替换指向 du 的指针，而不是更新它。改用类似的东西.=：

du .= (1. .- k.^2*(1. .+(im*a))).*u + T*( (1. .+(im*b)) .* abs.(invu).^2 .* invu)

否则你的导数只是 0。

于 2019-07-03T18:58:47.343 回答