茱莉亚的抛物型偏微分方程

Question

我试图使用Julia以数字方式求解抛物线偏微分方程,但我找不到任何可以帮助的可访问文档.

这是一个例子: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)

朱莉娅有可能这样做吗？在Matlab中使用pdepe可以解决这类问题.

Answer 1

现在,我们没有"完全停止"的PDE求解器,即你放入PDE的求解器.但是,PDE通过离散化到ODE来解决,因此为此编写一个完整的PDE求解器的方式如下.在BTW 博客文章中更详细地讨论了大部分内容.

拿你的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执行离散化的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

请注意,我没有做到最后,因为我不知道你想要什么边界条件.如果它的Dirichlet为零常数,那么你只需在终点处写下它,但是减去空间的值.而这里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表示.对于这些类型的方程式,如博客文章中所述CVODE_BDF,使用GMRESKrylov线性求解器的Sundials.jl 方法是一个不错的选择,因此我们这样做:

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

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

我们将在不久的将来为PDE自动执行此操作(任何订单的任何衍生产品),但它目前正在进行中,因此现在需要手动离散化(这在任何数值方法课程中都有教授) ,所以它不是太糟糕!).希望这可以帮助.如果您需要更多帮助,请查看我们的聊天频道,因为那里的大多数人都会有这种离散化的经验.