Комплекс ФДЭ (Гинзбург Ландау) в Юлии с помощью псевдоспектрального метода

Question

Я хочу научиться решать ПДД с Джулией, и сейчас я пытаюсь решить комплексное уравнение Гинзбурга Ландау (CGLE) с помощью псевдоспектрального метода в Джулии. Тем не менее, я борюсь с этим, и у меня есть немного идей, что попробовать.

CGLE гласит:

С помощью преобразования Фурье и его обратного я могу преобразовать в спектральную форму:

Это, например, также дано в этом старом скрипте, который я нашел (https://www.uni -muenster.de / Physik.TP / archive / fileadmin / lehre / NumMethoden / SoSe2009 / Skript / script.pdf ) Из того же источника я знаю, что альфа = 1, бета = 2 и начальные условия с небольшим шумом порядка 0,01 вокруг 0 ​​должны приводить к плоским волнам в качестве решений. Вот что я хочу проверить в первую очередь.

Следуя очень хорошему уроку Криса Ракауцкаса (https://youtu.be/okGybBmihOE),, я попытался использовать ApproxFun и diffrentialEquations следующим образом для решения этой проблемы:

РЕДАКТИРОВАТЬ: я исправил две ошибки из исходного поста: пропущенная точка и знак минус, но код все еще не дает правильных результатов

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 mathemitacaly, ни по их численной обработке, так что я до сих пор в основном работал с 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 оказался чрезвычайно медленным решателем для этого случая, просто использование по умолчанию хорошо работает для меня.

Любая помощь приветствуется.

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.