Julia中的抛物线偏微分方程答案

【问题标题】：Parabolic PDEs in juliaJulia中的抛物线偏微分方程
【发布时间】：2018-01-17 23:17:59
【问题描述】：

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

这是一个例子：t, x 是一维实数。我想解决 u(t,x)=[u1(t,x) u2(t,x)];你满足偏微分方程

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 解决。

【问题讨论】：

你查看过JuliaDiffEq的DifferentialEquations.jl吗？
嗨 Ismael，谢谢 - 是的，我已经看过这个，但我无法弄清楚如何在 PDE 中使用梯度项，即如何求解上述方程，其中 aij 不等于零。也没有提供示例，所以我不确定这个包是否可行。 EllipticFEM 包看起来很有希望，但文档非常非常稀疏（至少对我而言）
我建议也去看看 DiffEq 聊天室：gitter.im/JuliaDiffEq/Lobby

标签： julia differentialequations.jl

【解决方案1】：

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

带上你的 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 执行离散化拉普拉斯算子。

一阶导数项的处理方式类似，但在这种情况下，通常使用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

请注意，我没有做结尾，因为我不知道你想要什么边界条件。如果它是具有零常数的狄利克雷，那么您只需将其写在端点处，但删除超出空间的值。在这里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 表示。对于这些类型的方程，正如博文中提到的，带有GMRES Krylov 线性求解器的 Sundials.jl CVODE_BDF 方法是一个不错的选择，所以我们这样做了：

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

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

我们将在不久的将来为 PDE 自动执行此操作（任何阶的任何导数），但 it's currently a work-in-progress 所以现在必须进行手动离散化（这在任何数值方法课程中都有教授，所以这还不错！）。希望这可以帮助。如果您需要更多帮助，请check out our chat channel，因为那里的大多数人都会有这种离散化的经验。

【讨论】：

感谢您的深入回答——很清楚。