julia - 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 中有更深入的讨论顺便提一句。
带上你的 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是负数。然后当然要结合边界条件。然后你把你的 react 写成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 表示。对于这些类型的方程，如博客文章中所述，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 求解器的容差。
我们将在不久的将来为 PDE 自动执行此操作(任何订单的任何衍生品)，但 it's currently a work-in-progress所以现在必须进行手动离散化(这在任何数值方法类(class)中都有教授，所以还不错!)。希望这可以帮助。如果您需要更多帮助，请check out our chat channel因为那里的大多数人都会有这种离散化的经验。