使用 MATLAB 的 pdepe 的（非）耦合 PDE 的奇怪错误结果，时间加倍

Question

我正在尝试使用 pdpe 求解 1d 中的两个耦合反应扩散方程，即

$\partial_t u_1 = \nabla^2 u_1 + 2k(-u_1^2+u_2)$

$\partial_t u_2 = \nabla^2 u_1 + k(u_1^2-u_2)$

解在 $x\in[0,1]$ 域中，初始条件是两个以 $x=1/2$ 为中心的相同高斯分布。边界条件对两个分量都是吸收的，即 $u_1(0)=u_2(0)=u_1(1)=u_2(1)=0$。

pdepe给了我一个解决方案，没有提示任何错误。但是，我认为解决方案一定是错误的，因为当我将耦合设置为零时，即 $k=0$（并且如果我将其设置得非常小，比如 $k=0.001$），解决方案并不重合用简单扩散方程的解

$\partial_t u = \nabla^2 u$

从 pdepe 本身获得。

奇怪的是，耦合设置为零的 "coupled" 案例中的解 $u_1(t)=u_2(t)$ 以及非耦合案例的解如果我们设置 $t'=2t$，则通过构造 $u(t')$ 重合，也就是说，"coupled" 情况的解的演化速度是非耦合情况的解的两倍。

这是一个最小的工作示例：

偶联案例

    function [xmesh,tspan,sol] = coupled(k) %argument is the coupling k

    std=0.001; %width of initial gaussian
    center=1/2; %center of gaussian

    xmesh=linspace(0,1,10000);
    tspan=linspace(0,1,1000);

    sol = pdepe(0,@pdefun,@icfun,@bcfun,xmesh,tspan);

        function [c,f,s] = pdefun(x,t,u,dudx)
            c=ones(2,1);

            f=zeros(2,1);
            f(1) = dudx(1);
            f(2) = dudx(2);

            s=zeros(2,1);
            s(1) = 2*k*(u(2)-u(1)^2);
            s(2) = k*(u(1)^2-u(2));
        end

        function u0 = icfun(x)
            u0=ones(2,1);
            u0(1) = exp(-(x-center)^2/(2*std^2))/(sqrt(2*pi)*std);
            u0(2) = exp(-(x-center)^2/(2*std^2))/(sqrt(2*pi)*std);
        end

        function [pL,qL,pR,qR] = bcfun(xL,uL,xR,uR,t)
            pL=zeros(2,1);
            pL(1) = uL(1);
            pL(2) = uL(2);

            pR=zeros(2,1);
            pR(1) = uR(1);
            pR(2) = uR(2);

            qL = [0 0;0 0];
            qR = [0 0;0 0];
        end

    end

非耦合案例

    function [xmesh,tspan,sol] = uncoupled()

    std=0.001; %width of initial gaussian
    center=1/2; %center of gaussian

    xmesh=linspace(0,1,10000);
    tspan=linspace(0,1,1000);

    sol = pdepe(0,@pdefun,@icfun,@bcfun,xmesh,tspan);

        function [c,f,s] = pdefun(x,t,u,dudx)
            c=1;

            f = dudx;

            s=0;
        end

        function u0 = icfun(x)
            u0=exp(-(x-center)^2/(2*std^2))/(sqrt(2*pi)*std);
        end

        function [pL,qL,pR,qR] = bcfun(xL,uL,xR,uR,t)
            pL=uL;

            pR=uR;

            qL = 0;
            qR = 0;
        end

    end

现在，假设我们运行

    [xmesh,tspan,soluncoupled] = uncoupled();
    [xmesh,tspan,solcoupled] = coupled(0); %coupling k=0, i.e. uncoupled solutions

可以通过绘制任何时间索引 $it$ 的解来直接检查，即使它们应该相同，每个函数给出的解也不相同，例如

    hold all
    plot(xmesh,soluncoupled(it+1,:),'b')
    plot(xmesh,solcoupled(it+1,:,1),'r')
    plot(xmesh,solcoupled(it+1,:,2),'g')

另一方面，如果我们将非耦合解的时间加倍，则解是相同的

    hold all
    plot(xmesh,soluncoupled(2*it+1,:),'b')
    plot(xmesh,solcoupled(it+1,:,1),'r')
    plot(xmesh,solcoupled(it+1,:,2),'g')

情况$k=0$不是奇异的，可以将$k$设置得小但有限，与情况$k=0$的偏差最小，即解仍然快一倍作为非耦合解决方案。

我真的不明白这是怎么回事。我需要处理耦合情况，但如果 $k\to 0$ 时它没有给出正确的限制，我显然不相信结果。我不明白我可能在哪里犯错。会不会是bug？

Answer 1

我找到了错误的根源。问题在于 coupled() 函数的 bcfun 的 qL 和 qR 变量。 MATLAB 文档，请参阅 here and here，对于 q 应该是矩阵还是列向量，有些含糊不清。我用过矩阵

        qL = [0 0;0 0];
        qR = [0 0;0 0];

但实际上我应该使用列向量

        qL = [0;0];
        qR = [0;0];

令人惊讶的是，pdpe 没有抛出错误，只是给出了错误的结果。这也许应该由开发人员修复。

使用 MATLAB 的 pdepe 的（非）耦合 PDE 的奇怪错误结果，时间加倍

Strange wrong result for (un)coupled PDEs using MATLAB's pdepe, time is doubled

matlab

differential-equations

pde