julia 中的抛物线偏微分方程

Question

我正在尝试使用 Julia 以数值方式求解抛物线偏微分方程，但我找不到任何可以提供帮助的可访问文档。

这里有一个例子：t, x 是一维实数。我想解决 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)

在 Julia 中可以这样做吗？这种类型的问题可以在 Matlab 中使用 pdepe 解决。

Answer 1

现在，我们没有“完全停止”的 PDE 求解器，即您输入 PDE 并运行的求解器。但是，PDE 是通过离散化为 ODE 来求解的，因此为此编写句号 PDE 求解器的方式如下。这大部分是在 this blog post BTW 中进行了更深入的讨论。

带上你的偏微分方程。现在用 dx 离散化运算符。 LaPlacian 的二阶有限差分离散化是应用于我们状态的模板 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 执行离散化拉普拉斯算子。

一阶导数项的处理方法类似，但在本例中通常使用 upwinding scheme。您可以将 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

注意我没有做末端因为我不知道你想要什么边界条件。如果它是具有零常数的 Dirichlet，那么您只需将其写在端点处，但去掉 space 之外的值。这里 x[i] = x0 + dx*i.

现在您有一组 ODE，其中 u[i,j] = u_j(x_i)。因此，您将初始条件离散化为 u0[i,j] 并设置 ODE 问题：

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

For this, see the DiffEq documentation, specifically the ODE tutorial。现在您只需求解 PDE 的离散化 ODE 表示。对于这类方程，如博文post中提到的，Sundials.jl CVODE_BDF方法与GMRES Krylov线性求解器是一个不错的选择，所以我们这样做：

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

返回一个连续解，其中 sol(t)[i,j] 是 u_j(t,x_i) 的数值近似值。当然，较低的 dx 更准确，您应该根据需要调整 ODE 求解器的容差。

我们将在不久的将来为偏微分方程自动执行此操作（任何顺序的任何导数），但是it's currently a work-in-progress so for now manual discretization will have to do (and this is taught in any numerical methods course, so it's not too bad!). Hope this helps. If you need any more help, please check out our chat channel因为那里的大多数人都有进行这种离散化的经验。

julia 中的抛物线偏微分方程

Parabolic PDEs in julia

julia

differentialequations.jl