Параболические PDE в Юлии

Question

Я пытаюсь решить параболическое дифференциальное уравнение в частных производных численно, используя Джулию, но я не могу найти доступную документацию, которая может помочь.

Вот пример: t, x являются 1-мерным вещественным веществом. Я хочу решить для u(t,x)=[u1(t,x) u2(t,x)]; ты удовлетворяешь PDE

du1 / dt = d ^ 2u1 / dx ^ 2 + a11 (x, u) du1 / dx + a12 (x, u) du2 / dx + c1 (x, u)

du2 / dt = d ^ 2u2 / dx ^ 2 + a21 (x, u) du1 / dx + a22 (x, u) du2 / dx + c2 (x, u)

Возможно ли это сделать в Юлии? Проблема такого типа решаема с помощью pdepe в Matlab.

Answer 1

В настоящее время у нас нет "полноправных" решателей PDE, то есть решателей, куда вы вставляете PDE и отправляетесь. Тем не менее, PDE решаются путем дискретизации с ODE, поэтому способ, которым для этого был бы написан решатель PDE с полной остановкой, заключается в следующем. Большая часть этого обсуждается более подробно в этом сообщении в блоге.

Возьми свой PDE. Теперь дискретизируем операторов по dx, Конечно-разностная дискретизация лапласиана второго порядка - трафарет u[i-1] - 2 u[i] + u[i+1] применяется к нашему государству. Тогда, конечно, вы должны учитывать ваши граничные условия, когда вы дойдете до конца. Обычно хороший способ написать это в виде матрицы, поэтому:

const Mx = Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1])
# Do the reflections, different for x and y operators
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0

имеет Mx*u/dx^2 выполнить дискретизированный LaPlacian.

Производные члены первого порядка выполняются аналогично, но в этом случае обычно используется схема с намоткой. Вы можете взять свой du1/dx термин и заменить его ядром

a[i]*(u[i]-u[i-1])/dx

когда a положительный, или

a[i]*(u[i]-u[i+1])/dx

когда a отрицательно. Затем включите граничные условия, конечно. Тогда вы просто пишете свои реакции как c1(x[i],u[i]), Вместе это будет выглядеть (в нематричной форме:

function f(t,u,du)
    u1 = @view u[:,1]
    u2 = @view u[:,1]
    du1 = @view du[:,1]
    du2 = @view du[:,2]
    for i in 2:length(u)-1
        du1[i] = (u1[i-1] - 2u1[i] + u1[i+1])/dx^2 +
                a11(x[i],u1[i])*(u1[i]-u1[i-1])/dx +
                a12(x[i],u1[i])*(u1[i]-u1[i-1])/dx +
                c1(x1[i],u1[i])

        du2[i] = (u2[i-1] - 2u2[i] + u2[i+1])/dx^2 +
                a11(x[i],u2[i])*(u2[i]-u2[i-1])/dx +
                a12(x[i],u2[i])*(u2[i]-u2[i-1])/dx +
                c1(x1[i],u2[i])
    end

end

Обратите внимание, что я не делал концы, потому что я не знаю, какие граничные условия вы хотите. Если это Дирихле с нулевой константой, то вы просто пишете это в конечных точках, но отбрасываете значения, которые не заполнены пробелом. И здесь x[i] = x0 + dx*i,

Теперь у вас есть набор ODE, где u[i,j] = u_j(x_i), Таким образом, вы преобразуете свое первоначальное состояние в u0[i,j] и настройте проблему ODE:

using DifferentialEquations
prob = ODEProblem(f,u0,tspan)

Для этого смотрите документацию по DiffEq, в частности учебник по ODE. Теперь вы просто решаете дискретное представление ODE PDE. Для этих видов уравнений, как упоминалось в сообщении в блоге, Sundials.jl CVODE_BDF метод с GMRES Линейный решатель Крылова - хороший выбор, поэтому мы делаем:

sol = solve(prob,CVODE_BDF(linear_solver=:GMRES))

Это возвращает непрерывное решение, где sol(t)[i,j] численное приближение к u_j(t,x_i), Конечно ниже dx является более точным, и вы должны отрегулировать допуски решателя ODE, как считаете нужным.

В ближайшем будущем у нас будет возможность сделать это автоматически для PDE (любые производные по любому порядку), но в настоящее время это незавершенное производство, поэтому сейчас придется делать дискретизацию вручную (и этому учат в любом курсе численных методов). так что это не так уж плохо!). Надеюсь это поможет. Если вам нужна дополнительная помощь, пожалуйста, посетите наш канал чата, так как большинство людей там будут иметь опыт такого рода дискретизации.