7

我正在尝试使用 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 可以解决此类问题。

4

1 回答 1

15

目前，我们没有“完全停止”的 PDE 求解器，即您放入 PDE 并运行的求解器。但是，PDE 是通过离散化为 ODE 来求解的，因此为此编写一个完整的 PDE 求解器的方式如下。顺便说一句，这篇博文中对大部分内容进行了更深入的讨论。

带上你的 PDE。现在将运算符离散化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执行了离散化的拉普拉斯算子。

一阶导数项的处理方式类似，但在这种情况下，通常使用逆向方案。您可以使用您的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 教程。现在您只需求解 PDE 的离散 ODE 表示。对于这类方程，如博文中所述，使用Krylov 线性求解器的 Sundials.jlCVODE_BDF方法是一个不错的选择，因此我们这样做：GMRES

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

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

我们将在不久的将来为 PDE 自动执行此操作（任何顺序的任何导数），但它目前正在进行中，因此现在必须进行手动离散化（这在任何数值方法课程中都有讲授，所以还不错！）。希望这可以帮助。如果您需要更多帮助，请查看我们的聊天频道，因为那里的大多数人都会有这种离散化的经验。

于 2018-01-18T00:09:56.843 回答