julia - Julia 中的伪谱方法的复偏微分方程 (Ginzburg Landau)

Question

我想自学如何使用 Julia 求解偏微分方程，并且现在正在尝试使用 Julia 中的伪谱方法求解复杂的 Ginzburg Landau 方程 (CGLE)。然而，我对此很挣扎，我有点想尝试什么。

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 来解决此问题:

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

编辑2:发现我计算的波数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，但结果是相同的。我的问题是，我既不是偏微分方程数学专家，也不是数值处理方面的专家，到目前为止，我主要研究常微分方程。

EDIT2 现在看来这或多或少有效。但我仍然想知道一些事情

• 如何在指定域(例如 )上计算此值，我对 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:在这种情况下，罗达斯被证明是一个极其缓慢的求解器，仅使用默认值对我来说效果很好。

感谢任何帮助。

Answer 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。