julia - 伪谱方法在Julia中生成复杂PDE(Ginzburg Landau)

我想教自己如何使用Julia解决PDE，现在我正在尝试使用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.01的小噪声大约为0）会导致平面波作为解。那就是我要首先测试的。

遵循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的显式Laplace运算符，但结果是相同的。我在这里的问题是，我既不是真正的PDE数学专家，也不是其数值处理专家，到目前为止，我主要从事ODE的研究。

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：在这种情况下，Rodas确实是一个极慢的求解器，只是使用默认值对我来说效果很好。

任何帮助表示赞赏。

最佳答案

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。

关于julia - 伪谱方法在Julia中生成复杂PDE(Ginzburg Landau)，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/56874260/

approxfun

julia - 伪谱方法在Julia中生成复杂PDE(Ginzburg Landau)