用MATLAB的pdepe对(Un)耦合PDE的奇异错误结果，时间加倍

Question

我试图在一维中求解两个耦合反应扩散方程，使用pdpe，即

$\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\in0,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^2u$

从私人部门获得的信息。

奇怪的是，当我们设置$t'=2t$时，由耦合集“耦合”到零的解$u_1(t)=u_2(t)$与构造$u(t')$解偶联的解是一致的，即“耦合”情形的解演化速度是解偶联情形解的两倍。

下面是一个最小的工作示例：

耦合案例

    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$时没有给出正确的限制，那么我就不相信结果。我看不出我在哪里会犯错误。会不会是个虫子？

Answer 1

我找到了错误的根源。问题在于耦合()函数的qL和qR变量。MATLAB文档，参见这里和这里，对于q应该是矩阵还是列向量略有含糊。我用过矩阵

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

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

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

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