Strange wrong result for (un)coupled equations using pdepe, time is doubled

I am trying to solve two coupled reaction diffusion equations in 1d, using pdpe, namely

The solution is in the domain , with initial conditions being two identical Gaussian profiles centered at . The boundary conditions are absorbing for both components, i.e. .

Pdepe gives me a solution without prompting any errors. However, I think the solutions must be wrong, because when I set the coupling to zero, i.e. (and also if I set it to be very small, say ), the solutions do not coincide with the solution of the simple diffusion equation

as obtained from pdepe itself.

Strangely enough, the solutions from the "coupled" case with coupling set to zero, and the solution for the case uncoupled by construction coincide if we set , that is, the solution of the "coupled" case evolves twice as fast as the solution of the uncoupled case.

Here's a minimal working example:

Coupled case

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=zeros(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

Uncoupled case

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

Now, suppose we run

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

One can directly check by plotting the solutions for any time index that, even if they should be identical, the solutions given by each function are not identical, e.g.

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

On the other hand, if we double the time of the uncoupled solution, the solutions are identical

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

The case is not singular, one can set k to be small but finite, and the deviations from the case are minimal, i.e. the solution still goes twice as fast as the uncoupled solution.

I really don't understand what is going on. I need to work on the coupled case, but obviously I don't trust the results if it does not give the right limit when . I don't see where I could be making a mistake. Could it be a bug?

Answer 1

[m3tr0]

I found the source of the error. The problem lies in the qL and qR variables of bcfun for the coupled() function. The MATLAB documentation, see here and here, is slightly ambiguous on whether the q's should be matrices or column vectors. I had used matrices

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

but in reality I should have used column vectors

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

Amazingly, pdpe didn't throw an error, and simply gave wrong results. This should perhaps be fixed by the developers.