使用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本身获得。

奇怪的是，将耦合设置为零的“耦合”情况下的解$ u_1（t）= u_2（t）$，如果我们设置，则通过构造$ u（t'）$进行解耦合的情况的解是一致的$ t'= 2t $，也就是说，“耦合”情况的解的发展速度是非耦合情况的解的两倍。

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

耦合的情况

    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 \到0 $时没有给出正确的限制，则显然我不相信结果。我看不出我会在哪里犯错误。可能是错误吗？

[我正在尝试使用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）$解决方案...

Answer 1

我找到了错误的根源。问题出在bcfun的qL和qR变量对coupled（）函数。 MATLAB文档（请参见here和here）在q是矩阵还是列向量方面有些含糊。我使用过矩阵

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